# 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:

Decomposition | Factors | Rust methods |
---|---|---|

QR | $Q ~ R$ | `.qr()` |

LU with partial pivoting | $P^{-1} ~ L ~ U$ | `.lu()` |

LU with full pivoting | $P^{-1} ~ L ~ U ~ Q^{-1}$ | `.full_piv_lu()` |

Hessenberg | $P ~ H ~ P^*$ | `.hessenberg()` |

Cholesky | $L ~ L^*$ | `.cholesky()` |

Schur decomposition | $Q ~ U ~ Q^*$ | `.schur()` or `.try_schur(eps, max_iter)` |

Symmetric eigendecomposition | $U ~ \Lambda ~ U^*$ | `.hermitian_eigen()` or `.try_hermitian_eigen(eps, max_iter)` |

SVD | $U ~ \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 decompositionThe Cholesky decompositon of a `n × n`

Hermitian Definite Positive (SDP) matrix
$M$ is composed of a `n × n`

lower-triangular matrix $L$ such that $M = LL^*$.
Where $L^*$ designates the conjugate-transpose of $L$. 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.

Method | Effect |
---|---|

`.l()` | Retrieves the lower-triangular factor $L$ of this decomposition, setting its strictly upper-triangular part to 0. |

`.l_dirty()` | Retrieves reference to the lower-triangular factor $L$ 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 = b$ where $A$ is represented by `self` and $x$ the unknown. |

`.solve_mut(b)` | Overwrites `b` with the solution of the system $Ax = b$ where $A$ is represented by `self` and $x$ the unknown. |

`.unpack()` | Consumes `self` to return the lower-triangular factor $L$ of this decomposition, setting its strictly upper-triangular part to 0. |

`.unpack_dirty()` | Consumes `self` to return the lower-triangular factor $L$ of this decomposition. Its strictly upper-triangular part is not set to 0 and may contain garbage. |

## #

QRThe QR decomposition of a general `m × n`

matrix $M$ is composed of an ```
m ×
min(n, m)
```

unitary matrix $Q$, and a `min(n, m) × m`

upper triangular matrix
(or upper trapezoidal if $m < n$) $R$ such that $M = 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.

Method | Effect |
---|---|

`.q()` | Retrieves the unitary matrix $Q$ part of the decomposition. |

`.r()` | Retrieves the upper-triangular matrix $R$ part of the decomposition. |

`.q_tr_mul(rhs)` | Overwrites `rhs` with the result of `self.q() * rhs` . This is much more efficient than computing $Q$ explicitly. |

`.is_invertible()` | Determines if the decomposed matrix is invertible. |

`.inverse()` | Computes the inverse of the decomposed matrix. |

`.solve(b)` | Solves the system $Ax = b$ where $A$ is represented by `self` and $x$ the unknown. |

`.solve_mut(b)` | Overwrites `b` with the solution of the system $Ax = b$ where `A` is represented by `self` and $x$ the unknown. |

`.unpack()` | Consumes `self` to return the two factors $(Q, R)$ of this decomposition. |

## #

LU with partial or full pivotingThe LU decomposition of a general `m × n`

matrix is composed of a ```
m × min(n,
m)
```

lower triangular matrix $L$ with a diagonal filled with 1, and a ```
min(n, m)
× m
```

upper triangular matrix $U$ such that $M = 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 $P$ such that $PM = LU$. Implemented only for**square matrices**.`FullPivLU`

: decomposition with full (row and column) pivoting which computes additionally two permutation matrices $P$ and $Q$ such that $PMQ = 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.

Method | Effect |
---|---|

`.l()` | Retrieves the lower-triangular matrix $L$ part of the decomposition. |

`.u()` | Retrieves the upper-triangular matrix $U$ part of the decomposition. |

`.p()` | Computes the explicitly permutation matrix $P$ 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 = b$ where $A$ is represented by `self` and $x$ the unknown. |

`.solve_mut(b)` | Overwrites `b` with the solution of the system $Ax = b$ where $A$ is represented by `self` and $x$ the unknown. |

`.unpack()` | Consumes `self` to return the three factors $(P, L, U)$ of this decomposition. The four factors $(P, L, U, Q)$ are returned when using full pivoting. |

## #

Hessenberg decompositionThe hessenberg decomposition of a square matrix $M$ is composed of an orthogonal matrix $Q$ and an upper-Hessenberg matrix $H$ such that $M = QHQ^*$. The matrix $H$ 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.

Method | Effect |
---|---|

`.q()` | Retrieves the unitary matrix $Q$ part of the decomposition. |

`.r()` | Retrieves the Hessenberg matrix $H$ of this decomposition. |

`.unpack()` | Consumes `self` to return the two factors $(Q, H)$ of this decomposition. |

`.unpack_h()` | Consumes `self` to return the Hessenberg matrix $H$ of this decomposition. |

## #

Schur decompositionThe Schur decomposition of a general `m × n`

matrix $M$ is composed of an
`m × min(n, m)`

unitary matrix $Q$, and a `min(n, m) × m`

upper
quasi-upper-triangular matrix $T$ such that $M = 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 $\sigma_1$ and $\sigma_4$ and a conjugate pair of complex eigenvalues $\sigma_2$ and $\sigma_3$ 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 $(Q, T)$ of this decomposition. |

## #

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

The matrix $Q$ 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., $Q\Lambda{}Q^*$. Useful if some eigenvalues or eigenvectors have been manually modified. |

## #

Singular Value DecompositionThe Singular Value Decomposition (SVD) of a rectangular matrix is composed of two orthogonal matrices $U$ and $V$ 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 $V$.

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 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 = b$ where $A$ is the decomposed square matrix and $x$ the unknown. All singular value smaller or equal to `eps` are interpreted as zero. |

# #

Linear system resolutionAs a simple example the following example demonstrates creation of a general `4x4`

matrix $A$, a column
vector $b$, and the resolution of the column vector $x$ which satisfies the equation $Ax = 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 integrationThe **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:

Decomposition | Factors | Rust constructors |
---|---|---|

QR | $Q ~ R$ | `QR::new(matrix)` |

LU with partial pivoting | $P^{-1} ~ L ~ U$ | `LU::new(matrix)` |

Hessenberg | $P ~ H ~ P^*$ | `Hessenberg::new(matrix)` |

Cholesky | $L ~ L^*$ | `Cholesky::new(matrix)` |

Schur decomposition | $Q ~ U ~ Q^*$ | `Schur::new(matrix)` or `Schur::try_new(matrix)` |

Symmetric eigendecomposition | $U ~ \Lambda ~ U^*$ | `SymmetricEigen::new(matrix)` or `SymmetricEigen::try_new(matrix)` |

SVD | $U ~ \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.

**nalgebra-lapack**#

Setting up 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:

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.