Converting a matrix to the nearest positive definite matrix

Question

I have a matrix $A = \begin{bmatrix} 634.156 & 0 & 755912.06 \\ 0 & 1426.8604 & 598151.25\\ 755912.06 & 598151.25 & 1.1517e9\\ \end{bmatrix} $ with eigenvalues $\begin{bmatrix} 1.15179e9\\ 1254.2858\\ -1.52588e-5\\ \end{bmatrix}$. I want to calculate the Cholesky decomposition of the matrix A but A is not positive definite (the last eigenvalue of A is negative). How can I transform A into a nearest positive definite matrix inorder to calculate the Cholesky decomposition?

Answer 1

@ Muhammad Omer , if $A$ is your exact result (and not an approximation), then I think that your work is not serious for the following reasons:

1. The $a_{i,j}$ are known with $8,6$ or $5$ significant digits ; moreover the most important entry (considering the precision) is $a_{3,3}$ that is known with only $5$ digits.

2. The $<0$ eigenvalue of $A$ is $\approx -0.06$.

3. Since you know that the matrix is (in reality) SPD, then the principal question is: after a small modification of $A$, what will be the number of significant digits of the result $C$ s.t. $A=CC^T$ ?

For instance, if we replace $a_{3,3}$ with $1.151795$, then the $<0$ eigenvalue becomes $\approx 2.10^{-4}$. In other words, the digit $7$ is false and must be replaced with $8$ ; therefore, you have only $4$ significant digits. With this modified $A$ ($a_{3,3}=1.151795$), we obtain a matrix $C$ that can be written (with $10$ significant digits) $\begin{pmatrix}25.18245421& 0& 0\\0& 37.77380574& 0\\30017.41029& 15835.08038& 17.58435679\end{pmatrix}$ ; note that $||CC^T-A||\approx 0.37$, that implies that, if we keep $4$ digits for $C$, then the error will be huge (cf. the Steven post: " later computations might be numerically unstable, which may not have the desired effects.").

Conclusion: in my opinion, $C$ is known with $0$ significant digit.

Answer 2

If you have a square matrix there are several ways to find the nearest symmetric positive semidefinite matrix outlined by https://nhigham.com/2021/01/26/what-is-the-nearest-positive-semidefinite-matrix/

Then you can add a small perturbation learned from the eigenvalues to bring the final eigenvalues above 0 or a given machine precision or error level. https://nhigham.com/2021/02/16/diagonally-perturbing-a-symmetric-matrix-to-make-it-positive-definite/

Here is pseudocode:

The simplest for an input symmetric matrix A is:

   set a machine precision or error level above 0:
       eps = 1.e-17 or 1.e-11 etc.

   get an orthogonal matrix from eigen decomposition of A:
       evd = EVD.factorize(A).
       leftEigenVec = evd.left;

   replace the eigenvalues smaller than eps with eps:
       tau = evd.realEigenvalues
       for (i = 0; i < tau.length; ++i)
           if (tau[i] < eps)
               tau[i] = eps;
   then the positive semidefinite matrix is   
       aPSD = leftEV * tau * leftEV^T
   
If you have a square matrix that is not necessarily symmetric:

   calculate B, the the symmetric part of A:
       B = 0.5*(A + A^T)
 
   calculate the polar decomposition of B to get Q and H 
   (where H is formed from V, S, and V^T of SVD(B) 
    within the polar decomposition method):
       [q, h] = polarDecomposition(b);

   the symmetric positive semidefinite matrix is then:
       aPSD = (A + h)/2

Then a small perturbation may need to be applied 
to result in eigenvalues greater than or equal to eps:
  evd = EVD.factorize(aPSD);
  eig = evd.realEigenvalues;
  if any eig < eps:
      sort(eig);
      eigMin = max(-eig[0], eps);
      eigMin = min(eigMin, 0)
      D = MatrixUtil.zeros(A.length, A.length);
      for (i = 0; i < eig.length; ++i)
          D[i][i] = eigMin;
      aPSD = aPSD + D;