Sparse Linear Algebra

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 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 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-Mendelsohn, 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 = normest (A)
— Function File: n = normest (A, tol)
— Function File: [n, c] = normest (...)

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 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, w, iter] = onenormest (A, t)
— Function File: [est, v, w, iter] = onenormest (apply, apply_t, n, t)

Apply Higham and Tisseur's randomized block 1-norm estimator to matrix A using t test vectors. If t exceeds 5, then only 5 test vectors are used.
If the matrix is not explicit, e.g., when estimating the norm of inv (A) given an LU factorization, onenormest applies A and its conjugate transpose through a pair of functions apply and apply_t, respectively, to a dense matrix of size n by t. The implicit version requires an explicit dimension n.
Returns the norm estimate est, two vectors v and w related by norm (w, 1) = est * norm (v, 1), and the number of iterations iter. The number of iterations is limited to 10 and is at least 2.
References:

N.J. Higham and F. Tisseur, A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra. SIMAX vol 21, no 4, pp 1185-1201. http://dx.doi.org/10.1137/S0895479899356080
N.J. Higham and F. Tisseur, A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra. http://citeseer.ist.psu.edu/223007.html

See also: condest, norm, cond.

— Function File: condest (A)
— Function File: condest (A, t)
— Function File: [est, v] = condest (...)
— 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:

N.J. Higham and F. Tisseur, A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra. SIMAX vol 21, no 4, pp 1185-1201. http://dx.doi.org/10.1137/S0895479899356080
N.J. Higham and F. Tisseur, A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra. http://citeseer.ist.psu.edu/223007.html

See also: cond, norm, onenormest.

— Loadable Function: spparms ()
— Loadable Function: vals = spparms ()
— Loadable Function: [keys, vals] = spparms ()
— Loadable Function: val = spparms (key)
— Loadable Function: spparms (vals)
— Loadable Function: spparms ('defaults')
— Loadable Function: spparms ('tight')
— Loadable Function: spparms (key, val)

Query or set the parameters used by the sparse solvers and factorization functions. The first four calls above get information about the current settings, while the others change the current settings. The parameters are stored as pairs of keys and values, where the values are all floats and the keys are one of the following strings:

‘spumoni’
Printing level of debugging information of the solvers (default 0)
‘ths_rel’
Included for compatibility. Not used. (default 1)
‘ths_abs’
Included for compatibility. Not used. (default 1)
‘exact_d’
Included for compatibility. Not used. (default 0)
‘supernd’
Included for compatibility. Not used. (default 3)
‘rreduce’
Included for compatibility. Not used. (default 3)
‘wh_frac’
Included for compatibility. Not used. (default 0.5)
‘autommd’
Flag whether the LU/QR and the '\' and '/' operators will automatically use the sparsity preserving mmd functions (default 1)
‘autoamd’
Flag whether the LU and the '\' and '/' operators will automatically use the sparsity preserving amd functions (default 1)
‘piv_tol’
The pivot tolerance of the umfpack solvers (default 0.1)
‘sym_tol’
The pivot tolerance of the umfpack symmetric solvers (default 0.001)
‘bandden’
The density of non-zero elements in a banded matrix before it is treated by the lapack banded solvers (default 0.5)
‘umfpack’
Flag whether the umfpack or mmd solvers are used for the LU, '\' and '/' operations (default 1)

The value of individual keys can be set with spparms (key, val). The default values can be restored with the special keyword 'defaults'. The special keyword 'tight' can be used to set the mmd solvers to attempt a sparser solution at the potential cost of longer running time.

— Loadable Function: p = sprank (S)

Calculate the structural rank of the sparse matrix S. Note that only the structure of the matrix is used in this calculation based on a Dulmage-Mendelsohn permutation to block triangular form. As such the numerical rank of the matrix S is bounded by sprank (S) >= rank (S). Ignoring floating point errors sprank (S) == rank (S).
See also: dmperm.

— Loadable Function: [count, h, parent, post, r] = symbfact (S)
— Loadable Function: [...] = symbfact (S, typ)
— Loadable Function: [...] = symbfact (S, typ, mode)

Perform a symbolic factorization analysis on the sparse matrix S. Where

S
S is a complex or real sparse matrix.
typ
Is the type of the factorization and can be one of

