How to solve $Ax=b$ via backward and forward substitution on Matlab

Question

How can I solve $Ax=b$ in Matlab code via LU factorization. I know that the command [L,U]=LU(A) stores the LU factorization, but how do I use it for solving the system? If I just do $y=L\setminus b$, and then $x=U\setminus y$ does the Matlab backslash computes it efficiently or treats the system as full? Is there a specific command for doing it?

Answer 1

If you write y = L\b;, then MATLAB has no idea what you're original system is. It's just trying to solve $Ly = b$.

Whether doing

[L,U] = LU(A);
y = L\b;
x = U\y;

or

x = A\b;

is more efficient probably depends on A. I'm guessing that A\b is almost always going to be faster, as it may be possible to optimize better than LU factorization.

Answer 2

The MATLAB's backslash operator first determines the structure of the matrix $A$ in order to choose the robust and hopefully most efficient algorithm to solve $Ax=b$. In the documentation of mldivide, you can find the list of recognized matrix structures and what method is used for each of them. So yes, backslash uses the substitution for triangular matrices and can also efficiently solve systems with permuted triangular matrices.

If you know that the LU factorization will be the good method of choice (your $A$ is square, dense, non-symmetric, and non-singular), instead of using x = A \ b, you have the following options:

% Algo 1
[L, U] = lu(A);
y = L \ b;
x = U \ y;

% Algo 2
[L, U, p] = lu(A, 'vector');
y = L \ b(p);
x = U \ y;

% Algo 3
[L, U, p] = lu(A, 'vector');
y = linsolve(L, b(p), struct('LT', true));
x = linsolve(U, y,    struct('UT', true));

Algo 1 is your "basic" LU factorization solver. Backslash determines the structures of $L$ and $U$ and choses the proper substitution algorithm. However, if you do not specify the row permutation as an output argument of lu, the L-factor has permuted rows which, first, complicates the analysis performed by backslash and, second, leads to a less memory-efficient forward solve since we must pass through the rows of $L$ in a permuted order.

Algo 2 is hence an improved version of Algo 1, because $L$ (and $U$) is a truly triangular matrix and hence we have a slight gain in the performance of the forward solve. Of course, we must first permute $b$ accordingly.

Still, the backslash operator must detemine whether or not $L$ and $U$ are triangular and hence must check for zeros in all the appropriate positions. We can avoid this using the linsolve function and inform it directly, which matrix is lower and upper triangular. This gives Algo 3.

Here are some measurements I've made with 10 runs and a random $2000\times 2000$ system:

x = A \ b   18.855 sec
Algo 1      22.739 sec
Algo 2      17.652 sec
Algo 3      16.949 sec