Decompositions and Lapack

Matrix decomposition is a family of methods that aim to represent a matrix as the product of several matrices. Those factors can either allow more efficient operations like inversion or linear system resolution, and might provide some insight regarding intrinsic properties of some data to be analysed (e.g. by observing singular values, eigenvectors, etc.) Some decompositions are implemented in pure Rust or available as bindings to a Fortran Lapack implementation (refer to the section on nalgebra-lapack). In this section, we present decompositions implemented in pure Rust for real or complex matrices:

DecompositionFactorsRust methods
QRQ RQ ~ R.qr()
LU with partial pivotingP1 L UP^{-1} ~ L ~ U.lu()
LU with full pivotingP1 L U Q1P^{-1} ~ L ~ U ~ Q^{-1}.full_piv_lu()
HessenbergP H PP ~ H ~ P^*.hessenberg()
CholeskyL LL ~ L^*.cholesky()
Schur decompositionQ U QQ ~ U ~ Q^*.schur() or .try_schur(eps, max_iter)
Symmetric eigendecompositionU Λ UU ~ \Lambda ~ U^*.hermitian_eigen() or .try_hermitian_eigen(eps, max_iter)
SVDU Σ VU ~ \Sigma ~ V^*.svd(compute_u, compute_v) or .try_svd(compute_u, compute_v, eps, max_iter)

All those methods return a dedicated data structure representing the corresponding decomposition. Those structures themselves often offer interesting features. For example, the LU structure returned by the .lu(...) method allows the efficient resolution of multiple linear systems.

Methods prefixed by .try_ allow the customization of the error epsilon eps and of a hard limit of iteration number max_niter for iterative algorithms. None is returned if the given number of iterations is exceeded before convergence. By default, the relative and absolute error epsilons are equal to the floating-point epsilon (i.e. the difference between 1 and the least value greater than 1 that is representable).

In the following, all .unpack methods consume the decomposition data structure to produce the factors with as little allocations as possible.

Cholesky decomposition#

The Cholesky decomposition of a n × n Hermitian Definite Positive (SDP) matrix MM is composed of a n × n lower-triangular matrix LL such that M=LLM = LL^*. Where LL^* designates the conjugate-transpose of LL. If the input matrix is not SDP, such a decomposition does not exist and the matrix method .cholesky(...) returns None. Note that the input matrix is interpreted as a hermitian matrix. Only its lower-triangular part (including the diagonal) is read by the Cholesky decomposition (its strictly upper-triangular is never accessed in memory). See the wikipedia article for further details about the Cholesky decomposition.

Typical uses of the Cholesky decomposition include the inversion of SDP matrices and resolution of SDP linear systems.

Cholesky decompositpion of a 3x3 matrix.

MethodEffect
.l()Retrieves the lower-triangular factor LL of this decomposition, setting its strictly upper-triangular part to 0.
.l_dirty()Retrieves reference to the lower-triangular factor LL of this decomposition. Its strictly upper-triangular part is not set to 0 and may contain garbage.
.inverse()Computes the inverse of the decomposed matrix.
.solve(b)Solves the system Ax=bAx = b where AA is represented by self and xx the unknown.
.solve_mut(b)Overwrites b with the solution of the system Ax=bAx = b where AA is represented by self and xx the unknown.
.unpack()Consumes self to return the lower-triangular factor LL of this decomposition, setting its strictly upper-triangular part to 0.
.unpack_dirty()Consumes self to return the lower-triangular factor LL of this decomposition. Its strictly upper-triangular part is not set to 0 and may contain garbage.

QR#

The QR decomposition of a general m × n matrix MM is composed of an m × min(n, m) unitary matrix QQ, and a min(n, m) × m upper triangular matrix (or upper trapezoidal if m<nm < n) RR such that M=QRM = QR. This can be used to compute the inverse or a matrix or for solving linear systems. See also the wikipedia article for further details.

QR decomposition of a 3x3 matrix.

MethodEffect
.q()Retrieves the unitary matrix QQ part of the decomposition.
.r()Retrieves the upper-triangular matrix RR part of the decomposition.
.q_tr_mul(rhs)Overwrites rhs with the result of self.q() * rhs. This is much more efficient than computing QQ explicitly.
.is_invertible()Determines if the decomposed matrix is invertible.
.inverse()Computes the inverse of the decomposed matrix.
.solve(b)Solves the system Ax=bAx = b where AA is represented by self and xx the unknown.
.solve_mut(b)Overwrites b with the solution of the system Ax=bAx = b where A is represented by self and xx the unknown.
.unpack()Consumes self to return the two factors (Q,R)(Q, R) of this decomposition.

LU with partial or full pivoting#

