0
$\begingroup$

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...

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.