Mathematica Asked on December 20, 2020
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?
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP