pub struct Cholesky<T> { /* private fields */ }
Expand description
Cholesky decomposition.
Given a square, symmetric positive definite matrix A, there exists an invertible lower triangular matrix L such that
A = L LT.
This is called the Cholesky decomposition of A. For not too ill-conditioned A, the computation of the decomposition is very robust, and it takes about half the effort of an LU decomposition with partial pivoting.
Applications
The Cholesky decomposition can be thought of as a specialized LU decomposition for symmetric positive definite matrices, and so its applications are similar to that of LU.
The following example shows how to compute the Cholesky decomposition of a given matrix. In this example, we also unpack the decomposition to retrieve the L matrix, but in many practical applications we are not so concerned with the factor itself. Instead, we may wish to solve linear systems or compute the determinant or the inverse of a symmetric positive definite matrix. In this case, see the next subsections.
use rulinalg::matrix::decomposition::Cholesky;
// Need to import Decomposition if we want to unpack
use rulinalg::matrix::decomposition::Decomposition;
let x = matrix![ 1.0, 3.0, 1.0;
3.0, 13.0, 11.0;
1.0, 11.0, 21.0 ];
let cholesky = Cholesky::decompose(x)
.expect("Matrix is SPD.");
// Obtain the matrix factor L
let l = cholesky.unpack();
assert_matrix_eq!(l, matrix![1.0, 0.0, 0.0;
3.0, 2.0, 0.0;
1.0, 4.0, 2.0], comp = float);
Solving linear systems
After having decomposed the matrix, one may efficiently solve linear systems for different right-hand sides.
let b1 = vector![ 3.0, 2.0, 1.0];
let b2 = vector![-2.0, 1.0, 0.0];
let y1 = cholesky.solve(b1).expect("Matrix is invertible.");
let y2 = cholesky.solve(b2).expect("Matrix is invertible.");
assert_vector_eq!(y1, vector![ 23.25, -7.75, 3.0 ]);
assert_vector_eq!(y2, vector![-22.25, 7.75, -3.00 ]);
Computing the inverse of a matrix
While computing the inverse explicitly is rarely
the best solution to any given problem, it is sometimes
necessary. In this case, it is easily accessible
through the inverse()
method on Cholesky
.
Computing the determinant of a matrix
As with LU decomposition, the Cholesky
decomposition
exposes a method det
for computing the determinant
of the decomposed matrix. This is a very cheap operation.
Implementations
sourceimpl<T> Cholesky<T>where
T: 'static + Float,
impl<T> Cholesky<T>where
T: 'static + Float,
sourcepub fn decompose(matrix: Matrix<T>) -> Result<Self, Error>
pub fn decompose(matrix: Matrix<T>) -> Result<Self, Error>
Computes the Cholesky decomposition A = L LT for the given square, symmetric positive definite matrix.
Note that the implementation cannot reliably and efficiently verify that the matrix truly is symmetric positive definite matrix, so it is the responsibility of the user to make sure that this is the case. In particular, if the input matrix is not SPD, the returned decomposition may not be a valid decomposition for the input matrix.
Errors
- A diagonal entry is effectively zero to working precision.
- A diagonal entry is negative.
Panics
- The matrix must be square.
sourcepub fn det(&self) -> T
pub fn det(&self) -> T
Computes the determinant of the decomposed matrix.
Note that the determinant of an empty matrix is considered to be equal to 1.
sourcepub fn solve(&self, b: Vector<T>) -> Result<Vector<T>, Error>
pub fn solve(&self, b: Vector<T>) -> Result<Vector<T>, Error>
Solves the linear system Ax = b.
Here A is the decomposed matrix and b is the supplied vector.
Errors
If the matrix is sufficiently ill-conditioned, it is possible that the solution cannot be obtained.
Panics
- The supplied right-hand side vector must be dimensionally compatible with the supplied matrix.