numpy - Ensuring positive definite covariance matrix -

the outputs of neural network act entries of covariance matrix. however, 1 one corresponde between outputs , entries results in not positive definite covariance matrices.

thus, read https://www.quora.com/when-carrying-out-the-em-algorithm-how-do-i-ensure-that-the-covariance-matrix-is-positive-definite-at-all-times-avoiding-rounding-issues , https://en.wikipedia.org/wiki/cholesky_decomposition, more specificially "when has real entries, l has real entries , factorization may written a = ll^t".

now outputs corresponds entries of l matrix , generate covariance matrix multiplying transpose.

however, still have error not positive definite matrix. how possible?

i found matrix produces error, see

print l.shape print sigma.shape  s = sigma[1,18,:,:] # matrix gives error l_ = l[1,18,:,:] print l_ s = np.dot(l_,np.transpose(l_)) print s chol = np.linalg.cholesky(s) 

gives output:

(3, 20, 2, 2) (3, 20, 2, 2) [[ -1.69684255e+00   0.00000000e+00]  [ -1.50235415e+00   1.73807144e-04]] [[ 2.87927461  2.54925847]  [ 2.54925847  2.25706792]] ..... linalgerror: matrix not positive definite 

however, code copying values works fine (but not exact same values because not decimals printed)

b = np.array([[-1.69684255e+00, 0.00000000e+00], [-1.50235415e+00, 1.73807144e-04]]) = np.dot(b,b.t) chol_a = np.linalg.cholesky(a) 

so questions are:

• is method of using sigma = ll' correct (with ' transpose)?
• if yes, why getting error? due rounding issues?

edit: computed eigenvalues

print np.linalg.eigvalsh(s) [ -7.89378944432428397703915834426880e-08    5.13634252548217773437500000000000e+00] 

and second case

print np.linalg.eigvalsh(a) [  1.69341869415973178547574207186699e-08    5.13634263409323210680668125860393e+00] 

so there slight negative eigenvalue first case, declares non positive definiteness. how solve this?

this looks numerical issue, in general not true ll' positive definite. example take l matrix each column [1 0 0 0 ... 0] (or more extreme - take l 0 matrix of arbitrary dimensionality), ll' won't pd. in general recommend doing

s = ll' + eps 

which takes care of both problems (for small eps), , 'regularized' covariance estimate. can go "optimal" (under assumtpions) value of eps using ledoit-wolf estimator.