# Contradictions in PositiveDefiniteQ, Det and Inverse in Mathematica

A few words before the working example. Below we find 2 functions which should theoretically transform any matrix to its closest Symmetric Positive Definite (SPD) matrix: positivizeMatrix, and positivizeMatrix2. The difference between both is that the second one has the extra condition Det[res]<=0. Any SPD is also invertible.

   positivizeMatrix[X_] := Module[{res},
res = 0.5*(X + Transpose[X]);

(*This method only works for Hermitian matrices*)

If[Not[PositiveDefiniteMatrixQ[res]],

auxsystDP = Eigensystem[res];

duDP = DiagonalMatrix[auxsystDP[[1]]];

wDP = ReplacePart[duDP, {i_, i_} /; duDP[[i, i]] <= 0 :> 0];

res = Transpose[auxsystDP[[2]]].wDP.auxsystDP[[2]];
(*Eigensystem dá o eigenvectors em linhas e não colunas*)

lenRes = Length[res];

powerDP = 20;

While[Not[PositiveDefiniteMatrixQ[res]],

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

res = 0.5*(res + Transpose[res]);
powerDP--;
]; ,
res
];

res
];


The examples use 3 by 3 matrices where each component is between 0 and 50. I assumed it wouldn't be at the edge of machine precision. Now, for the results:

In[48]:= SeedRandom[1234];
list = RandomReal[{0, 50}, {10000, 3, 3}];

In[53]:= res1 = ParallelMap[positivizeMatrix[#] &, list];
res2 = ParallelMap[positivizeMatrix2[#] &, list];

In[55]:= ParallelMap[PositiveDefiniteMatrixQ, res1] // Total

Out[55]= 10000 True

In[58]:= ParallelMap[Det[#] > 0 &, res1] // Total

Out[58]= 1749 False + 8251 True

In[59]:= ParallelMap[Det[#] > 0 &, res2] // Total

Out[59]= 10000 True


This means I need to add the extra condition of Det[res]<=0.

Now the problem with Inverse, with matrices which have positive Det:

In[63]:= ParallelMap[Quiet[Inverse[#]] &, res2] // Total;

During evaluation of In[63]:= Inverse::sing: Matrix {{39.1928,33.6862,14.4105},{33.6862,<<18>>,<<18>>},{14.4105,32.7791,25.9402}} is singular.

During evaluation of In[63]:= Inverse::sing: Matrix {{21.7121,15.7919,28.2016},{15.7919,<<18>>,<<18>>},{28.2016,27.5316,40.2261}} is singular.

During evaluation of In[63]:= Inverse::sing: Matrix {{29.7926,31.3398,32.2753},{31.3398,<<18>>,<<18>>},{32.2753,33.9515,34.9649}} is singular.

During evaluation of In[63]:= General::stop: Further output of Inverse::sing will be suppressed during this calculation.


My point is: Which criteria can I use in my function positivizeMatrix such that I can be sure, whatever matrix I obtain, it can be used for Inverse, and any other in-built Mathematica function which theoretically should work with SPD matrices?

P.S: Even if the components remain between 0 and 10, the problem still subsists. I've just tried it. Since the results are similar, I won't change the text above. However,it makes me wonder if it really is a machine-precision problem...