GNU Octave - Sparse Linear Algebra

GNU Octave Manual Version 3
by John W. Eaton, David Bateman, Søren Hauberg
Paperback (6"x9"), 568 pages
ISBN 095461206X
RRP £24.95 ($39.95)

Get a printed copy>>>

GNU Octave Manual Version 3

Buy the book here >>>
$39.95

20.2 Linear Algebra on Sparse Matrices

Octave includes a polymorphic solver for sparse matrices, where the exact solver used to factorize the matrix depends on the properties of the sparse matrix itself. Generally, the cost of determining the matrix type is small relative to the cost of factorizing the matrix, and the type is cached once it is calculated, so that it is not re-determined each time the matrix is used.

The selection tree for how the linear equation is solved is:

If the matrix is diagonal, solve directly and goto 8
If the matrix is a permuted diagonal, solve directly taking into account the permutations. Goto 8
If the matrix is square, banded and if the band density is less than that given by spparms ("bandden") continue, else goto 4.
1. If the matrix is tridiagonal and the right-hand side is not sparse continue, else goto 3b.
  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 8.
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 8.
If the matrix is upper or lower triangular perform a sparse forward or backward substitution, and goto 8
If the matrix is a upper triangular matrix with column permutations or lower triangular matrix with row permutations, perform a sparse forward or backward substitution, and goto 8
If the matrix is square, hermitian with a real positive diagonal, attempt sparse Cholesky factorization using cholmod.
If the sparse Cholesky factorization failed or the matrix is not hermitian with a real positive diagonal, and the matrix is square, factorize using umfpack.
If the matrix is not square, or any of the previous solvers flags a singular or near singular matrix, find a minimum norm solution using cxsparse⁽¹⁹⁾.

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

The QR solver factorizes the problem with a Dulmage-Mendhelsohn decomposition, to separate the problem into blocks that can be treated as over-determined, multiple well determined blocks, and a final over-determined block. For matrices with blocks of strongly connected nodes this is a big win as LU decomposition can be used for many blocks. It also significantly improves the chance of finding a solution to over-determined problems rather than just returning a vector of NaN's.

All of the solvers above can 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, triangular or diagonal 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 in these cases, and Octave relies on simpler techniques to detect singular matrices or the underlying lapack code in the case of banded matrices.

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 that identifying the type of the matrix incorrectly will lead to unpredictable results, and so matrix_type should be used with care.

Function File: [n, c] = normest (a, tol)

Estimate the 2-norm of the matrix a using a power series analysis. This is typically used for large matrices, where the cost of calculating the norm (a) is prohibitive and an approximation to the 2-norm is acceptable.

tol is the tolerance to which the 2-norm is calculated. By default tol is 1e-6. c returns the number of iterations needed for normest to converge.

Function File: [est, v] = condest (a, t)

Function File: [est, v] = condest (a, solve, solve_t, t)

Function File: [est, v] = condest (apply, apply_t, solve, solve_t, n, t)

Estimate the 1-norm condition number of a matrix A using t test vectors using a randomized 1-norm estimator. If t exceeds 5, then only 5 test vectors are used.

If the matrix is not explicit, e.g. when estimating the condition number of a given an LU factorization, condest uses the following functions:

apply: A*x for a matrix x of size n by t.
apply_t: A'*x for a matrix x of size n by t.
solve: A \ b for a matrix b of size n by t.
solve_t: A' \ b for a matrix b of size n by t.

The implicit version requires an explicit dimension n.

condest uses a randomized algorithm to approximate the 1-norms.

condest returns the 1-norm condition estimate est and a vector v satisfying norm (A*v, 1) == norm (A, 1) * norm (v, 1) / est. When est is large, v is an approximate null vector.

References:

Nicholas J. Higham and Françoise Tisseur, “A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra.” SIMAX vol 21, no 4, pp 1185-1201.⁽²⁰⁾