Matrix Utilities

driada.utils.matrix.nearestPD(A)[source]

Find the nearest positive-definite matrix to input.

A Python/Numpy port of John D’Errico’s nearestSPD MATLAB code, implementing Higham’s algorithm for computing the nearest symmetric positive semidefinite matrix in the Frobenius norm.

Parameters:: A (array_like) – Input matrix, must be square. Can be non-symmetric.
Returns:: The nearest positive-definite matrix to A.
Return type:: numpy.ndarray

Notes

The algorithm first symmetrizes the input matrix, then uses eigenvalue decomposition and iterative adjustment to ensure positive-definiteness. The function handles numerical precision issues that arise with matrices having eigenvalues near zero.

References

https://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd

N.J. Higham, “Computing a nearest symmetric positive semidefinite matrix” (1988): https://doi.org/10.1016/0024-3795(88)90223-6

Examples

>>> import numpy as np
>>> A = np.array([[1, -1], [-1, 0]])  # Not positive-definite
>>> A_pd = nearestPD(A)
>>> eigvals = np.linalg.eigvals(A_pd)
>>> np.all(eigvals >= 0)  # All eigenvalues are non-negative
True
>>> is_positive_definite(A_pd)  # Result is positive-definite
True

See also

is_positive_definite: Check if a matrix is positive-definite.

driada.utils.matrix.is_positive_definite(B)[source]

Returns true when input is positive-definite, via Cholesky.

Tests whether a matrix is positive-definite by attempting Cholesky decomposition. A matrix is positive-definite if and only if it has a Cholesky decomposition.

Parameters:: B (array_like) – Square matrix to test for positive-definiteness.
Returns:: True if the matrix is positive-definite, False otherwise.
Return type:: bool

Notes

This function uses the Cholesky decomposition as a numerical test for positive-definiteness. The Cholesky decomposition exists if and only if the matrix is symmetric (or Hermitian) and positive-definite.

Examples

>>> import numpy as np
>>> A = np.array([[2, -1], [-1, 2]])  # Positive-definite
>>> is_positive_definite(A)
True

>>> B = np.array([[1, 2], [2, 1]])  # Not positive-definite
>>> is_positive_definite(B)
False