How to Compute Cholesky Decomposition of a Matrix?

Cholesky Decomposition

In this post, we will learn about the Cholesky decomposition, definitions of the terms related to this decomposition, the syntax for the function used to implement Cholesky decomposition, and a few examples.

What is Matrix Decomposition?

In linear algebra, matrix decomposition or matrix factorization refers to breaking the matrix down into two or more constituent parts, which, when multiplied together, yield the original matrix.

When the matrix is decomposed into smaller parts, it allows for efficiency and data analysis. Matrix decomposition also optimizes some algorithms, such as least square estimations.


Cholesky Decomposition

Cholesky Decomposition is the type of decomposition of a matrix into its lower triangular matrix and conjugate transpose.

Mathematically, it is defined as A=LL* where A is the original matrix, L is a lower triangular matrix, and L* is its conjugate transpose.

There are three conditions for a matrix to be decomposed by Cholesky.

  • The matrix should be symmetric
  • In the case of complex matrices, it should be Hermitian
  • The elements of the matrix should be positive-definite

Let us define each of these terms in detail.

Square matrix

A matrix is considered a square matrix if the number of rows in the matrix is equal to the number of columns.

One example of such a matrix is given below.

Square Matrix
Square Matrix

Lower Triangular Matrix

A square matrix in which the elements above the principal diagonal are zero is called a lower triangular matrix.

Lower Triangular Matrix
Lower Triangular Matrix

Symmetric Matrix

A square matrix is symmetric if it is equal to its transpose. Interchanging rows and columns obtain a transpose of a matrix.

An example is given below.

Symmetric Matrix
Symmetric Matrix

Conjugate Transpose and Hermitian Matrix

Suppose there is a matrix with complex entries (a+bi) form, then the conjugate of that matrix is obtained by replacing the complex conjugate of all the elements.

Conjugate
Conjugate Matrix

The conjugate transpose of a matrix is obtained by interchanging the rows and columns of the conjugate matrix.

Conjugate Transpose
Conjugate Transpose

Now a Hermitian matrix is a square matrix equal to its conjugate transpose matrix.

Hermitian Matrix
Hermitian Matrix

Positive-Definite Matrix

A symmetric matrix is said to be positive-definite if all the eigen values of the matrix are positive.

The eigenvalues of a matrix can be computed using the np.linalg.eigenvalsh function of the NumPy library.

Let us see an example to check if the given matrix is positive definite.

#to check if the matrix is positive definite
import numpy as np
A = np.array([[2, 1], [1, 2]])
def is_pos_definite(A):
    eig_vals = np.linalg.eigvalsh(A)
    return np.all(eig_vals > 0)
print(is_pos_definite(A))  

import numpy as np: We are importing the NumPy library as np to create arrays and use its functions.

In the next line, we create a 2×2 matrix and store it in a variable called A.

Next, we create a function called is_pos_definite to check if the matrix is positive definite. The matrix A is passed as an array to this function. The def keyword is used to define a new function.

Inside this function, we have a new variable eig_vals, which is used to store the eigenvalues of the matrix A. The np.linalg.eigvalsh is used to compute the eigenvalues of matrix A.

In the next line, we are checking if the eigenvalues of matrix A are all positive(greater than zero) using np.all.

Lastly, we are printing the result of this function. If the eigenvalues of matrix A are positive, this line prints True. Else it prints False.

The output is given below.

Positive Definite
Positive Definite

Applications of Cholesky Decomposition

The main usage of Cholesky’s Decomposition is in Monte Carlo Simulation.

If you are unfamiliar with the Monte Carlo simulation, Check out this article on Monte Carlo in Python.

Apart from higher-level computations, Cholesky Decomposition is also used to solve systems of linear equations, find out matrix inversions, and also perform principal component analysis(PCA).

Check out this article on PCA if you want to know more about data reduction.


Exploring the syntax of linalg.cholesky decomposition

The NumPy implementation of Cholesky decomposition only takes a Symmetric matrix(real-valued) or Hermitian matrix(complex-valued), but in both cases, the matrix should be positive definite.

The np.linalg.cholesky takes a square matrix A and returns the LL* of the matrix.

There is no facility to check if the matrix is Hermitian.

The lower triangular matrix of A is used to compute the Cholesky decomposition.

The syntax of this method is as follows.

linalg.cholesky(a)

There is only one parameter in the syntax.

a: The matrix passed as an argument should be array-like. Hermitian and positive-definite matrices are accepted.

Return Type: Returns a lower triangular matrix.

This function may raise an exception when the matrix is not a positive-definite matrix. This exception is called LinAlgError.

Check out this article on how to solve matrix-related errors-LinAlgError.

Now that we have seen the method’s basic functionality let us dive into the examples!


Computing Cholesky Decomposition for Symmetric Matrix

A symmetric matrix is a square matrix equal to its transpose.

Being symmetric is not the only criterion for a matrix to be decomposed. The matrix should also be positive-definite.

Here is the code for the decomposition of a symmetric matrix which is also positive-definite.

# Example1 2x2 symmetric matrix
import numpy as np
B=np.array([[7,2],[2,1]])
print(B)  

In this code snippet, we just import the NumPy library as np, create a two-dimensional array(2×2 matrix), and store it in the variable B.

print(B) is used to print the matrix residing in B.

