In this tutorial, we will learn how to perform Cholesky decomposition in Python using SciPy. We will start with a simple 3×3 matrix with elements hand coded and show how to perform Cholesky decomposition using linalg.cholesky() function in SciPy.
What is Cholesky decompositon
Cholesky decomposition one of the common methods for decomposing a positive-definite matrix into the product of two matrices, a lower triangular matrix and its conjugate transpose.
If the matrix is symmetric and positive definite, then we can decompose it into two matrices as shown below.
A = L * L^T
where L is a lower triangular matrix and L^T is the transpose of L.
A matrix is symmetric if it is equal to its transpose (i.e., if A = A^T). And a matrix is positive definite if all its eigenvalues are positive. If the matrix is not symmetric or positive definite, then it cannot be decomposed using Cholesky decomposition.
Cholesky decomposition is at the heart of many applications including finding efficient numerical solutions. One of the attractive features about Cholesky decomposition is that it is easy to implement and it is numerically stable, which means that it is less likely to produce large errors when used to solve systems of linear equations.
How to perform Cholesky decomposition in Python with SciPy
In Python, we can perform Cholesky decomposition using SciPy’s cholesky() function from the linalg module. The basic syntax of cholesky() function in SciPy is as follows.
scipy.linalg.cholesky(a,
lower=False,
overwrite_a=False,
check_finite=True)
Cholesky() function takes in input matrix a as argument and returns the lower triangular matrix L of Cholesky decomposition.
Note, we can also use Numpy to perform Cholesky decomposition using linalg.cholesky() in Numpy.
Let us see an example of decomposing a 3×3 matrix using SciPy’s cholesky() function.
First, import cholesky from SciPy’s linalg module and Numpy.
from scipy.linalg import cholesky import numpy as np
Here is our matrix to be decomposed using cholesky() function.
A = np.array([[25, 15, -5],
[15, 18, 0],
[-5, 0, 11]])
A
array([[25, 15, -5],
[15, 18, 0],
[-5, 0, 11]])
Note that the matrix is symmetric as we need. Now we can use cholesky() function from SciPy. Here we use the argument lower=TRUE to get lower triangular matrix as output.
cholesky(A, lower=True)
array([[ 5., 0., 0.],
[ 3., 3., 0.],
[-1., 1., 3.]])
By specifying, lower=False we get upper triangular matrix as output.
cholesky(A, lower=False)
array([[ 5., 3., -1.],
[ 0., 3., 1.],
[ 0., 0., 3.]])
Now let us verify if we can reconstruct our input matrix by using the product of L and L^T.
np.dot(L, L.T)
array([[25., 15., -5.],
[15., 18., 0.],
[-5., 0., 11.]])
We can see that the reconstructed matrix matches with our original input matrix.
Cholesky decomposition using random matrix
In the second example showing how to do Cholesky decomposition, we use randomly generated matrix of dimension 4×4 as input and factorize it using cholesky() function SciPy’s linalg module.
First, let generate a matrix sampled from integers using Numpy’s Random Generator object and integers() function.
rng = np.random.default_rng(seed=2022)
a = rng.integers(low=1, high=9, size=(4, 4))
a
array([[6, 2, 6, 1],
[2, 5, 8, 1],
[1, 6, 7, 7],
[7, 1, 5, 1]])
Note our random matrix is not symmetrical. Therefore we can not decompose using Cholesky decomposition as it is.
Let us convert our random matrix into a symmetric matrix by multiplying with its transpose.
# making it symmetric
A = np.dot(a, a.T)
A
array([[ 77, 71, 67, 75],
[ 71, 94, 95, 60],
[ 67, 95, 135, 55],
[ 75, 60, 55, 76]])
Now we have symmetric matrix and this can be used as input to cholesky() function to decompose it.
L = cholesky(A, lower=True)
The lower triangular matrix from cholesky decomposition looks like this.
L
array([[ 8.77496439, 0. , 0. , 0. ],
[ 8.09120093, 5.34157912, 0. , 0. ],
[ 7.63535862, 6.21928056, 6.16618585, 0. ],
[ 8.54704323, -1.71407068, 0.0649585 , 0.07611864]])
Further, we can use L to reconstruct our original matrix.
np.dot(L, L.T)
array([[ 77., 71., 67., 75.],
[ 71., 94., 95., 60.],
[ 67., 95., 135., 55.],
[ 75., 60., 55., 76.]])
If you are interested in using Numpy for Cholesky decomposition, check out this post for using Cholesky decomposition using NumPy’s linalg.cholesky(). Both the posts uses the same dataset for decomposing and as expeted results are identical.
