# How to avoid unnecessary checks when inverting this LU decomposition

Background for the question

I am currently working on a Matlab code in which the systems of linear equations $$Ax_1 = b_1$$, $$Ax_2 = b_2$$, ... have to be solved. As the matrix $$A$$ is constant during the simulation, I LU-decompose it once and for all at the beginning and solve the $$n$$th system as

xn = U\(L\bn);


At this point it is probably worth mentioning that I only know the $$b$$-vectors one at a time, meaning that i cannnot solve the equations all at once.

As can be seen in the figures below, the matrices $$L$$ and $$U$$ very sparse and I therefore store them in sparse format. Moreover, they are upper an lower triangular. By looking at the documentation for the operator "\" it can be seen that Matlab will end up using a triangular solver to solve the systems. Before it gets to that, it however has to see if $$L$$ and $$U$$ are square, compute their bandwidth, see if they are diagonal, check if they look triangular and then finally check that they actually are triangular before using the triangular solver.

The actual question

My question is whether there is a way of bypassing all these checks and go directly to the triangular solver?

What have I tried so far

So far I have looked into the functions "linsolve" and creating upper and lower triangular solvers using "dsp.LowerTriangularSolver" and "dsp.UpperTriangularSolver". Unfortunately, all of these only work on full matrices.

• If your matrix is triangular, why do you compute an LU decomposition at all? If it is lower triangular, then $L=A$ and $U=I$. – Wolfgang Bangerth Sep 20 '19 at 16:11
• The matrix $A$, which I LU-decompose, is not triangular. – Mathias Klahn Sep 20 '19 at 16:17
• So what do the pictures correspond to, then? – Wolfgang Bangerth Sep 23 '19 at 22:01
• The pictures show the sparsity patterns of $L$ and $U$. – Mathias Klahn Sep 24 '19 at 12:29

Looks like using decomposition would simplify your code, and eliminate this problem as a side effect.

dA = decomposition(A);
for k = 1:n
b = compute_b(i);
x(:,k) = dA \ b;
end


In any case, this looks a lot like premature optimization. Why do you think these checks have a large impact on performance?

• First of all thanks for your answer! I don't think that the checks have a large impact on performance, but I do think that they are redundant and therefore wanted to get rid of them. Do you by "premature optimization" mean something like that it is not given that this way of solving the linear system is the fastest, and therefore other strategies should be investigated before optimizing this particular strategy? – Mathias Klahn Sep 23 '19 at 12:40
• "Premature optimization" is a term from computer science; it means spending a lot of time in optimizing a part of code that ends up having a very small impact on the total runtime of the code; see en.wikipedia.org/wiki/Program_optimization#When_to_optimize . In this case, I am not sure that checking if the matrix is triangular is going to give you a big speed-up. Matlab does it this way precisely because usually it's decently cheap. – Federico Poloni Sep 23 '19 at 12:47