The LU decomposition of a general m × n matrix is composed of a m × min(n, m) lower triangular matrix LL with a diagonal filled with 1, and a min(n, m) × m upper triangular matrix UU such that M=LUM = LU. This decomposition is typically used for solving linear systems, compute determinants, matrix inverse, and matrix rank. Two versions of the decomposition are implemented in nalgebra:

  • LU decomposition with partial (row) pivoting which computes additionally only one permutation matrix PP such that PM=LUPM = LU. Implemented only for square matrices.
  • FullPivLU: decomposition with full (row and column) pivoting which computes additionally two permutation matrices PP and QQ such that PMQ=LUPMQ = LU.

Partial pivoting should provide good results in general. Is used internally to compute the determinant, inversibility of a general matrix. Full pivoting is less efficient but more numerically stable. See also the wikipedia article for further details.

LU decomposition of a 3x3 matrix.

MethodEffect
.l()Retrieves the lower-triangular matrix LL part of the decomposition.
.u()Retrieves the upper-triangular matrix UU part of the decomposition.
.p()Computes the explicitly permutation matrix PP that made the decomposition possible.
.is_invertible()Determines if the decomposed matrix is invertible.
.inverse()Computes the inverse of the decomposed matrix.
.determinant()Computes the determinant of the decomposed matrix.
.solve(b)Solves the system Ax=bAx = b where AA is represented by self and xx the unknown.
.solve_mut(b)Overwrites b with the solution of the system Ax=bAx = b where AA is represented by self and xx the unknown.
.unpack()Consumes self to return the three factors (P,L,U)(P, L, U) of this decomposition. The four factors (P,L,U,Q)(P, L, U, Q) are returned when using full pivoting.

Hessenberg decomposition#

The hessenberg decomposition of a square matrix MM is composed of an orthogonal matrix QQ and an upper-Hessenberg matrix HH such that M=QHQM = QHQ^*. The matrix HH being upper-Hessenberg means that all components below its first subdiagonal are zero. See also the wikipedia article for further details.

The hessenberg decomposition is typically used as an intermediate representation of a wide variety of algorithms that can benefit from its structure close to the structure of an upper-triangular matrix.

Hessenberg decomposition of a 3x3 matrix.

MethodEffect
.q()Retrieves the unitary matrix QQ part of the decomposition.
.r()Retrieves the Hessenberg matrix HH of this decomposition.
.unpack()Consumes self to return the two factors (Q,H)(Q, H) of this decomposition.
.unpack_h()Consumes self to return the Hessenberg matrix HH of this decomposition.

Schur decomposition#

The Schur decomposition of a general m × n matrix MM is composed of an m × min(n, m) unitary matrix QQ, and a min(n, m) × m upper quasi-upper-triangular matrix TT such that M=QTQM = QTQ^*. A quasi-upper-triangular matrix is a matrix which is upper-triangular except for some 2x2 blocks on its diagonal (i.e. some of its subdiagonal elements are sometimes non-zero and two consecutive diagonal elements cannot be zero simultaneously).

It is noteworthy that the diagonal elements of the quasi-upper-triangular matrix are the eigenvalues of the decomposed matrix. For real matrices, complex eigenvalues of the 2x2 blocks on the diagonal corresponds to a conjugate pair of complex eigenvalues. In the following example, the decomposed 4x4 real matrix has two real eigenvalues σ1\sigma_1 and σ4\sigma_4 and a conjugate pair of complex eigenvalues σ2\sigma_2 and σ3\sigma_3 equal to the complex eigenvalues of the 2x2 diagonal block in the middle.

Schur decomposition of a 4x4 matrix.

MethodEffect
.eigenvalues()Retrieves the eigenvalues of the decomposed matrix. None if the decomposed matrix is real but some of its eigenvalues should be complex.
.complex_eigenvalues()Retrieves all the eigenvalues of the decomposed matrix returned as complex numbers.
.unpack()Consumes self to return the two factors (Q,T)(Q, T) of this decomposition.

Eigendecomposition of a hermitian matrix#

The eigendecomposition of a square hermitian matrix MM is composed of an unitary matrix QQ and a real diagonal matrix Λ\Lambda such that M=QΛQM = Q\Lambda Q^*. The columns of QQ are called the eigenvectors of MM and the diagonal elements of Λ\Lambda its eigenvalues.

The matrix QQ and the eigenvalues of the decomposed matrix are respectively accessible as public the fields eigenvectors and eigenvalues of the SymmetricEigen structure.

Eigendecomposition of a 3x3 hermitian matrix.

MethodEffect
.recompose()Recomputes the original matrix, i.e., QΛQQ\Lambda{}Q^*. Useful if some eigenvalues or eigenvectors have been manually modified.

Singular Value Decomposition#

The Singular Value Decomposition (SVD) of a rectangular matrix is composed of two orthogonal matrices UU and VV and a diagonal matrix Σ\Sigma with positive (or zero) components. Typical uses of the SVD are the pseudo-inverse, rank computation, and the resolution of least-square problems.

