0
$\begingroup$

I use the following functions to receive a matrix, and make it symmetric positive definite:

isPD[mat_] := Module[{norm},
   norm = Norm[mat];
   Apply[And, Map[(# > 0) &, Eigenvalues[mat/norm]*norm]]];

positivizeMatrix[X_] := Module[{res},
   res = X;
   res = 0.5*(res + Transpose[res]);
   norm = Norm[res];
   
   If[Not[isPD[res]],
    auxsystDP = Eigensystem[res/norm];
    duDP = DiagonalMatrix[auxsystDP[[1]]]*norm;
    wDP = ReplacePart[duDP, {i_, i_} /; duDP[[i, i]] < 0 :> 0];
    lenRes = Length[res];
    powerDP = 30;
    While[Not[isPD[res]],
     res = 
      Transpose[
        auxsystDP[[2]]].(wDP + 
         IdentityMatrix[
           lenRes]*$MinMachineNumber*2^(-powerDP)).auxsystDP[[2]];
     res = 0.5*(res + Transpose[res]);
     powerDP--;
     ]; 
    ];
   res
   ];

Next, I have a function that simply does v.LinearSolve[mat, v, Method -> "Cholesky"];, where v is a constant vector, and mat is 'SPDized' matrix using the above positivizeMatrix.

Here an example where they don't work together...

LinearSolve[{{0.2767887221630757`, -0.20435493943670538`, 
   0.17471789061898896`}, {-0.20435493943670538`, 
   0.25577623175273967`, -0.243257996494102`}, {0.17471789061898896`, \
-0.243257996494102`, 
   0.23474884689385866`}}, {-0.7008470974262857`, \
-0.9457472725629061`, -1.654496598836136`}, Method -> "Cholesky"]

returns LinearSolve::npdef: The matrix {{0.276789,-0.204355,0.174718},{<<1>>},{0.174718,-0.243258,0.234749}} is not positive definite, whereas isPD returns True.

How can I create a positivize function that will work for any matrix dim, and be consistently with the linear solve, using the Cholesky method?

Improve this question
$\endgroup$
2
  • $\begingroup$ You use a modified (regularized) definition of positive definiteness. I can imagine a case where Norm is almost zero. The matrix would not be positive definite (eigenvalues essentially zero to machine precision). But if multiplied by a big number (1/Norm) the eigenvalues would become positive. $\endgroup$
    – yarchik
    Dec 1, 2020 at 21:05
  • $\begingroup$ @yarchik even if I do norm=Max[Norm[mat],1], the problem still persists... $\endgroup$
    – An old man in the sea.
    Dec 1, 2020 at 21:34

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.