The Cholesky factorization of a nonsingular symmetric matrix $A$ is $$A = LDL^T,$$ where $L$ is a lower triangular matrix with unit diagonal, and $D$ is a diagonal matrix.
In?[1]:
import numpy as np
In?[2]:
def cholesky(A):
# Find the dimension of the matrix
n = np.shape(A)[0]
# Initialize empty matrices L and D
L, D = np.zeros((n, n)), np.zeros((n, n))
# Initialize A_0 = A
A_k = A
for k in range(0, n):
# Extract colums of L and diagonal elements of D
L[:, k] = A_k[:, k] / A_k[k, k]
D[k, k] = A_k[k, k]
# Compute the next iteration of A_k
A_k = A_k - D[k, k] * (L[:, k : k + 1] @ L[:, k : k + 1].T)
return L, D
Let us first test the algorithm.
In?[3]:
N = 10
a = np.random.randint(-100, 100, size=(N, N))
A = (a + a.T) / 2
L, D = cholesky(A)
print(np.linalg.norm((L @ D @ L.T - A), ord=2))
3.4245319277847024e-14
One could add code to make sure that $A$ is indeed a symmetric $n \times n$ nonsingular matrix.
In?[4]:
def cholesky(A):
n = np.shape(A)[0]
if n != np.shape(A)[1]:
raise Exception("Matrix is not square")
if not np.all(np.abs(A - A.T) < 1e-10):
raise Exception("Matrix is not symmetric")
L, D = np.zeros((n, n)), np.zeros((n, n))
A_k = A
for k in range(0, n):
if np.abs(A_k[k, k]) <= 1e-10:
raise Exception("Matrix is singular")
L[:, k] = A_k[:, k] / A_k[k, k]
D[k, k] = A_k[k, k]
A_k = A_k - D[k, k] * (L[:, k : k + 1] @ L[:, k : k + 1].T)
return L, D
We can use our algorithm to compute the Cholesky factorization of the matrix from exercise 2.7.
In?[5]:
A = np.array(
[
[1, 1, 0, 0, 0, 0],
[1, 2, 1, 0, 0, 0],
[0, 1, 2, 1, 0, 0],
[0, 0, 1, 3, 1, 0],
[0, 0, 0, 1, 3, 1],
[0, 0, 0, 0, 1, np.pi],
]
)
L, D = cholesky(A)
In?[6]:
print(L)
[[1. 0. 0. 0. 0. 0. ] [1. 1. 0. 0. 0. 0. ] [0. 1. 1. 0. 0. 0. ] [0. 0. 1. 1. 0. 0. ] [0. 0. 0. 0.5 1. 0. ] [0. 0. 0. 0. 0.4 1. ]]
In?[7]:
print(D)
[[1. 0. 0. 0. 0. 0. ] [0. 1. 0. 0. 0. 0. ] [0. 0. 1. 0. 0. 0. ] [0. 0. 0. 2. 0. 0. ] [0. 0. 0. 0. 2.5 0. ] [0. 0. 0. 0. 0. 2.74159265]]
In?[8]:
print(np.pi - (2/5))
2.741592653589793
In?[9]:
A = np.array(
[
[1, 1, 0, 0, 0, 0],
[1, 2, 1, 0, 0, 0],
[0, 1, 2, 1, 0, 0],
[0, 0, 1, 3, 1, 0],
[0, 0, 0, 1, 3, 1],
[0, 0, 0, 0, 1, (2/5)],
]
)
L, D = cholesky(A)
--------------------------------------------------------------------------- Exception Traceback (most recent call last) Cell In[9], line 12 1 A = np.array( 2 [ 3 [1, 1, 0, 0, 0, 0], (...) 9 ] 10 ) ---> 12 L, D = cholesky(A) Cell In[4], line 15, in cholesky(A) 13 for k in range(0, n): 14 if np.abs(A_k[k, k]) <= 1e-10: ---> 15 raise Exception("Matrix is singular") 17 L[:, k] = A_k[:, k] / A_k[k, k] 18 D[k, k] = A_k[k, k] Exception: Matrix is singular
In?[?]: