Cholesky decomposition for SPD systems

Cholesky factors a symmetric positive-definite A = L L^T at roughly half the cost of LU, and the factorization succeeding at all (every sqrt argument staying positive) is itself a certificate of positive-definiteness.

import numpy as np

from mathematicskit.linalg import cholesky_decompose, cholesky_solve, is_symmetric_positive_definite, random_spd_matrix

Factor and solve

A = random_spd_matrix(5, seed=1)
result = cholesky_decompose(A)
print("L @ L.T == A:", np.allclose(result.L @ result.L.T, A))

x_true = np.arange(1.0, 6.0)
b = A @ x_true
x = cholesky_solve(result, b)
print("recovered x:", np.round(x, 8))
print("matches x_true:", np.allclose(x, x_true, atol=1e-8))

L @ L.T == A: True
recovered x: [1. 2. 3. 4. 5.]
matches x_true: True

Positive-definiteness as a byproduct of the factorization succeeding

print("is_symmetric_positive_definite(A):", is_symmetric_positive_definite(A))
print("is_symmetric_positive_definite(non-SPD):", is_symmetric_positive_definite(np.array([[1.0, 2.0], [2.0, 1.0]])))

is_symmetric_positive_definite(A): True
is_symmetric_positive_definite(non-SPD): False