The singular values, left singular vectors, and right singular vectors are respectively stored on the public fields singular_values, u and v_t. Note that v_t represents the adjoint (i.e. conjugate-transpose) of the matrix VV.

SVD decomposition of a 3x3 matrix.

MethodEffect
.recompose()Reconstructs the matrix from its decomposition. Useful if some singular values or singular vectors have been manually modified.
.pseudo_inverse(eps)Computes the pseudo-inverse of the decomposed matrix. All singular values below eps will be interpreted as equal to zero.
.rank(eps)Computes the rank of the decomposed matrix, i.e., the number of singular values strictly greater than eps.
.solve(b, eps)Solves the linear system Ax=bAx = b where AA is the decomposed square matrix and xx the unknown. All singular value smaller or equal to eps are interpreted as zero.

Linear system resolution#

As a simple example the following example demonstrates creation of a general 4x4 matrix AA, a column vector bb, and the resolution of the column vector xx which satisfies the equation Ax=bAx = b.

let a = Matrix4::new(
1.0, 1.0, 2.0, -5.0,
2.0, 5.0, -1.0, -9.0,
2.0, 1.0, -1.0, 3.0,
1.0, 3.0, 2.0, 7.0,
);
let mut b = Vector4::new(3.0, -3.0, -11.0, -5.0);
let decomp = a.lu();
let x = decomp.solve(&b).expect("Linear resolution failed.");
assert_relative_eq!(a * x, b);
/*
* It is possible to perform the resolution in-place.
* This is particularly useful to avoid allocations when
* `b` is a `DVector` or a `DMatrix`.
*/
assert!(decomp.solve_mut(&mut b), "Linear resolution failed.");
assert_relative_eq!(x, b);
/*
* It is possible to solve multiple systems
* simultaneously by using a matrix for `b`.
*/
let b = Matrix4x3::new(
3.0, 2.0, 0.0,
-3.0, 0.0, 0.0,
-11.0, 5.0, -3.0,
-5.0, 10.0, 4.0,
);
let x = decomp.solve(&b).expect("Linear resolution failed.");
assert_relative_eq!(a * x, b);

The example above relies on a LU decomposition of the matrix m. Other decompositions can be used as well, depending on what properties the matrix involved in the linear solve has:

  • The cholesky decomposition will be more efficient if the matrix is known to be hermitian-positive-definite.
  • The SVD decomposition will be useful if the matrix is known to be singular or near-singular. This will allow resolution of systems, ignoring dimensions with near-zero singular values.
  • The LU and QR are good for invertible matrix with no specific properties. Right now, only resolution of invertible, square systems are implemented.
info

If the given b is a matrix then .solve(&b) and .solve_mut(&mut b) will still solve Ax = b, where x is also a matrix. This is equivalent to solving several systems with different right-hand-sides (each being one column of the b matrix).

Lapack integration#

The nalgebra-lapack crate is based on bindings to C or Fortran implementation of the Linear Algebra PACKage, aka. LAPACK. The factorizations supported by nalgebra-lapack are the same as those supported by the pure-Rust version. They are all computed by the constructor of a Rust structure:

DecompositionFactorsRust constructors
QRQ RQ ~ RQR::new(matrix)
LU with partial pivotingP1 L UP^{-1} ~ L ~ ULU::new(matrix)
HessenbergP H PP ~ H ~ P^*Hessenberg::new(matrix)
CholeskyL LL ~ L^*Cholesky::new(matrix)
Schur decompositionQ U QQ ~ U ~ Q^*Schur::new(matrix) or Schur::try_new(matrix)
Symmetric eigendecompositionU Λ UU ~ \Lambda ~ U^*SymmetricEigen::new(matrix) or SymmetricEigen::try_new(matrix)
SVDU Σ VU ~ \Sigma ~ V^*SVD::new(matrix) or SVD::try_new(matrix)

The ::try_new constructors return None if the decomposition fails while ::new constructiors panic.

The next subsections describe how to select the desired Lapack implementation, and provide more details regarding each decomposition.

Setting up nalgebra-lapack#

Several implementations of Lapack exist. The desired one should be selected on your Cargo.toml file by enabling the related feature for your nalgebra-lapack dependency. The currently supported implementations are:

  • OpenBLAS enabled by the openblas feature.
  • netlib enabled by the netlib feature.
  • Accelerate enabled by the accelerate feature on Mac OS only.

The openblas feature is enabled by default. The following examples shows how to enable Accelerate instead:

[dependencies.nalgebra_lapack]
version = "*" # Replace the * by the latest version number.
default-features = false
features = [ "Accelerate" ]

Note that enabling two such features simultaneously will lead to compilation errors inside of the nalgebra-lapack crate itself. Thus, specifying default-features = false is extremely important when selecting an implementation other than the default.