Lapack: Cholesky matrix factorization problems

Question

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]])

Answer 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.

Answer 2

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.