22.3 Linear Algebra on Sparse Matrices

Octave includes a poly-morphic solver for sparse matrices, where the exact solver used to factorize the matrix, depends on the properties of the sparse matrix itself. The cost of determining the matrix type is small relative to the cost of factorizing the matrix itself, but in any case the matrix type is cached once it is calculated, so that it is not re-determined each time it is used in a linear equation.

The selection tree for how the linear equation is solve is

1. If the matrix is not square go to 9.
2. If the matrix is diagonal, solve directly and goto 9
3. If the matrix is a permuted diagonal, solve directly taking into account the permutations. Goto 9
4. If the matrix is banded and if the band density is less than that given by spparms ("bandden") continue, else goto 5.
   1. If the matrix is tridiagonal and the right-hand side is not sparse continue, else goto 4b.
      1. If the matrix is hermitian, with a positive real diagonal, attempt Cholesky factorization using Lapack xPTSV.
      2. If the above failed or the matrix is not hermitian with a positive real diagonal use Gaussian elimination with pivoting using Lapack xGTSV, and goto 9.
   2. If the matrix is hermitian with a positive real diagonal, attempt Cholesky factorization using Lapack xPBTRF.
   3. if the above failed or the matrix is not hermitian with a positive real diagonal use Gaussian elimination with pivoting using Lapack xGBTRF, and goto 9.
5. If the matrix is upper or lower triangular perform a sparse forward or backward subsitution, and goto 9
6. If the matrix is a upper triangular matrix with column permutations or lower triangular matrix with row permutations, perform a sparse forward or backward subsitution, and goto 9
7. If the matrix is hermitian with a real positive diagonal, attempt sparse Cholesky factorization using CHOLMOD.
8. If the sparse Cholesky factorization failed or the matrix is not hermitian with a real positive diagonal, factorize using UMFPACK.
9. If the matrix is not square, or any of the previous solvers flags a singular or near singular matrix, find a minimum norm solution

   FIXME: QR solvers not yet written.

The band density is defined as the number of non-zero values in the matrix divided by the number of non-zero values in the matrix. The banded matrix solvers can be entirely disabled by using spparms to set bandden to 1 (i.e. spparms ("bandden", 1)).

All of the solvers above, expect the banded solvers, calculate an estimate of the condition number. This can be used to detect numerical stability problems in the solution and force a minimum norm solution to be used. However, for narrow banded matrices, the cost of calculating the condition number is significant, and can in fact exceed the cost of factoring the matrix. Therefore the condition number is not calculated for banded matrices, and therefore unless the factorization is exactly singular, these numerical instabilities won't be detected. In cases where, this might be a problem the user is recommended to disable the banded solvers as above, at a significant cost in terms of speed.

The user can force the type of the matrix with the matrix_type function. This overcomes the cost of discovering the type of the matrix. However, it should be noted incorrectly identifying the type of the matrix will lead to unpredictable results, and so matrix_type should be used with care.