Cholesky Decomposition

Cholesky decomposition is to express a symmetric positive definite matrix as a decomposition of the product of a lower triangular matrix and its transpose. It requires that all eigenvalues of the matrix must be positive, so the diagonal elements of the decomposed lower triangle are also positive. The Cholesky decomposition, also known as the square root method, is a modification of the LU decomposition when matrix A is a symmetric positive definite matrix.

(1) LDLT decomposition

If \(A\) is symmetric, we would expect some sort of symmetry in the LU decomposition.

\(U \ne L^T\) due to diagonals in \(U\) and \(L\) not the same. But if we take \(u_{ii}\) out of \(U\) and store in diagonal matrix \(D\) we get \(U = DL^T\), and \(A=LDL^T\) called the LDLT decomposition.

It's roughly half the cost of LU decomposition as we don’t need to compute the upper diagonal part of \(U\).

(2) Cholesky decomposition

When \(A\) is symmetric, positive definite (\(d_{ii} < 0\) in \(D\)), the LDLT decomposition becomes Cholesky decomposition (Cholesky is stable - no pivoting needed): \[A=LD^{1/2}D^{1/2}L^T=GG^T\]

It’s possible to calculate Cholesky directly without forming \(L\) and \(D\). Use \[ A=\left(\begin{array}{ccc} a_{11} & \cdots & a_{1 n} \\ \vdots & \ddots & \vdots \\ a_{n 1} & \cdots & a_{n n} \end{array}\right)=\left(\begin{array}{ccc} g_{11} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ g_{n 1} & \cdots & g_{n n} \end{array}\right)\left(\begin{array}{ccc} g_{11} & \cdots & g_{n 1} \\ \vdots & \ddots & \vdots \\ 0 & \cdots & g_{n n} \end{array}\right) \]

(3) How to determine matrix is pos. def

Definition \(x^{T}Ax > 0\) is not very useful, \(\lambda > 0\) works but it's too expensive to compute eigenvalues.

If not:

1. Check if symmetric

2. Check if all diagonal elements are positive (this is just a sign of pos. def.)

3. Try Cholesky, and if Cholesky fail, exit and perform standard LU

(4) Cholesky decomposition in Python

import numpy as np
from scipy.linalg import ldl
from scipy.linalg import cholesky
from scipy.linalg import solve
A=np.array([[9,3,3],[3,10,7],[3,7,6]]);
# LDLT-decomposition
L, D, P = ldl(A,lower=1)
# Cholesky-decomposition
G = cholesky(A, lower=1)
# Note: scipy.solve can solve systems using
# Cholesky. Number of operations halved.
# assume_a='pos' => ldlt-solution
b = np.array([[8], [-1], [-4]])
x = solve(A,b, assume_a='pos')