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.