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;endtoc;tic;[L, U] = lu(A);for i = 1:r_rounds    y = L\b;endtoc;tic;for i = 1:r_rounds    c2 = U\y;endtoc;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 getA\bElapsed time is 13.380155 seconds.lu(A)Elapsed time is 0.002662 seconds.y=L\bElapsed time is 10.557989 seconds.c2=U\yElapsed 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.
Why is the speed of the parts of the LU-decomposition so different?
linear solver;sparse;matrix;performance
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: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. 