Cholesky Decomposition - mikera/core.matrix GitHub Wiki

Cholesky Decomposition is the factorizaiton A = LL^T, where A is an n n symmetric positive matrix, L is an n x n lower triangular matrix with real and positive diagonal entries and L^T is the conjugate transpose of L.
Every symmetric positive definite matrix can be factorized by cholesky decomposition and this cholesky factorization is unique for every such matrix.

A symmetrix n x n real matrix is a positive definite matrix if for every z, where z is a non zero vector of size n, z^TMz > 0. To put more intuitively, if all eigen values of a positive symmetric matrix are real and positive, then it is positive definite.

Just like all positive numbers, all positive symmetric matrices have "square roots". When a positive symmetric matrix is decomposed into LL^T, we may think of L as a matrix square root of A. Matrix square roots of a matrix are not unique.

consider the cholesky decomposition of a 3x3 matrix A,

    ⌈a₁₁ a₁₂ a₁₃⌉
A = |a₂₁ a₂₂ a₂₃|
	⌊a₃₁ a₃₂ a₃₃⌋

    ⌈l₁₁  0   0 ⌉ ⌈l₁₁ l₂₁ l₃₁⌉
  =	|l₂₁ l₂₂  0 | | 0  l₂₂ l₃₂| ≡ LL* 
	⌊l₃₁ l₃₂ l₃₃⌋ ⌊ 0   0  l₃₃⌋

    ⌈l₁₁l₁₁     l₂₁l₁₁             l₃₁l₁₁       ⌉
  =	|l₂₁l₁₁  l₂₁l₂₁+l₂₂l₂₂      l₃₁l₂₁+l₃₂l₂₂   |
	⌊l₃₁l₂₁  l₃₁l₂₁+l₃₂l₂₂  l₃₁l₃₁+l₃₂l₃₂+l₃₃l₃₃⌋  

from here, we can generalize for diagonal elements

and for elements in L below the diagonal

to perform cholesky decomposition in core.matrix, we use the linear/cholesky function. Like all decomposition functions in this library, we can specify which matrix we want out of the resultant matrix by passing a :return parameter in the options map. If :return is not passed, both L and L^T will be returned.

(let [{:keys [L L*]} (cholesky M)] ....)
(let [{:keys [L*]} (cholesky M {:return [:L*]})] ....)

Cholesky decomposition can sometimes be used instead of LU decomposition when positive definite matrices are involved

If we have a to solve an equation Ax = b, where A is positive definite, we factor A as A = LL^T using cholesky decomposition and then solve LL^Tx = b

• forward substitution: Lz = b
• backward substitution: L^Tx = z