Matrix decompositions§
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 nalgebralapack). In this section, we present decompositions implemented in pure Rust for real or complex matrices:
Decomposition  Factors  Rust methods 

QR  .qr() 

LU with partial pivoting  .lu() 

LU with full pivoting  .full_piv_lu() 

Hessenberg  .hessenberg() 

Cholesky  .cholesky() 

Schur decomposition  .schur() or .try_schur(eps, max_iter) 

Symmetric eigendecomposition  .hermitian_eigen() or .try_hermitian_eigen(eps, max_iter) 

SVD  .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 floatingpoint 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 decompositon of a n × n
Symmetric Definite Positive (SDP) matrix
is composed of a n × n
lowertriangular matrix such that .
Where designates the conjugatetranspose of . 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 lowertriangular part (including the diagonal) is
read by the Cholesky decomposition (its strictly uppertriangular 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.
Method  Effect 

.l() 
Retrieves the lowertriangular factor of this decomposition, setting its strictly uppertriangular part to 0. 
.l_dirty() 
Retrieves reference to the lowertriangular factor of this decomposition. Its strictly uppertriangular part is not set to 0 and may contain garbage. 
.inverse() 
Computes the inverse of the decomposed matrix. 
.solve(b) 
Solves the system where is represented by self and the unknown. 
.solve_mut(b) 
Overwrites b with the solution of the system where is represented by self and the unknown. 
.unpack() 
Consumes self to return the lowertriangular factor of this decomposition, setting its strictly uppertriangular part to 0. 
.unpack_dirty() 
Consumes self to return the lowertriangular factor of this decomposition. Its strictly uppertriangular part is not set to 0 and may contain garbage. 
QR§
The QR decomposition of a general m × n
matrix is composed of an m ×
min(n, m)
unitary matrix , and a min(n, m) × m
upper triangular matrix
(or upper trapezoidal if ) such that . This can be used to
compute the inverse or a matrix or for solving linear systems. See also the
wikipedia article for further
details.
Method  Effect 

.q() 
Retrieves the unitary matrix part of the decomposition. 
.r() 
Retrieves the uppertriangular matrix part of the decomposition. 
.q_tr_mul(rhs) 
Overwrites rhs with the result of self.q() * rhs . This is much more efficient than computing explicitly. 
.is_invertible() 
Determines if the decomposed matrix is invertible. 
.inverse() 
Computes the inverse of the decomposed matrix. 
.solve(b) 
Solves the system where is represented by self and the unknown. 
.solve_mut(b) 
Overwrites b with the solution of the system where A is represented by self and the unknown. 
.unpack() 
Consumes self to return the two factors 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 with a diagonal filled with 1, and a min(n, m)
× m
upper triangular matrix such that . 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 such that . Implemented only for square matrices.FullPivLU
: decomposition with full (row and column) pivoting which computes additionally two permutation matrices and such that .
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.
Method  Effect 

.l() 
Retrieves the lowertriangular matrix part of the decomposition. 
.u() 
Retrieves the uppertriangular matrix part of the decomposition. 
.p() 
Computes the explicitly permutation matrix 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 where is represented by self and the unknown. 
.solve_mut(b) 
Overwrites b with the solution of the system where is represented by self and the unknown. 
.unpack() 
Consumes self to return the three factors of this decomposition. The four factors are returned when using full pivoting. 
Hessenberg decomposition§
The hessenberg decomposition of a square matrix is composed of an orthogonal matrix and an upperHessenberg matrix such that . The matrix being upperHessenberg 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 uppertriangular matrix.
Method  Effect 

.q() 
Retrieves the unitary matrix part of the decomposition. 
.r() 
Retrieves the Hessenberg matrix of this decomposition. 
.unpack() 
Consumes self to return the two factors of this decomposition. 
.unpack_h() 
Consumes self to return the Hessenberg matrix of this decomposition. 
Schur decomposition§
The Schur decomposition of a general m × n
matrix is composed of an
m × min(n, m)
unitary matrix , and a min(n, m) × m
upper
quasiuppertriangular matrix such that . A
quasiuppertriangular matrix is a matrix which is uppertriangular except for
some 2x2 blocks on its diagonal (i.e. some of its subdiagonal elements are
sometimes nonzero and two consecutive diagonal elements cannot be zero
simultaneously).
It is noteworthy that the diagonal elements of the quasiuppertriangular 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 and and a conjugate pair of complex eigenvalues and equal to the complex eigenvalues of the 2x2 diagonal block in the middle.
Method  Effect 

.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 of this decomposition. 
Eigendecomposition of a hermitian matrix§
The eigendecomposition of a square hermitian matrix is composed of an unitary matrix and a real diagonal matrix such that . The columns of are called the eigenvectors of and the diagonal elements of its eigenvalues.
The matrix and the eigenvalues of the decomposed matrix are respectively
accessible as public the fields eigenvectors
and eigenvalues
of the
SymmetricEigen
structure.
Method  Effect 

.recompose() 
Recomputes the original matrix, i.e., . 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 and and a diagonal matrix with positive (or zero) components. Typical uses of the SVD are the pseudoinverse, rank computation, and the resolution of leastsquare 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. conjugatetranspose) of the matrix .
Method  Effect 

.recompose() 
Reconstructs the matrix from its decomposition. Useful if some singular values or singular vectors have been manually modified. 
.pseudo_inverse(eps) 
Computes the pseudoinverse 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 where is the decomposed square matrix and 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 , a column
vector , and the resolution of the column vector which satisfies the equation .
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 inplace.
* 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 hermitianpositivedefinite.
 The SVD decomposition will be useful if the matrix is known to be singular or nearsingular. This will allow resolution of systems, ignoring dimensions with nearzero 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.
Note
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 righthandsidesß (each being one column
of the b
matrix).
Lapack integration§
The nalgebralapack crate is based on bindings to C or Fortran implementation of the Linear Algebra PACKage, aka. LAPACK. The factorizations supported by nalgebralapack are the same as those supported by the pureRust version. They are all computed by the constructor of a Rust structure:
Decomposition  Factors  Rust constructors 

QR  QR::new(matrix) 

LU with partial pivoting  LU::new(matrix) 

Hessenberg  Hessenberg::new(matrix) 

Cholesky  Cholesky::new(matrix) 

Schur decomposition  Schur::new(matrix) or Schur::try_new(matrix) 

Symmetric eigendecomposition  SymmetricEigen::new(matrix) or SymmetricEigen::try_new(matrix) 

SVD  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 nalgebralapack§
Several implementations of Lapack exist. The desired one should be selected on
your Cargo.toml
file by enabling the related feature for your
nalgebralapack 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 = "0.5"
defaultfeatures = false
features = [ "Accelerate" ]
Note that enabling two such features simultaneously will lead to compilation
errors inside of the nalgebralapack crate itself. Thus, specifying
defaultfeatures = false
is extremely important when selecting an
implementation other than the default.