Speeding a function that given a matrix, returns a close by Positive Definite Matrix

Question

I'm dealing with several matrices, that due to machine-precision considerations, they are supposed to be symmetric Positive definite (SPD), but aren't.

From matrix $X$, we proceed as:

1. Find $Y=\frac12 (X+X^\top)$, the closest symmetric matrix to $X$.
2. Take an eigendecomposition $Y=QDQ^\intercal$, and form the diagonal matrix $D_+=\max(D,0)$ (elementwise maximum).
3. Find the smallest $\epsilon>0$ such that the computer recognizes the following as a PD matrix:$Z=Q(D_++\epsilon I)Q^\intercal$. Theoretically, any positive value should be valid. However, due to machine precision, there is a lower bound.

Then we just use matrix $Z$ as your closest PD matrix.

For this method, I constructed the following function:

positivizeMatrix[X_] := Block[{res},
  
  res = X;
  
  res = 0.5*(res + Transpose[res]);
  (*This method only works for Hermitian matrices*)
  
  If[Not[PositiveDefiniteMatrixQ[res]] || 
    MatrixRank[res] != Length[res],
   
   auxsystDP = Eigensystem[res];
   
   duDP = DiagonalMatrix[auxsystDP[[1]]];
   
   wDP = ReplacePart[duDP, {i_, i_} /; duDP[[i, i]] <= 0 :> 0];
   
   lenRes = Length[res];
   
   powerDP = 0;
   
   While[Not[PositiveDefiniteMatrixQ[res]] || 
     MatrixRank[res] != Length[res],
    
    res = 
     Transpose[
       auxsystDP[[2]]].(wDP + 
        IdentityMatrix[lenRes]*$MinMachineNumber*10^
          powerDP).auxsystDP[[2]];
    
    powerDP++; 
    ]; 
   ];
  res
  ]

Since I'm working with matrices supposed to be SPD, and because Mathematica's SymmetricMatrixQ and PositiveDefiniteQ often fail with 'my' matrices (see for example this old question of mine), the condition Not[PositiveDefiniteMatrixQ[res]] || MatrixRank[res] != Length[res] was the only one which I managed to find so that in posterior computations, when I'm drawing from an inverse wishart distribution (or using the Cholesky Method in LinearSolve), no error message pops up.

I think the main bottleneck is that, even though the matrix is close to being SPD, we need to iterate through the while cycle several times, each iteration being very time consuming.

The objective is to speed up this function.

I may consider answer that diverges from the numerical algorithm that I explained above, as long as the algorithm is as general as this one.

P.S.: I'll include a bonus for this question (at the end) of at least 100 points. Depending on the quality of the answer, it may be more.

Edit: here's simple criteria that your programme must pass...

mat = {{313.95758323354806, 443.4761851932803, 
    274.8711923640176}, {443.4761851932803, 671.32551818649, 
    417.77978277045105}, {274.8711923640176, 417.7797827704511, 
    260.0520353920102}};
mat4 = positivizeMatrix[mat]
MatrixRank[mat4]

The output must be 3, otherwise, it most likely won't work on my big programme where this small function will be called hundreds of thousands of times.

Answer 1

Not sure why you use the While loop. A faster strategy might be to do the change proportional to the average of the positive eigenvalues. $MinMachineNumber is a very small number; maybe it is better to use MachineEpsilon instead. That is the smallest number will have an actual effect when added to $1$. I am also not sure why you require calling PositiveDefiniteMatrixQ. Shouldn't the numerically obtained eigenvalues all being greater than $0$ be enough of a certificate for you?

f[X_] := Module[{\[Lambda], U},
  {\[Lambda], U} = Eigensystem[0.5*(X + Transpose[X])];
  Transpose[U].Times[Clip[\[Lambda], {Mean[Ramp[\[Lambda]]] $MachineEpsilon, \[Infinity]}], U]
  ]

And if you really insist on doing the While loop for finding the smallest regularization, then do it with a binary search like this:

g[X_] := Module[{\[Lambda], U, a, b, c, iter},
  {\[Lambda], U} = Eigensystem[0.5*(X + Transpose[X])];
  \[Lambda] = Ramp[\[Lambda]];
  a = $MinMachineNumber;
  b = Max[\[Lambda]];
  c = Sqrt[a] Sqrt[b];
  iter = 0;
  While[a < 0.5 b,
   iter++;
   If[PositiveDefiniteMatrixQ[Transpose[U].((\[Lambda] + c) U)],
    b = c,
    a = c
    ];
   c = Sqrt[a] Sqrt[b];
   ];
  Print[iter, " iterations required, \[Epsilon] = ", b];
  Transpose[U].((\[Lambda] + b) U)
  ]