Example Matrix 1
Example Matrix 1

Let us find the Cholesky decomposition of this matrix.

L=np.linalg.cholesky(B)
L

We are calling the linalg.cholesky function to decompose matrix B.

This result is stored in L.

The output is a lower triangular matrix, as shown below.

Cholesky Decomposition For Symmetiric Matrix
Cholesky Decomposition For Symmetric Matrix

Computing Cholesky Decomposition for Hermitian matrix

A matrix is said to Hermitian if it equals its own conjugate transpose.

Again, being a Hermitian matrix is not enough for decomposition. The matrix should also be a positive-definite matrix.

Here is the code for the Cholesky decomposition of a Hermitian matrix

#2x2 hermitian matrix
import numpy as np
C=np.array([[2,1-2j],[1+2j,3]])
print(C)

In the first line, we are importing the numpy library.

Next, we create a new variable called C for storing the 2×2 Hermitian matrix with complex numbers.

Next, we print the elements of the matrix to the screen.

The output is given below.

Example Matrix 2
Example Matrix 2

Here is the code for Cholesky’s decomposition.

L=np.linalg.cholesky(C)
L

The decomposed matrix is stored in the variable called L, which is obtained by the method linalg.cholesky.

In the next line, we are printing the elements of the matrix L, which is a lower triangular matrix.

Cholesky Decomposition For 2x2 Hermitian Matrix
Cholesky Decomposition For 2×2 Hermitian Matrix

Computing Cholesky Decomposition for 3×3 Hermitian matrix

Let us do something different in this example. Let us create a hermitian matrix and define a new function to check if the matrix is positive definite with the help of linalg.eigvalsh. Then pass the matrix as an argument to the Cholesky method.

The code to create a 3×3 Hermitian matrix is given below.

#3x3 hermitian matrix
D=np.array([[10, -2 + 3j, 1],
     [-2 - 3j, 6, 2 + 2j],
     [1, 2 - 2j, 7]])
print(D)
def is_pos_definite(D):
    eig_vals = np.linalg.eigvalsh(D)
    return np.all(eig_vals > 0)
print(is_pos_definite(D))  


All we have done in the above snippet is create a three-dimensional array and store it in a variable called D.

Then, we defined a function called is_pos_definite to check if the matrix is positive-definite.

Inside this function, we have also called an instance of linalg.einvalsh to which the newly created matrix is passed as an argument.

The eig_vals is an object used to return the eigenvalues of matrix D.

np.all is a method of the Numpy library that is used to check if all the eigenvalues of the matrix are greater than zero.

Lastly, we are printing the result of the above computation. The last print statement returns True if the eigenvalues are positive else, it returns False.

What do you think the output will be?

Example Matrix 3
Example Matrix 3

Well, the output is true because the matrix passed as an argument is a positive-definite matrix.

Next, we compute the Cholesky Decomposition.

L=np.linalg.cholesky(D)
L

And the output is as follows.

Cholesky Decomposition For 3x3 Hermitian Matrix
Cholesky Decomposition For 3×3 Hermitian Matrix

The output is a bit unclear, so here is the formulated lower triangular matrix.

Lower Triangular Matrix Example 3
Lower Triangular Matrix Example 3

Illustration of LinAlgError Exception

In the previous examples, we have seen the decomposition of positive-definite matrices. But now let us see what happens if we try to use this method on a matrix that is not a positive-definite matrix.

#example of LinAlgError exception
#3x3 hermitian matrix
E=np.array([[1,2-3j,3+4j],[2+3j,0,4-5j],[3-4j,4+5j,2]])
print(E)
def is_pos_definite(E):
    eig_vals = np.linalg.eigvalsh(E)
    return np.all(eig_vals > 0)
print(is_pos_definite(E)) 

In this example, a non-positive-definite Hermitian matrix is stored in a variable called E.

Next, we are using the previously created function is_pos_definite to check if E has positive eigenvalues.

The output is given below.

Example Matrix 4
Example Matrix 4

The output clearly shows that E is not a positive-definite matrix. What will happen if we pass E as an argument to the linalg.cholesky method?

L=np.linalg.cholesky(E)
L
LinAlgError Exception
LinAlgError Exception

As seen in the above image, the linalgerror is raised. This means that the matrix E is not positive-definite.

Conclusion

To summarize everything we have discussed in this post, we first started with what is a matrix decomposition, then Cholesky decomposition, some basic terms related to linear algebra of matrices such as symmetric matrix, transpose of a matrix, Hermitian matrix, Positive-definite matrix, and so on.

We have also seen the conditions a matrix must satisfy to be decomposed by the Cholesky approach, one of the main conditions being the matrix should be positive-definite.

We have also seen the Numpy approach to check if the matrix is a positive-definite matrix.

Next, we have seen the applications of Cholesky decomposition.

We have also seen the syntax of the method, its parameters, and the return object, which is a lower triangular matrix.

Coming to the examples, we have seen the decomposition of a symmetric matrix and also a Hermitian matrix.

Lastly, we have seen how the LinAlgError exception is raised if the matrix we are trying to decompose is not a positive-definite matrix.

References

To know more about this method and examples, refer to the NumPy documentation.

Check out this StackOverflow answer to find out how to eliminate the exception we encountered in the last example.