3
$\begingroup$

I know that an easy way to solve the matrix problem $$A\cdot x=b$$ is the LU decomposition $$\begin{split} L,\,U&=\text{lu}(A)\\ y&=\text{solve}(L,\,b)\\ x&=\text{solve}(U,\,y) \end{split}$$ Thus I implemented the following code in matlab:

r_num = 10000;
r_rounds = 10000;

A = -30*sparse(diag(ones(r_num,1),0))+16*sparse(diag(ones(r_num-1,1),1))+16*sparse(diag(ones(r_num-1,1),-1))-sparse(diag(ones(r_num-2,1),-2))-sparse(diag(ones(r_num-2,1),2));

b = rand(r_num,1);

tic;
for i = 1:r_rounds
    c1 = A\b;
end
toc;

tic;
[L, U] = lu(A);
for i = 1:r_rounds
    y = L\b;

end
toc;
tic;
for i = 1:r_rounds
    c2 = U\y;
end
toc;

I already know that $A\setminus b$ is doing the LU-decomposition in the background, but I still was interested in the result. Now my script returned that the calculation of $L\setminus b$ was nearly as expensive as the direct calculation $A\setminus b$, while the calculation of $U\setminus y$ was extremely fast. But where does the difference come from? Or is my measurement method flawed?

When timing the LU-decomposition separately, I get

A\b
Elapsed time is 13.380155 seconds.
lu(A)
Elapsed time is 0.002662 seconds.
y=L\b
Elapsed time is 10.557989 seconds.
c2=U\y
Elapsed time is 1.229101 seconds.

The reason for timing the LU-decomposition separately (outside of the loop) is because for my later application I will have to deal with a constant $A$, but a large amount of changing $b$-values. Thus I try to keep the amount of elements in the loop as low as possible.

Improve this question
$\endgroup$
  • 1
    $\begingroup$ It can be expected that the forward substitution will be faster than the back substitution, due to this paper. Not sure how that applies to Matlab though. $\endgroup$ – Jakub Klinkovský May 26 '17 at 9:08
  • $\begingroup$ @ChristianClason: I added a timing round with separate timing for the decomposition, it does not change the result (the decomposition is quite fast). The cholesky decomposition is not possible, after my target matrix (what I want to benchmark for) is not symmetric anymore. $\endgroup$ – arc_lupus May 26 '17 at 9:10
  • 1
    $\begingroup$ @jik Doesn't Matlab use a column compressed sparse format, which would reverse that result? That would fit with the idea that it's the U solve which has fewer cache misses and is faster here. $\endgroup$ – origimbo May 26 '17 at 10:24
  • $\begingroup$ @origimbo It apparently doesn't matter here, but anyway: I don't think that using CSC instead of CSR is an improvement for the backward substitution. Although you could process each column from top to bottom, the column-wise traversal completely destroys the per-row aggregated access to the solution vector (or right hand side if you modify it in-place). For some "nice" matrices you might get away with it, but I think using CSR (with the modification proposed in the paper) is generally faster. $\endgroup$ – Jakub Klinkovský May 26 '17 at 16:04
6
$\begingroup$

First, don't forget to also time the LU decomposition in a loop! Otherwise it's not really a fair comparison. If I do that, I get the following timings:

A\b:   Elapsed time is 10.294272 seconds.
lu(A): Elapsed time is 26.675160 seconds.
L\b:   Elapsed time is 9.934748 seconds.
U\y:   Elapsed time is 1.469325 seconds.

You can see that both A\b and lu(A) are of the same order of magnitude, which is expected. Since the matrix in this example is in fact symmetric, you'd also expect that Matlab will not do an LU decomposition. Replacing lu by chol gives a timing of 10.067633 seconds -- very close to Matlab's backslash.

(In fact, you can ask Matlab what it's doing by setting spparms('spumoni',2), which will tell you on A\b:

sp\: bandwidth = 2+1+2.
sp\: is A diagonal? no.
sp\: is band density (1) > bandden (0.5) to try banded solver? yes.
sp\: is LAPACK's banded solver successful? yes.

So it's in fact using a specialized solver.)

To see why L\b is so much slower than U\y, it's instructive to look at the sparsity pattern of both via spy:

sparsity pattern of L sparsity pattern of U

You can see that L has a much more complicated sparsity pattern than U (which is upper tridiagonal, corresponding to the sparsity pattern of A), so the corresponding substitution step will be much slower. In fact, asking Matlab again what it's doing on L\b:

sp\: bandwidth = 2+1+707.
sp\: is A diagonal? no.
sp\: is band density (0.0043793) > bandden (0.5) to try banded solver? no.
sp\: is A triangular? no.
sp\: is A morally triangular? yes.
sp\: permute and solve.

So you see that extra work is being done here.

Improve this answer
$\endgroup$
  • $\begingroup$ Why is the matrix $L$ shown as upper-triangular? Is it in fact $L^T$? Also, assuming that $A$ is symmetric positive-definite, then it has a unique decomposition $A = L L^T$, so the factors $L$ and $U$ should have the same sparsity pattern (up to the transposition). $\endgroup$ – Jakub Klinkovský May 26 '17 at 9:46
  • $\begingroup$ @jlk Good question, but this is what the OP's code produces. (Note that lu, different than chol, produces a permuted lower triangular matrix. The documentation explicitly states Return value L is a product of lower triangular and permutation matrices.) $\endgroup$ – Christian Clason May 26 '17 at 10:04
  • $\begingroup$ Hmm, so a permuted triangular matrix is generally not triangular and L\b is probably doing yet another factorization. Still, I don't see why the pattern of U is nearly diagonal... $\endgroup$ – Jakub Klinkovský May 26 '17 at 10:18
  • $\begingroup$ @jlk Good point! (The pattern of U -- is derived from A, which is pentadiagonal; that's a consequence of how LU decomposition is derived from Gaussian elimination.) $\endgroup$ – Christian Clason May 26 '17 at 10:26
  • $\begingroup$ If you want you can explicitly do the permutation by calling [L,U,P]=lu(A); y=L\(P*b); x=U\y; -- here, L will in fact be lower tridiagonal. $\endgroup$ – Christian Clason May 26 '17 at 10:29

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.