Numerically find the nearest positive semi definite matrix to a symmetric matrix

Question

I have a symmetric matrix $M$ which I want to numerically project onto the positive semi definite cone.

To do so, I decompose it into $M = QDQ^T$ and transform all negative eigenvalues to zero. (according to this post for example How to find the nearest/a near positive definite from a given matrix?)

When I numerically do this (double precision), if M is quite large (say 100*100), the matrix I obtain is not PSD, (according to me, due to numerical imprecision) and I'm obliged to repeat the process a long time to finally get a PSD matrix. My question is : is it a normal side-effect ? If so, is there a trick I missed or a better way to do it ?

Thank you very much for the help.

Answer 1

After you compute $Q$ and $D$, form $D'=\max(D,0)$, and compute $A'=QD'Q^\top$, the algorithms involved in multiplying those matrices do not promise that $A'$ will be exactly $QD'Q^\top$. Most commonly, they are backward stable, and promise that the actual floating-point output will be $(Q+\delta Q)(D'+\delta D')(Q+\delta Q)^\top$, for some small perturbations $\delta Q,\delta D'$. The actual perturbations don't matter here, but they should be small.

So when you replace the diagonal entries of $D$ with zeros, backward stability may still allow roundoff errors to make it look like those zeros were replaced by very small "random" numbers. In effect you can think of it as computing $$Q\max(D, \epsilon\times\mathrm{randn})Q^\top.$$ If you replace the zeros of $D$ with a small number $x>0$ instead, then perturbing $x$ by $\epsilon$ will no longer change its sign and the eigenvalues will come out all positive. This is a common issue with basically any floating point comparison where you compare $u>0$ and $u$ is computed approximately with roundoff errors.