Repository navigation
np.linalg.cholesky error on symmetric positive-definite well-conditioned matrix? #8385
Description
I can feed np.linalg.cholesky a symmetric, positive-definite, reasonably conditioned matrix (Sigma_true below) and numpy returns an error message: matrix is not positive definite... Any ideas, folks? I'm a newbie to contributing so helpful suggestions are welcome. If this is indeed expected behavior, then the error message could be more informative. Thanks everyone!
## steps to reproduce issue (Sigma_true does not cholesky decompose)
from __future__ import division
import numpy as np
import scipy.stats
def create_kappa_delta_matrices(kappa_mean, delta_mean, num_samp, random_state1, random_state2):
## normal distribution
assert 1/2*(1 + delta_mean) > kappa_mean
assert kappa_mean> delta_mean
kappas=scipy.stats.norm.rvs(loc= kappa_mean, scale=np.sqrt(kappa_mean*(1-kappa_mean )), size=np.int((num_samp**2 -num_samp)/2 ), random_state=random_state1)
deltas=scipy.stats.norm.rvs(loc= delta_mean, scale=np.sqrt(delta_mean*(1-delta_mean )), size=np.int((num_samp**2 -num_samp)/2 ), random_state=random_state2)
iu1 = np.triu_indices(num_samp, 1) ## indices for upper triangular
Kappa=np.zeros((num_samp,num_samp))
Delta=np.zeros((num_samp,num_samp))
Kappa[iu1]= kappas
Delta[iu1]=deltas
Kappa = Kappa + Kappa.T
Delta = Delta + Delta.T
for i in range(num_samp):
Kappa[i, i]=1
Delta[i, i]=1
Kappa=Kappa.reshape(1, num_samp, num_samp)
Delta=Delta.reshape(1, num_samp, num_samp)
return(Kappa, Delta)
## this is the code that creates the positive-definite well-conditioned matrix (apologies that it is a bit involved):
num_samp=200
kappa_mean=.02
delta_mean= kappa_mean**2 +.001
(Kappa, Delta) = create_kappa_delta_matrices(kappa_mean, delta_mean, num_samp, 1239812, 23425)
Iden = np.diag(np.ones(num_samp)).reshape(1, num_samp, num_samp)
theta = np.concatenate((Kappa, Delta, Iden), axis=0)
sigma2 = np.array((.2, .6, .2)).reshape(3, 1)
Sigma_true = np.sum([theta[i] * sigma2[i] for i in range(theta.shape[0])], axis=0)
## checking that Sigma_true is symmetric positive-definite and well-conditioned:
assert num_samp==np.sum(np.linalg.svd(Sigma_true)[1]>0)
assert (Sigma_true==Sigma_true.T).all()
print(np.linalg.cond(Sigma_true))
## trying to Cholesky-decompose:
chol=np.linalg.cholesky(Sigma_true)## The expected output is a cholesky decomposition of Sigma_true.
## The actual output is:
>>> print(np.linalg.cond(Sigma_true))
43.769107972
>>>chol=np.linalg.cholesky(Sigma_true)
Traceback (most recent call last):
File "<stdin>", line 2, in <module>
File "/home/*****/anaconda3/envs/reml/lib/python3.5/site-packages/numpy/linalg/linalg.py", line 612, in cholesky
r = gufunc(a, signature=signature, extobj=extobj)
File "/home/*****/anaconda3/envs/reml/lib/python3.5/site-packages/numpy/linalg/linalg.py", line 93, in _raise_linalgerror_nonposdef
raise LinAlgError("Matrix is not positive definite")
numpy.linalg.linalg.LinAlgError: Matrix is not positive definiteHere is the version history:
>>> import platform; print(platform.platform())
Linux-3.13.0-92-generic-x86_64-with-debian-jessie-sid
>>> import sys; print("Python", sys.version)
Python 3.5.2 |Continuum Analytics, Inc.| (default, Jul 2 2016, 17:53:06)
[GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
>>> import numpy; print("NumPy", numpy.__version__)
NumPy 1.11.2
Activity
charris commentedon Dec 16, 2016
It's not positive definite.
srmcc commentedon Dec 16, 2016
That's interesting @charris. Unless I missed something silly, the plot thickens. Sigma_true is a real symmetric matrix, therefore it is normal (https://en.wikipedia.org/wiki/Normal_matrix#Special_cases). For normal matrices, the singular values of the SVD are the eigenvalues (https://en.wikipedia.org/wiki/Normal_matrix#Consequences, 2nd proposition). According to the SVD, Sigma_true is positive definite.
I had checked that above with this assertion:
charris commentedon Dec 16, 2016
Singular values are always positive by definition, they are the absolute value of the eigenvalues. All the Hermitean symmetry guarantees is real eigenvalues.
charris commentedon Dec 16, 2016
For instance:
Note that
ais Hermitean with eigenvalues[1, -1], but the singular values are[1, 1].srmcc commentedon Dec 16, 2016
Indeed-- many thanks!