Search code examples
matlablinear-algebrasparse-matrixtriangular

Solve *sparse* upper triangular system


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?


Solution

  • 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.