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

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 = 
       auxsystDP((2))).(wDP + 

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.