‘sym’
Factorize S. This is the default.
‘col’
Factorize S' * S.
‘row’
Factorize S * S'.
‘lo’
Factorize S'

mode
The default is to return the Cholesky factorization for r, and if mode is 'L', the conjugate transpose of the Cholesky factorization is returned. The conjugate transpose version is faster and uses less memory, but returns the same values for count, h, parent and post outputs.

The output variables are

count
The row counts of the Cholesky factorization as determined by typ.
h
The height of the elimination tree.
parent
The elimination tree itself.
post
A sparse boolean matrix whose structure is that of the Cholesky factorization as determined by typ.

For non square matrices, the user can also utilize the spaugment function to find a least squares solution to a linear equation.

— Function File: s = spaugment (A, c)

Create the augmented matrix of A. This is given by
          [c * eye(m, m),A; A', zeros(n,
          n)]
This is related to the least squares solution of A \\ b, by
          s * [ r / c; x] = [b, zeros(n,
          columns(b)]
where r is the residual error
          r = b - A * x
As the matrix s is symmetric indefinite it can be factorized with lu, and the minimum norm solution can therefore be found without the need for a qr factorization. As the residual error will be zeros (m, m) for under determined problems, and example can be
          m = 11; n = 10; mn = max(m ,n);
          A = spdiags ([ones(mn,1), 10*ones(mn,1), -ones(mn,1)],
                       [-1, 0, 1], m, n);
          x0 = A \ ones (m,1);
          s = spaugment (A);
          [L, U, P, Q] = lu (s);
          x1 = Q * (U \ (L \ (P  * [ones(m,1); zeros(n,1)])));
          x1 = x1(end - n + 1 : end);
To find the solution of an overdetermined problem needs an estimate of the residual error r and so it is more complex to formulate a minimum norm solution using the spaugment function.
In general the left division operator is more stable and faster than using the spaugment function.

Finally, the function eigs can be used to calculate a limited number of eigenvalues and eigenvectors based on a selection criteria and likewise for svds which calculates a limited number of singular values and vectors.

— Loadable Function: d = eigs (A)
— Loadable Function: d = eigs (A, k)
— Loadable Function: d = eigs (A, k, sigma)
— Loadable Function: d = eigs (A, k, sigma, opts)
— Loadable Function: d = eigs (A, B)
— Loadable Function: d = eigs (A, B, k)
— Loadable Function: d = eigs (A, B, k, sigma)
— Loadable Function: d = eigs (A, B, k, sigma, opts)
— Loadable Function: d = eigs (af, n)
— Loadable Function: d = eigs (af, n, B)
— Loadable Function: d = eigs (af, n, k)
— Loadable Function: d = eigs (af, n, B, k)
— Loadable Function: d = eigs (af, n, k, sigma)
— Loadable Function: d = eigs (af, n, B, k, sigma)
— Loadable Function: d = eigs (af, n, k, sigma, opts)
— Loadable Function: d = eigs (af, n, B, k, sigma, opts)
— Loadable Function: [V, d] = eigs (A, ...)
— Loadable Function: [V, d] = eigs (af, n, ...)
— Loadable Function: [V, d, flag] = eigs (A, ...)
— Loadable Function: [V, d, flag] = eigs (af, n, ...)

Calculate a limited number of eigenvalues and eigenvectors of A, based on a selection criteria. The number of eigenvalues and eigenvectors to calculate is given by k and defaults to 6.
By default, eigs solve the equation where is the corresponding eigenvector. If given the positive definite matrix B then eigs solves the general eigenvalue equation
The argument sigma determines which eigenvalues are returned. sigma can be either a scalar or a string. When sigma is a scalar, the k eigenvalues closest to sigma are returned. If sigma is a string, it must have one of the following values.

'lm'
Largest Magnitude (default).
'sm'
Smallest Magnitude.
'la'
Largest Algebraic (valid only for real symmetric problems).
'sa'
Smallest Algebraic (valid only for real symmetric problems).
'be'
Both Ends, with one more from the high-end if k is odd (valid only for real symmetric problems).
'lr'
Largest Real part (valid only for complex or unsymmetric problems).
'sr'
Smallest Real part (valid only for complex or unsymmetric problems).
'li'
Largest Imaginary part (valid only for complex or unsymmetric problems).
'si'
Smallest Imaginary part (valid only for complex or unsymmetric problems).

If opts is given, it is a structure defining possible options that eigs should use. The fields of the opts structure are:

issym
If af is given, then flags whether the function af defines a symmetric problem. It is ignored if A is given. The default is false.
isreal
If af is given, then flags whether the function af defines a real problem. It is ignored if A is given. The default is true.
tol
Defines the required convergence tolerance, calculated as tol * norm (A). The default is eps.
maxit
The maximum number of iterations. The default is 300.
p
The number of Lanzcos basis vectors to use. More vectors will result in faster convergence, but a greater use of memory. The optimal value of p is problem dependent and should be in the range k to n. The default value is 2 * k.
v0
The starting vector for the algorithm. An initial vector close to the final vector will speed up convergence. The default is for arpack to randomly generate a starting vector. If specified, v0 must be an n-by-1 vector where n = rows (A)
disp
The level of diagnostic printout (0|1|2). If disp is 0 then diagnostics are disabled. The default value is 0.
cholB
Flag if chol (B) is passed rather than B. The default is false.
permB
The permutation vector of the Cholesky factorization of B if cholB is true. That is chol (B(permB, permB)). The default is 1:n.
It is also possible to represent A by a function denoted af. af must be followed by a scalar argument n defining the length of the vector argument accepted by af. af can be a function handle, an inline function, or a string. When af is a string it holds the name of the function to use.
af is a function of the form y = af (x) where the required return value of af is determined by the value of sigma. The four possible forms are

A * x
if sigma is not given or is a string other than 'sm'.
A \ x
if sigma is 0 or 'sm'.
(A - sigma * I) \ x
for the standard eigenvalue problem, where I is the identity matrix of the same size as A.
(A - sigma * B) \ x
for the general eigenvalue problem.

The return arguments of eigs depend on the number of return arguments requested. With a single return argument, a vector d of length k is returned containing the k eigenvalues that have been found. With two return arguments, V is a n-by-k matrix whose columns are the k eigenvectors corresponding to the returned eigenvalues. The eigenvalues themselves are returned in d in the form of a n-by-k matrix, where the elements on the diagonal are the eigenvalues.
Given a third return argument flag, eigs returns the status of the convergence. If flag is 0 then all eigenvalues have converged. Any other value indicates a failure to converge.
This function is based on the arpack package, written by R. Lehoucq, K. Maschhoff, D. Sorensen, and C. Yang. For more information see http://www.caam.rice.edu/software/ARPACK/.
See also: eig, svds.

— Function File: s = svds (A)
— Function File: s = svds (A, k)
— Function File: s = svds (A, k, sigma)
— Function File: s = svds (A, k, sigma, opts)
— Function File: [u, s, v] = svds (...)
— Function File: [u, s, v, flag] = svds (...)

Find a few singular values of the matrix A. The singular values are calculated using
          [m, n] = size(A)
          s = eigs([sparse(m, m), A;
                              A', sparse(n, n)])
The eigenvalues returned by eigs correspond to the singular values of A. The number of singular values to calculate is given by k and defaults to 6.
The argument sigma specifies which singular values to find. When sigma is the string 'L', the default, the largest singular values of A are found. Otherwise, sigma must be a real scalar and the singular values closest to sigma are found. As a corollary, sigma = 0 finds the smallest singular values. Note that for relatively small values of sigma, there is a chance that the requested number of singular values will not be found. In that case sigma should be increased.
opts is a structure defining options that svds will pass to eigs. The possible fields of this structure are documented in eigs. By default, svds sets the following three fields:

tol
The required convergence tolerance for the singular values. The default value is 1e-10. eigs is passed tol / sqrt(2).
maxit
The maximum number of iterations. The default is 300.
disp
The level of diagnostic printout (0|1|2). If disp is 0 then diagnostics are disabled. The default value is 0.

If more than one output is requested then svds will return an approximation of the singular value decomposition of A
          A_approx = u*s*v'
where A_approx is a matrix of size A but only rank k.
flag returns 0 if the algorithm has succesfully converged, and 1 otherwise. The test for convergence is
          norm (A*v - u*s, 1) <= tol * norm (A, 1)
svds is best for finding only a few singular values from a large sparse matrix. Otherwise, svd (full(A)) will likely be more efficient.

22.2 Linear Algebra on Sparse Matrices

Footnotes