2.8. Cholesky Decomposition for Positive Definite Systems¶
Systems of linear equations of the form \(\mathbf{A}\,\bm{x} = \bm{b}\) are often symmetric and positive definite, especially equations derived from physical applications such as electrical, mechanical, and structural systems. We will see in Over-determined Systems and Vector Projections that the solutions to over-determined systems, such as found with least squares regression problems (Least Squares Regression), may make use of positive definite matrices of the form \(\mathbf{A}^T\mathbf{A}\). The diagonal values of a positive definite matrix are always positive. The special symmetry of positive definite matrices allows for a simple and numerically accurate factoring called the Cholesky decomposition that returns an upper triangular matrix, \(\bf{R}\). The transpose of \(\bf{R}\) gives us a lower triangular matrix. The product of the lower and upper triangular matrices is the decomposition of \(\bf{A}\) that we will use as we did with LU decomposition to find the solution using a triangular solver.
The Cholesky decomposition can only be used for positive definite matrices. A positive definite matrix \(\bf{A}\) is real, symmetric, and satisfies the property
The left side of equation (2.15) reduces to a dot product that yields a real, scalar number. Unfortunately, applying the test from equation (2.15) is difficult because it allows for an infinite set of \(\bm{x}\) vectors. An alternate criterion for evaluating whether a matrix is positive definite consists of three tests.
Is the matrix symmetric?
Are all elements along the matrix’s diagonal positive?
Are all of the eigenvalues of the matrix positive?
The procedure for finding the eigenvalues of a matrix is addressed in Application of Eigenvalues and Eigenvectors, but we will not use the third test. A faster final test to see if the Cholesky decomposition can be used is to attempt to use it. The decomposition fails if the algorithm needs to take the square root of a number that is zero or negative. So we can attempt to use the Cholesky decomposition algorithm if the matrix is symmetric and has only positive values along the diagonal.
The development of the Cholesky decomposition algorithm comes from the algebra of how \(\bf{A}\) is computed as the product \(\mathbf{R}^T\, \bf{R}\).
The first row is simple because each element of \(\bf{A}\) comes from only one product.
In the subsequent rows, we find the value of the diagonal element using previously found values in the column of the diagonal element.
The algorithm uses equation (2.17) to check if the matrix might still be positive definite. If the argument to the square root is zero or a negative number, the algorithm is finished because the matrix is not positive definite. The remaining values of each row are found after finding the diagonal element.
The chol_decomp function uses equation (2.16),
equation (2.17), and equation (2.18)
to find the Cholesky decomposition.
It raises an exception if the matrix is not positive definite because it
is not square, symmetric, or has non-positive values on the diagonal.
It returns a boolean flag indicating if the algorithm could complete the
decomposition, which also tells us if the matrix is positive definite or
not.
A matrix may pass the obvious tests for positive definite pass, but may still
not be positive definite.
# -*- coding: utf-8 -*-
"""
File: chol_decomp.py
"""
import numpy as np
from scipy.linalg import issymmetric
def chol_decomp(A: np.array) -> tuple[np.array, np.array]:
# Cholesky decomposition. If A is posive definite, A = R'*R.
# isPositiveDefinite is a boolean variable. Use R only if
# isPositiveDefinite is true.
m,n = A.shape
if m != n:
raise Exception('Matrix is not square')
else:
if not issymmetric(A):
raise Exception('Matrix is not symmetric')
if np.any(np.diag(A <= 0)):
msg = 'Matrix is not positive definite\n'
msg = msg + 'Diagonal not all positive'
raise Exception(msg)
B = A.copy() # A is mutable, don't change it.
isPositiveDefinite = True
for i in np.arange(0,n):
for k in np.arange(0,i): # skip i = 0
B[i,i] = B[i,i] - B[k,i]**2
if A[i,i] <= 0:
isPositiveDefinite = False
R = np.array([])
return R,isPositiveDefinite
B[i,i] = np.sqrt(B[i,i])
for j in np.arange(i+1,n):
for k in np.arange(0,i):
B[i,j] = B[i,j] - B[k,i] * B[k,j]
B[i,j] = B[i,j] / B[i,i]
R = np.triu(B)
return R,isPositiveDefinite
if __name__ == "__main__":
A = np.random.rand(3,3)
try:
R,_ = chol_decomp(A)
except Exception as e:
print(f"Expected error occurred: {e}")
try:
B = A.T @ A
R,isPD = chol_decomp(B)
print(isPD)
print(np.allclose(B, (R.T @ R)))
except Exception as e:
print(f"An error occurred: {e}")
try:
B[1,1] = -1
R,isPD = chol_decomp(B)
print(isPD)
except Exception as e:
print(f"Expected error occurred: {e}")
The Python scipy.linalg module has a few functions that make
use of the Cholesky decomposition. The cholesky function computes
the decomposition. It raises an exception if the decomposition fails.
A triangular solver can be used when the decomposition is successful.
An alternative for when the factor is not needed is cho_factor
and cho_solve. The cho_factor returns the Cholesky factor and
Boolen variable indicating if the factor is lower or upper triangular.
# -*- coding: utf-8 -*-
"""
File: cholesky_test.py
Script to test our Cholesky decomposition and similar functions
in the scipy.linalg module.
"""
import numpy as np
from scipy.linalg import cholesky, cho_factor, cho_solve, solve_triangular
from chol_decomp import chol_decomp
A = np.random.rand(4,4)
# not positive definite case
try:
R,_ = chol_decomp(A)
except Exception as e:
print(f"Expected error occurred: {e}")
try:
R = cholesky(A)
except Exception as e:
print(f"Expected error occurred: {e}")
# positive definite case
b = np.random.rand(4,1)
B = A.T @ A
try:
R = cholesky(B)
R,isPD = chol_decomp(B)
print(isPD)
print(np.allclose(B, (R.T @ R)))
if isPD:
x = solve_triangular(R, solve_triangular(R.T, b, lower=True))
# verify correct
print(np.allclose((B@x - b), np.zeros((4, 1))))
except Exception as e:
print(f"An error occurred: {e}")
# Use cho_factor and cho_solve
try:
R,lower = cho_factor(B)
x = cho_solve((R,lower), b)
# verify correct
print(np.allclose((B@x - b), np.zeros((4, 1))))
except Exception as e:
print(f"An error occurred: {e}")
