Cholesky decomposition in positive semi-definite matrix

Question

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.

Answer 1

If $\mathbf{A}$ is semi-positive, we have that the sequence $\{\mathbf{A}_k\} = \{\mathbf{A}+\frac{1}{k}\mathbf{I}\}$ consists of positive matrices¹ . 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.

---

¹ proof sketch on wikipedia

Answer 2

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