Solve *sparse* upper triangular system

Question

If I want to solve a full upper triangular system, I can call linsolve(A,b,'UT'). However this currently is not supported for sparse matrices. How can I overcome this?

Answer 1

Edit Since what you need is a triangular solve procedure, also called backward/forward substitution, you can use ordinary MATLAB backslash \ operator for that:

x = U\b

As mentioned in the original answer, MATLAB will recognise the fact that your matrix is triangular. To be sure of that, you can compare the performance to cs_usolve procedure found in SuiteSparse. It is a mex function implemented in C that computes sparse triangular solve for upper-triangular sparse matrix (there are similar functions there too: cs_lsolve, cs_utsolve and cs_ltsolve).

You can have a look at a performance comparison of native MATLAB and cs_l(t)solve in the context of sparse Cholesky factorization. Essentially, MATLAB performance is good. The only pitfall is if you want to solve a transposed system

x = U'\b

MATLAB does not recognize that and explicitly creates a transpose of U. In that case you should call cs_utsolve explicitly.

Original answer If your system is symmetric and you only store the upper triangular matrix part (that is how I understood full in your question), and if Cholesky decomposition is suitable for you, chol handles symmetric matrices, if your matrix is positive definite. For indefinite matrices you can use ldl. Both handle sparse storage and work on the symmetric matrix parts.

Newer matlab versions use cholmod and suitesparse for that. That is by far the best performing Cholesky factorization I know of. In matlab it is also parallelised usin parallel BALS.

The factor you obtain from the above functions is upper triangular matrix L such that

A=LL'

All you need to do now is perform forward and backward substitution, which is simple and cheap. In matlab this is automatically done in tha backslash operator

x=L'\(L\b)

the matrix can be sparse, and matlab will recognise that it is upper/lower triangular. You would also use this call together with forward substitution for factors obtained using the cholesky factorization.

Answer 2

UT and LT systems are amongst the easiest systems to solve. Have a read on the wiki about it. Knowing this, it is easy to write your own UT or LT solver:

%# some example data
A = sparse( triu(rand(100)) );
b = rand(100,1);

%# solve UT system by back substitution    
x = zeros(size(b));
for n = size(A,1):-1:1    
    x(n) = (b(n) - A(n,n+1:end)*x(n+1:end) ) / A(n,n);    
end

The procedure is quite similar for LT systems.

Having said that, it is generally much easier and faster to use Matlab's backslash operator:

x = A\b

which also works for spares matrices, as nate already indicated.

Note that this operator also solves UT systems which have non-square A or if A has some elements equal to zero (or < eps) on the main diagonal. It solves these cases in a least-squares sense, which may or may not be desirable for you. You could check for these cases before carrying out the solve:

if size(A,1)==size(A,2) && all(abs(diag(A)) > eps)
    x = A\b;
else
    %# error, warning, whatever you want
end

Read more about the (back)slash operator by typing

>> help \

or

>> help slash

on the Matlab command prompt.

Answer 3

you can use MLDIVIDE( \ ) or MRDIVIDE( / ) operators on your sparse matrices...