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

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.

• How did you compute the decomposition $QDQ^T$? Did you check if the resulting matrix $Q$ is sufficiently orthogonal, i.e., that $Q^TQ$ is numerically close to the identity matrix?
– cthl
Jan 25, 2019 at 17:28
• I computed it with np.linalg.eigh(). I made a verification and it seems fairly close to the identity : I got a frobenius norm of about $2$x$10^{-14}$ between the two Jan 25, 2019 at 18:24
• Did you try using a small positive number instead of zero for the replaced eigenvalues? Jan 26, 2019 at 3:01
• What's your operational definition of "positive semidefinite"? Jan 26, 2019 at 5:16
• Adding a small positive seems indeed to have solved the issue. Thank you very much ! Jan 26, 2019 at 10:55

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.