Cholesky decomposition in positive semi-definite matrix

While trying to apply the algorithm described in the article:

Robust adaptative metropolis algorithm with coerced acceptance rate (2011), Matti Vihola

I used the a Cholesky decomposition to find $S_n$:

$S_n S_n^T = S_{n-1} D_n S_{n-1}^T$

The article ensures that the (known) matrix in the right hand side will always be positive definite. However, after a lot of iterations, it seems that floating point imprecision is causing it to be only positive semi-definite. I would like to know if it is possible (and how?) to get $S_n$ using a robust Cholesky decomposition of a matrix with pivoting which works when standard Cholesky fails.

The problem is when using it, the result is obviously not the $S_n$ matrix, but rather an $X_n$ of the form:

$P^T X_n Y X_n^T P = S_{n-1} D_n S_{n-1}^T$

From what @J.M. has said here, it seems possible:

The Cholesky decomposition might fail in floating point when given a symmetric positive semidefinite matrix. However, one can modify Cholesky to do symmetric pivoting so that the matrix is factored for "as long as the matrix seems positive definite". You'll have to modify your Kalman formula if you adopt this, though.

So, it seems that I need to do some algebra here to get $S_n$, I just don't know how since I'm pretty n00b in such subject.

• The question might get a better response at scicomp.stackexchange.com. – Chris Godsil Jun 13 '13 at 2:05
• I didn't knew this one, looks promising. Even so, I think the case here is more a matter of algebra than a computational problem. – random_user Jun 13 '13 at 2:22

I ran into a similar problem while implementing an algorithm from following paper.

"Online and Batch Learning of Pseudo-Metrics" Shai Shalev-Shwartz, Yoram Singer and Andrew Y. Ng, ICML 2004.

In this case also the algorithm guarantees that the matrix being computed will be positive definite. To me, the issue seems to be computational. So I decided to find the nearest matrix which will allow me to continue the computation. Following paper outlines how this can be done.

"Computing a nearest symmetric positive semidefinite matrix," Nicholas J. Higham, Linear Algebra and its Applications, Volume 103, May 1988, Pages 103-118

1. Symmetric part of $A$ given by $B=(A+A^T)/2$
2. Eigendecomposition of $B$ gives eigenvectors $Z$ and eigen values $\lambda_i$.
3. Positive approximant of $A$ is given as $Z~\text{diag}(\text{max}(0,\lambda_i))Z^T$

Equivalent Matlab code:

 [Z,L] = eig(0.5*(A+A'));
A_out = Z*max(L,0)*Z';


Code available at following website which solves similar problem uses same technique as well.

OASIS: Large scale online Learning of Image Similarity through ranking

• Thanks, do you mean eigenvectors Z on the second item? – random_user May 30 '14 at 2:26
• Also, I tried this approach a while ago, but without the step 1, since in my problem, as far as I remember, there wasn't a problem of symmetry. – random_user May 30 '14 at 2:32
• @random_user Thanks, I corrected the typo. While I came across the OASIS code first, I thought this was some ad-hoc process. But with Higham's paper, I was assured that this was proper thing to do. – Ninad Thakoor May 30 '14 at 18:59

If $\mathbf{A}$ is semi-positive, we have that the sequence $\{\mathbf{A}_k\} = \{\mathbf{A}+\frac{1}{k}\mathbf{I}\}$ consists of positive matrices1 . What's more, $\mathbf{A}_k \to \mathbf{A}$ in operator norm, and with $\mathbf{A}_k = \mathbf{L}_k\mathbf{L}_k^T$, we have that $\mathbf{L}_k$ converges to a limit we name $\mathbf{L}$. Finally, $\mathbf{L}$ has the desired properties, that is $\mathbf{A} = \mathbf{L}\mathbf{L}^T$.

Therefore, an efficient algorithm would be to add a small diagonal $\varepsilon \mathbf{I}$ to $\mathbf{A}$ so that $\mathbf{A} + \varepsilon \mathbf{I}$ is numerically positive definite, and then perform the Cholesky decomposition. We finally remark that computing the eigenvalues of $\mathbf{A}$ asks for more floating point operations.

1 proof sketch on wikipedia