0

Issue 1

Can somebody recommend a less awkward way of doing a Cholesky factorization in python? Particularly the last line bugs me.

SigmaSqrt = matrix(Sigma)
cvxopt.lapack.potrf(SigmaSqrt)
SigmaSqrt = matrix(np.tril(SigmaSqrt))

Issue 2

I have the problem, that one entire line & column (e.g. all elements in first row and all elements in first column) are zero, lapack fails with the error that the matrix is not positive definite. What is the best way of dealing with this?

Currently I am doing this: (which seems uber-awkward...)

try:
    SigmaSqrt = matrix(Sigma)
    cvxopt.lapack.potrf(SigmaSqrt)
    SigmaSqrt = matrix(np.tril(SigmaSqrt))
except ArithmeticError:
    SigmaSqrt = matrix(Sigma.ix[1:,1:])
    cvxopt.lapack.potrf(SigmaSqrt)
    SigmaSqrt = matrix(np.tril(SigmaSqrt))
    SigmaSqrt = sparse([[v0],[v0[1:].T, SigmaSqrt]])
ARF
  • 7,420
  • 8
  • 45
  • 72

2 Answers2

1

You could just use numpy.linalg.cholesky. Also if all of one column or all of one row are zeros, the matrix will be singular, have at least on eigenvalue that will be zero and therefore, not be positive definite. Since Cholesky is only defined for matrices that are "Hermitian (symmetric if real-valued) and positive-definite" it would not work for it.

EDIT: to "deal with" your problem depends on what you want. Anything you do to make it work would yeild a cholesky that will not be the Cholesky of the original matrix. If you are doing an iterative process and it can be fudge a little, if it is already symmetric then use numpy.linalg.eigvalsh to find the eigenvalues of the matrix. Let d be the most negative eigenvalue. Then set A += (abs(d) + 1e-4) * np.indentity(len(A)). this will make it positive definite.

EDIT: It is a trick used in the Levenberg–Marquardt algorithm. This links to a Wikipedia article on Newton's Method that mentions it as the article on Levenberg–Marquard is doesn't go into this.Also, here is a paper on it as well. Basically this will shift all the eigenvalues by (abs(d) + 1e-4) which will thereby make them all positive, which is a sufficient condition to the matrix being positive definite.

Chris Hagmann
  • 1,086
  • 8
  • 14
  • Thanks for the help. What is the "name" of the procedure you describe in your edit, so I can read brush up my linear algebra? – ARF Apr 16 '14 at 21:06
  • @ArikRaffaelFunke Added more information. It is a trick used in the Levenberg–Marquardt algorithm. – Chris Hagmann Apr 17 '14 at 02:30
  • Also, there is `scipy.linalg.cho_factor` and `scipy.linalg.cho_solve` that can be used in tandem to solve `Ax=b` using the Cholesky. – Chris Hagmann Apr 17 '14 at 14:36
1

Another option is the chompack module: chompack homepage chompack.cholesky I use it myself in combination with the cvxopt module. It works perfectly with the (sparse-) matrices from cvxopt.

Angela
  • 53
  • 4
  • Do you know by any chance whether there is a speed difference for dense matrices between `chompack.cholesky` and `numpy.linalg.cholesky`? Or are you using `chompack` simply because it is convenient with cvxopt matrices? – ARF Oct 07 '14 at 17:37