How to Calculate Cholesky Decomposition in Python
Cholesky decomposition is a widely used matrix factorization method for hermitian positive-definite matrices. It provides many computational benefits for varying algorithms, such as solving linear systems.
We go through how to calculate Cholesky decomposition using the essential scientific computation libraries for Python: NumPy & SciPy. Additionally, we go show you a custom implementation for Cholesky factorization without any external dependencies.
What is Cholesky Decomposition
In simple terms, Cholesky decomposition is matrix factorization for any Hermitian positive-definite matrix \(\mathbf{A}\). It decomposes the matrix into form \[\mathbf{A} = \mathbf{L}\mathbf{L}^\text{H},\]
where \(\mathbf{L}\) is lower triangular matrix and \(\mathbf{L}^\text{H}\) is the conjugate transpose of \(\mathbf{L}\). The diagonal entries of \(\mathbf{L}\) are real and positive.
Before processing, let's clarify a few terms in the Cholesky decomposition definition.
- What is the Hermitian matrix? Matrix \(\mathbf{A}\) is Hermitian if \(\mathbf{A} = \mathbf{A}^\text{H}\).
- What is a positive-definite matrix? A Hermitian matrix \(\mathbf{A}\) is positive-definite if \(\mathbf{x}^\mathbf{H}\mathbf{A}\mathbf{x}\) is positive real value for any non-zero vector \(\mathbf{x}\).
There are a few nice implications of positive-definite matrices:
1. All positive-definite matrices are invertible.
2. The eigenvalues are real and positive.
Basically what we are looking for is a lower triangular matrix \(\mathbf{L}\) that satisfies \(\mathbf{A} = \mathbf{L}\mathbf{L}^\text{H}\) with the aforementioned assumptions. It's not more complex than that.
What are the benefits of Cholesky decomposition?
Most of us are not doing this just for fun 😊. We are of course interested in the benefits of finding the matrix \(\mathbf{L}\).
The main application of Cholesky decomposition is to simplify the solution of linear equations of the form \(\mathbf{A}\mathbf{x} = \mathbf{b}\). Let's assume that we have the matrix factorization available. That is, we know \(\mathbf{L}\) from \(\mathbf{A} = \mathbf{L}\mathbf{L}^\text{H}\). Then, we can solve the linear system in two computationally efficient steps
- Solve \(\mathbf{y}\) from \(\mathbf{L}\mathbf{y} = \mathbf{b}\) by forward substitution.
- Solve \(\mathbf{x}\) from \(\mathbf{L}^\text{H}\mathbf{x} = \mathbf{y}\) by back substitution.
Other applications include
- Matrix inversion: Noticing that \(\mathbf{A}^{-1} = (\mathbf{L}^\text{H})^{-1}(\mathbf{L})^{-1}\), we can see that is enough to invert the lower triangular matrix \(\mathbf{L}\). This is a lot easier to do than inverting a general matrix. Numerical solutions for fixed-size matrices can be found analytically.
- Linear least squares: Not too surprisingly, \(\mathbf{A}\mathbf{x} = \mathbf{b}\) arises in linear least squares problems. As such, we can also utilize Cholesky decomposition to solve it. This should work fine for smaller problems, but for large problems numerically more stable methods should be used.
- Kalman filters: Cholesky decomposition is used as a utility to in tracking the true covariance in Kalman filters. Covariance matrices are inherently Hermitian and positive-definite so the decomposition always exists.
- Non-linear optimization: Many non-linear optimization methods utilize Cholesky decomposition. For example, the quasi-Newton method can use the decomposition on the Hessian matrix to reduce the memory requirements and simplify the search direction computation.
- Monte Carlo simulation: Cholesky decomposition is used in Monte Carlo simulations to generate correlated random variables with given statistics. This is basically done by taking the decomposition of the derived covariance matrix. Then, by multiplying the matrix \(\mathbf{L}\) from the factorization by uncorrelated random samples, the resulting samples are correlated with the given covariance.
Calculating Cholesky Decomposition in Python
Now, we are ready to see how to calculate Cholesky decomposition in practice. We start by creating a common data set to test different methods. Then, we go through matrix factorization using NumPy and SciPy. Finally, we show you how to implement Cholesky–Banachiewicz algorithm in Python to calculate the decomposition.
Test data set
Before going into the details on how to calculate Cholesky decomposition, we will create a common test data set. We'll use the same data for all the methods. This makes it easy to compare the results and verify that everything goes as planned.
We will make a 5x5 test matrix \(\mathbf{A}\).
import random
random.seed(0)
N = 5 # We create a 5x5 test matrix
A = []
for i in range(N):
A.append(([random.random() for j in range(N)]))
# Make it symmetric
for i in range(N):
for j in range(N):
A[i][j] = A[j][i]
# Make it positive semidefinite by making diagonals dominant
for i in range(N):
A[i][i] *= 10This will give us the following test array
[[8.44421852, 0.40493414, 0.90811289, 0.25050634, 0.31014757],
[0.40493414, 7.83798589, 0.50468686, 0.90974626, 0.72983175],
[0.90811289, 0.50468686, 2.81837844, 0.98278548, 0.89883829],
[0.25050634, 0.90974626, 0.98278548, 8.10217236, 0.68398393],
[0.31014757, 0.72983175, 0.89883829, 0.68398393, 4.72142715]]We know that our matrix is Hermitian as it is real and symmetric, but we need to test that our test matrix satisfies the positive-definite condition. For this, we rely on the NumPy library. Namely, we test that all eigenvalues are positive and real.
import numpy as np
D, U = np.linalg.eig(A)
print(np.all(D > 0) and np.all(np.imag(D) == 0))If everything goes according to the plan, this will print out True.
NumPy
NumPy is our first bet, whenever we need to do scientific computation or engineering in Python. It has great support for various linear algebra operations. Also, Cholesky composition is readily supported. It can be accessed via np.linalg.cholesky. Using it is simple, you just pass the matrix you want to factorize and it returns the lower triangular matrix from the decomposition
import numpy as np
L = np.linalg.cholesky(A)This gives us the lower triangular matrix \(\mathbf{L}\)
array([[2.90589375, 0. , 0. , 0. , 0. ],
[0.13934926, 2.79617018, 0. , 0. , 0. ],
[0.31250726, 0.16491815, 1.64119458, 0. , 0. ],
[0.0862063 , 0.32105824, 0.55014625, 2.77290489, 0. ],
[0.10673053, 0.25569221, 0.5016565 , 0.11421472, 2.09283372]])So far so good. Clearly, the matrix is lower triangular and the diagonal values are real and positive. Finally, we can verify that \(\mathbf{L}\mathbf{L}^\text{H} = \mathbf{A}\) by evaluating L @ L.T
array([[8.44421852, 0.40493414, 0.90811289, 0.25050634, 0.31014757],
[0.40493414, 7.83798589, 0.50468686, 0.90974626, 0.72983175],
[0.90811289, 0.50468686, 2.81837844, 0.98278548, 0.89883829],
[0.25050634, 0.90974626, 0.98278548, 8.10217236, 0.68398393],
[0.31014757, 0.72983175, 0.89883829, 0.68398393, 4.72142715]])This matches pretty much exactly the test matrix, which concludes our verification.
SciPy
Whenever something is missing from NumPy, the second place to check is SciPy. SciPy also provides an interface to a good number of linear algebra operations. Cholesky decomposition can be found from scipy.linalg.cholesky. The interface is almost the same as with NumPy with the expectation that by default SciPy implementation returns the upper triangular matrix. To get the lower triangular matrix, you need to explicitly pass lower=True to the method.
import scipy.linalg
L = scipy.linalg.cholesky(A, lower=True)
This gives us the following lower triangular matrix \(\mathbf{L}\)
array([[2.90589375, 0. , 0. , 0. , 0. ],
[0.13934926, 2.79617018, 0. , 0. , 0. ],
[0.31250726, 0.16491815, 1.64119458, 0. , 0. ],
[0.0862063 , 0.32105824, 0.55014625, 2.77290489, 0. ],
[0.10673053, 0.25569221, 0.5016565 , 0.11421472, 2.09283372]])
Check with L @ L.T gives us the test matrix \(\mathbf{A}\) as expect
array([[8.44421852, 0.40493414, 0.90811289, 0.25050634, 0.31014757],
[0.40493414, 7.83798589, 0.50468686, 0.90974626, 0.72983175],
[0.90811289, 0.50468686, 2.81837844, 0.98278548, 0.89883829],
[0.25050634, 0.90974626, 0.98278548, 8.10217236, 0.68398393],
[0.31014757, 0.72983175, 0.89883829, 0.68398393, 4.72142715]])Custom algorithm
If you find yourself in a situation, where you cannot rely on external libraries to calculate the Cholesky decomposition. Don't worry. It is fairly straightforward to calculate the Cholesky decomposition in Python without NumPy or Scipy.
We use Cholesky–Banachiewicz algorithm to calculate the Cholesky decomposition. This is a reasonably simple algorithm with low memory overhead.
import math
L = [ [0] * N for i in range(N) ]
for i in range(N):
for j in range(i+1):
sum = 0
for k in range(j+1):
sum += L[i][k] * L[j][k]
if i == j:
L[i][i] = math.sqrt(A[i][i] - sum)
else:
L[i][j] = (1.0 / L[j][j]) * (A[i][j] - sum)This gives us the following result for \(\mathbf{L}\). We can immediately, see that this clearly matches the results from NumPy and SciPy.
[[2.90589375, 0. , 0. , 0. , 0. ],
[0.13934926, 2.79617018, 0. , 0. , 0. ],
[0.31250726, 0.16491815, 1.64119458, 0. , 0. ],
[0.0862063 , 0.32105824, 0.55014625, 2.77290489, 0. ],
[0.10673053, 0.25569221, 0.5016565 , 0.11421472, 2.09283372]]Summary
Calculating Cholesky decomposition is not hard to do in Python. When you're familiar with the essential libraries, you can easily perform the matrix factorization with a single function call. Even if you are working in an environment where external linear algebra libraries are not accessible, creating your own implementation of the decomposition is not hard.