SciPy - Cholesky Decomposition

---

---

Cholesky Decomposition in SciPy

Cholesky decomposition is a numerical technique used to decompose a positive-definite matrix into the product of a lower triangular matrix and its transpose. This is particularly useful in numerical computations such as solving systems of linear equations, optimizing algorithms or performing Monte Carlo simulations.

In SciPy the Cholesky decomposition can be performed using the scipy.linalg.cholesky() function. We can define Cholesky decomposition for a lower triangular matrix L, for a given symmeytric positive - definite matrix A as follows −

A = LL^T

Where −

• L is a lower triangular matrix.
• L^T is the transpose of L.

Alternatively we can also compute Cholesky decomposition for the upper triangular matrix U as follows −

A = U^TU

Where −

• U is a Upper triangular matrix.
• U^T is the transpose of U.

Syntax of scipy.linalg.cholesky()

Following is the syntax of using the scipy.linalg.cholesky() function to compute Cholesky Decomposition −

scipy.linalg.cholesky(a, lower=False, overwrite_a=False, check_finite=True)

Where −

• a(array_like): The matrix A to decompose. It must be symmetric and positive-definite.
• lower(bool, optional): If True then computes the lower triangular matrix L and if False(default) then computes the upper triangular matrix U.
• overwrite_a(bool, optional): If True then allows modification of the input matrix to save memory. Default value is False.
• check_finite(bool, optional): If True(default) then checks if the input matrix contains only finite numbers. This may be skipped for performance reasons.

The scipy.linalg.cholesky() function returns the cholesky factor of the input matrix A. If this factor is lower = True then it is lower triangular otherwise it is upper triangular.

Basic Cholesky Decomposition

Following is the example of computing the basic Cholesky Decomposition using the scipy.linalg.cholesky() function in scipy −

import numpy as np
from scipy.linalg import cholesky

# Define a symmetric, positive-definite matrix
A = np.array([[4, 2], 
              [2, 3]])

# Perform Cholesky decomposition to get the lower triangular matrix
L = cholesky(A, lower=True)

# Print the results
print("Matrix A:")
print(A)

print("\nLower Triangular Matrix L:")
print(L)

# Verify: Reconstruct the original matrix using L
A_reconstructed = L @ L.T
print("\nReconstructed Matrix A (L @ L.T):")
print(A_reconstructed)

Here is the output of the basic cholesky decomposition computed using the scipy.linalg.cholesky() function −

Matrix A:
[[4 2]
 [2 3]]

Lower Triangular Matrix L:
[[2.         0.        ]
 [1.         1.41421356]]

Reconstructed Matrix A (L @ L.T):
[[4. 2.]
 [2. 3.]]

Setting lower parameter to True

When the lower parameter in scipy.linalg.cholesky() is set to True then the function computes the lower triangular matrix L. Following is the example of copmuting it −

import numpy as np
from scipy.linalg import cholesky

# Define a symmetric, positive-definite matrix
A = np.array([[6, 3, 4], 
              [3, 6, 5], 
              [4, 5, 10]])

# Perform Cholesky decomposition with lower=True
L = cholesky(A, lower=True)

# Print the results
print("Matrix A:")
print(A)

print("\nLower Triangular Matrix L (lower=True):")
print(L)

# Verify: Reconstruct the original matrix using L
A_reconstructed = L @ L.T
print("\nReconstructed Matrix A (L @ L.T):")
print(A_reconstructed)

Here is the output of the basic cholesky decomposition computed using the scipy.linalg.cholesky() function −

Matrix A:
[[4 2]
 [2 3]]

Lower Triangular Matrix L:
[[2.         0.        ]
 [1.         1.41421356]]

Reconstructed Matrix A (L @ L.T):
[[4. 2.]
 [2. 3.]]
PS D:\Tutorialspoint> python sample.py
Matrix A:
[[ 6  3  4]
 [ 3  6  5]
 [ 4  5 10]]

Lower Triangular Matrix L (lower=True):
[[2.44948974 0.         0.        ]
 [1.22474487 2.12132034 0.        ]
 [1.63299316 1.41421356 2.30940108]]

Reconstructed Matrix A (L @ L.T):
[[ 6.  3.  4.]
 [ 3.  6.  5.]
 [ 4.  5. 10.]]

Key Points

Here are the key points of the Cholesky Decomposition in scipy −

• Matrix Requirements: The matrix A must be symmetric and positive-definite. If A is not positive-definite then the decomposition will fail.
• Efficiency: Cholesky decomposition is faster and more stable than other methods like LU decomposition for positive-definite matrices because it takes advantage of the symmetry and positive-definiteness.
• Error Handling: If the matrix is not positive-definite then a LinAlgError is raised.

Applications

Cholesky decomposition is a powerful numerical tool with applications in various fields particularly in computational mathematics, physics, engineering and machine learning. Below are the key areas where it is commonly used −

• Solving Systems of Linear Equations: For symmetric positive-definite matrices the Cholesky decomposition provides a more efficient way to solve Ax=b compared to methods like LU decomposition.
• Optimization Problems: Many optimization problems involve minimizing quadratic objective functions subject to linear constraints. Cholesky decomposition is used to efficiently solve the systems arising in QP and also used in Machine Learning to optimization algorithms for training models like ridge regression and kernel methods.
• Monte Carlo Simulations: To generate correlated random variables the Cholesky decomposition is used to transform uncorrelated normal random variables.
• Numerical Stability: Cholesky decomposition is more numerically stable than LU decomposition for symmetric positive-definite matrices because it avoids pivoting and exploits matrix symmetry and used in numerical libraries and algorithms to improve the accuracy of results.