Suppose I want to solve the 2D Poisson equation with Neumann boundary conditions. The solution is non-unique up to an additive constant.

I have previously asked a related question here for the 1D case, which may provide some context for this question: Numerically Solving a Poisson Equation with Neumann Boundary Conditions

There are two problems, which I'll use different notation for:

  • The "Original Equation" $A x = b$, where $A$ is $m \times m$ and has rank $m-1$. This equation is singular because its solution is unique only up to an additive constant, which this equation can not resolve.
  • The "Modified Equation" $C y = d$, where $C$ is $(m+1) \times m$ and has rank $m$. This equation adds uniqueness constraint to the original equation, making $C$ full-rank and the solution unique.

In both cases, $x$ and $y$ should be identical, to machine precision.

This problem can be uniquely solved by specifying a uniqueness constraint. This is done differently for each approach (MATLAB notation):

% Generate A as an m-by-m matrix

% Generate b as an m-by-1 column vector

%% Original Equation
% Solve A*x==b for x
xp = A \ b;   % "Primary" solution for x
              % xp isn't unique, however:
              % The uniqueness constraint must 
              % be applied before x==y.

err_x = norm( A*xp - b, 2 )

% Impose the uniqueness constraint x(4) == 3.14159
x = xp - xp(4) + 3.14159;   % Now, x should equal y (to be calculated)

%% Modified Equation
% Add the constraint x(4) == 3.14159
extraRow = zeros(1,m);
extraRow(4) = 1.0;
C = [A; extraRow];    % Add to the matrix A
d = [b; 3.14159];     % Add to the RHS vector, b

% Solve C*y == d for y
y = C \ d;

err_y = norm( C*y - d, 2 )

I have tried to solve these in MATLAB using the backslash operator (\ or mldivide()) which evaluates the matrix to be solved, then chooses an optimal algorithm to solve it.

In my own tests, MATLAB uses LU decomposition to solve the Original Equation and QR decomposition to solve the Modified Equation.

Test Calculation

I performed the above calculation for an example 2D problem to solve for $\Phi(x,z)$, with a problem size of $Nx=Nz=150$.

2D PDE: $\nabla^2 \Phi = C_1 \frac{\partial f(x)}{\partial{x}}$

Boundary conditions for x-boundaries: $\frac{\partial \Phi}{\partial x} = C_1 f(x)$

Boundary conditions for z-boundaries: $\frac{\partial \Phi}{\partial z} = 0$

Given the form of the source terms, the problem has an analytic, 1D solution for particular $f(x)$: $\frac{\partial \Phi}{\partial x} = C_1 f(x)$, or $\Phi(x,z) = C_1 \int f(x) dx$.


I was surprised to find that the LU decomposition approach yielded far less error than the QR decomposition! (Specifically, err_x $\sim 10^{-11}$ was several orders of magnitude less than err_y $\sim 10^{-8}$.)


For a problem sizes of order $Nx = Nz \sim 100-400$, the modified approach (y = C\d), using QR decomposition, takes roughly twice as long as the original approach (xp = A\b), which uses LU decomposition?

My Question


What's going on for each approach? Is there a compelling reason that LU decomposition out-performs QR decomposition for this type of problem? If not, under what conditions would LU decomposition out-perform QR decomposition, or vice-versa?

(I'm curious how Gaussian Elimination with/without partial pivoting would compare, but that doesn't need to be part of this discussion.)

This question is definitely relevant (but not identical): https://scicomp.stackexchange.com/questions/1026/when-do-orthogonal-transformations-outperform-gaussian-elimination

  • $\begingroup$ How do you quantify 'far less error'? $\endgroup$ – copper.hat Mar 30 '16 at 20:23
  • $\begingroup$ @copper.hat I found that norm(C*y-d,2) is greater than norm(A*x-b,2) by several orders of magnitude for a modest problem size. I'll add this to the code momentarily. $\endgroup$ – jvriesem Mar 30 '16 at 22:03
  • $\begingroup$ @David Thanks for noticing! What is it? (Re-reading the question section, I am not seeing it.) $\endgroup$ – jvriesem Mar 31 '16 at 7:33
  • $\begingroup$ @David Fixed -- thanks! $\endgroup$ – jvriesem Mar 31 '16 at 7:38
  • $\begingroup$ I cleaned up my comments. What are the condition numbers of the matrices. I don't really have much experience with this, but I'm interested in the question. $\endgroup$ – David Mar 31 '16 at 7:39

I recognize that I'm probably far too late for my answer to be of much use to you, but I'll add an answer here for posterity.

First, regarding the choice of method: MATLAB has a very clear methodology by which it selects an particular algorithm to solve this type of equation, as is noted in the documentation. The choice of algorithm is essentially based on how much structure your matrix has that can be exploited to achieve better performance. In your case, the boundary conditions your problem imposes prevent your matrix from being e.g. symmetric, which would allow you to use a faster method. As a result, you end up using LU decomposition, which is one of the slower methods for solving $Ax = b$ when $A$ is square. In the Modified equation case, the imposition of the additional constraint makes your matrix non-square. The majority of common solution methods for systems of linear equations (including LU factorization) do not work for such matrices; in MATLAB, the fallback solution for these types of equations is the QR decomposition. This is why the two different problems you pose end up being solved using different methods.

Second, the speed issue is due precisely to the fact that LU factorization, despite being one of the slower methods available for square matrices, is still faster than QR factorization. The question you linked to actually explains this as well, but what it boils down to is that the number of operations for LU factorization is proportional to $\frac{2}{3}m^3$ (where m is the size of the matrix), whereas for QR factorization the operation count is $~\frac{4}{3}m^3$. These operation counts can be derived by walking through the algorithms for performing these factorizations and seeing how many computations are required for a matrix of size $m$. If you're interested in more detail, you can find more detailed explanations in e.g. Trefethen and Bau's Numerical Linear Algebra.

On the final topic of the error, I'm a little less certain what is going on. The potential culprit that immediately stands out is the fact that in the first case you are applying the constraint after solving, whereas in the second case you do so before solving. This makes it not quite an apples to apples comparison. It would not surprise me if somehow the matrix multiplication involved with solving for y after computing the QR decomposition led to more error propagation. While it is possible that it has to do with the condition numbers of the matrix, I think it's unlikely that this would lead to such a drastic difference between the LU and the QR results. The only reason I think that this argument might have legs is that in QR factorization you need to compute an orthogonal matrix, and numerically it is very easy for that orthogonality to be very slightly violated in ways that might affect the rest of your solution.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.