# 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 =
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. • Why not just use the lower triangular portion of the matrix? – Asim May 19, 2021 at 17:23 • @Asim Maybe I could... I'm not sure of the theoretical implications of just using the lower triangular portion of the matrix though... Do you have some bibliography? thanks ;) May 19, 2021 at 18:46 • It depends upon what you want to do with the PD matrix. If it is to draw from a Inverse Wishart distribution, then just taking the lower triangular portion and then replicating this on the upper triangle should work. If the matrices are small, you could just write your own code for the Inverse Wishart rather than using the Mathematica function. This is likely to be faster and will give less headache. – Asim May 19, 2021 at 19:28 • @Asim I'm not following you... Do you mean that I should take the lower triangular part of the almost SPD matrix, and then flip it and add it with itself (L+U) to create a SPD matrix? May 20, 2021 at 9:47 • What is your distance measure when you mention "close by" PDM? – A.G. May 20, 2021 at 11:53 ## 1 Answer 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)
]

• Henrik, thanks for your answer. I've tried both of your programmes, and they both fail when incorporate them on in a bigger programme that I have... I've edited the original question to include a sufficient criteria that it should meet. When I try with your functions, the rank of the returned matrix is still 2, and this will cause problems later on. May 20, 2021 at 14:30