CODEWITHSUNDEEP
numpyscipylinear-algebranumerical-methods
Ohad
Ohad

Reputation: 211

Cholesky decomposition of a positive-semidefinite matrix using some Python library

I am looking for a built-in function that implements a Cholesky decomposition of a positive semidefinite matrix in Python.

There exist implementations of NumPy (numpy.linalg.cholesky) and SciPy (scipy.linalg.cholesky), but they work only for positive definite matrices.

I am looking for an extension to positive semidefinite matrices, and didn't find one.

Upvotes: 1

Views: 279

Answers (1)

Matt Haberland
Matt Haberland

Reputation: 3738

You can use scipy.linalg.lapack.dpstrf, which is a pretty-bones wrapper of LAPACK's dpstrf.

import numpy as np
import scipy.linalg

# Generate a positive semidefinite matrix
n, m = 10, 7
A = np.random.rand(n, m)
A = A @ A.T

L, piv, rank, info = scipy.linalg.lapack.dpstrf(A, lower=1)

# For clarity, generate explicit permutation matrix based on definition in [2]
# (For efficiency, we would permute the matrix directly)
P = np.zeros((n, n))
P[piv-1, np.arange(n)] = 1

# We've set `lower=1`; by default it would be upper-triangular
L = np.tril(L)

# Information about complete-pivot Cholesky suggests that
# we should end up with all-zero columns, but that's not
# what we see in LAPACK's output. Apparently this is needed
# if you want to check the decomposition below.
L[:, rank:] = 0

np.testing.assert_equal(rank, m)
np.testing.assert_allclose(P.T @ A @ P, L @ L.T)

Upvotes: 1

Related Questions