How to do Cholesky Matrix Decomposition with Numpy

Choleski decomposition, also known as Choleski factorization, is one of the commonly used matrix decomposition methods that factorises a Hermitian, positive-definite matrix into the product of a lower triangular matrix and its conjugate transpose. This decomposition is widely used in scientific and engineering applications, such as linear algebra, machine learning, and optimization.

In this blog post, we will learn how to perform Choleski decomposition in Python using the NumPy library. We will also go through some examples to understand the steps and applications of this decomposition.

To start with, let’s import NumPy.

import numpy as np

Cholesky decomposition with Numpy linalg.cholesky()

To show how to use Numpy linalg.cholesky() function to do Cholesky decomposition we will start with an input matrix that is squre and symmetric. Note the off diagonal upper triangular matrix is the same as lower triangular matrix.

A = np.array([[25, 15, -5],
             [15, 18, 0],
             [-5, 0, 11]])
A

array([[25, 15, -5],
       [15, 18,  0],
       [-5,  0, 11]])

We can use symmetric square matrix as argument to linalg.cholesky() to get the decomposition. Numpy’s linalg.cholesky() function returns the factorized lower triangular matrix L as the other matrix is the transpose of it,

# Cholesky factorization in Numpy's linalg
L  = np.linalg.cholesky(A)
L

array([[ 5.,  0.,  0.],
       [ 3.,  3.,  0.],
       [-1.,  1.,  3.]])

Let us see if we can reconstruct the original matrix by using L and its transpose. The transpose of L looks like this

# transpose of matrix L
L.T

array([[ 5.,  3., -1.],
       [ 0.,  3.,  1.],
       [ 0.,  0.,  3.]])

And dot product of L and its transpose can be perform with Numpy’s np.dot() function and the result from the dot product is our original matrix.

# get original matrix from L and L.T
np.dot(L, L.T)

array([[25., 15., -5.],
       [15., 18.,  0.],
       [-5.,  0., 11.]])

We can further verify that by using Numpy’s np.isclose() function to find if each element of the original matrix is close to the reconstructed matrix from L.

# verify of the result is close to original matrix
np.isclose(A, np.dot(L, L.T))

array([[ True,  True,  True],
       [ True,  True,  True],
       [ True,  True,  True]])

Cholesky decomposition on random matrix with Numpy

Let us consider another example of performing Cholesky decomposition in Python using Numpy. This time we will start with a matrix with randomly generated integers. To generate random integers we will use Numpy’s Random generator class with integers() function.

In the example below we have generated a random matrix with integers from 1 to 9 of dimension 4×4.

# gerate a random matrix with integers
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]])

As we can see the above matrix is not symmetric and therefore not positive definite. Let us make it symmetric by doing dot product of the matrix and 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 our matrix is squre and symmetrix. We can apply Numpy’s linalg.cholesky() function to decompose.

L  = np.linalg.cholesky(A)
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]])

We get the lower triangular matrix L as our results.

As before, we can further test the results by reconstructing the original matrix with dot product

np.dot(L, L.T)

array([[ 77.,  71.,  67.,  75.],
       [ 71.,  94.,  95.,  60.],
       [ 67.,  95., 135.,  55.],
       [ 75.,  60.,  55.,  76.]])

And by checking if each reconstructed element is close to the original matrix’s element.

np.isclose(A, np.dot(L, L.T))

array([[ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True]])