1
$\begingroup$

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.

The sparsity pattern of U] The sparsity pattern of L

$\endgroup$
  • $\begingroup$ 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$. $\endgroup$ – Wolfgang Bangerth Sep 20 at 16:11
  • $\begingroup$ The matrix $A$, which I LU-decompose, is not triangular. $\endgroup$ – Mathias Klahn Sep 20 at 16:17
  • $\begingroup$ So what do the pictures correspond to, then? $\endgroup$ – Wolfgang Bangerth Sep 23 at 22:01
  • $\begingroup$ The pictures show the sparsity patterns of $L$ and $U$. $\endgroup$ – Mathias Klahn Sep 24 at 12:29
1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Mathias Klahn Sep 23 at 12:40
  • $\begingroup$ "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. $\endgroup$ – Federico Poloni Sep 23 at 12:47

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.