• API
• FAQ
• Tools
• Archive
SHARE
TWEET # Untitled a guest Jun 25th, 2019 59 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. In:= positivizeMatrix[X_] := Module[{res},
2.    res = 0.5*(X + Transpose[X]);
3.
4.    (*This method only works for Hermitian matrices*)
5.
6.    If[Not[PositiveDefiniteMatrixQ[res]],
7.
8.     auxsystDP = Eigensystem[res];
9.
10.     duDP = DiagonalMatrix[auxsystDP[]];
11.
12.     wDP = ReplacePart[duDP, {i_, i_} /; duDP[[i, i]] <= 0 :> 0];
13.
14.     res = Transpose[auxsystDP[]].wDP.auxsystDP[];
15.     (*Eigensystem dá o eigenvectors em linhas e não colunas*)
16.
17.     lenRes = Length[res];
18.
19.     powerDP = 20;
20.
21.     While[Not[PositiveDefiniteMatrixQ[res]],
22.
23.      res = res + IdentityMatrix[lenRes]*\$MachineEpsilon*2^(-powerDP);
24.      (*Adding the smallest nugget possible.
25.      This may not be a good idea. It still creates machine-
26.      precision problems... See below. *)
27.
28.      powerDP--;
29.      ]; ,
30.     res
31.     ];
32.
33.    (*Theoretically the next block of code shouldn't be needed,
34.    but due to numerical machine accuracy, we sometimes need this*)
35.
36.    If[Not[SymmetricMatrixQ[res]],
37.     res = 0.5*(res + Transpose[res]);
38.     ];
39.
40.    res
41.    ];
42.
43. In:= positivizeMatrix2[X_] := Module[{res},
44.    res = 0.5*(X + Transpose[X]);
45.    (*This method only works for Hermitian matrices*)
46.
47.    If[Not[PositiveDefiniteMatrixQ[res]] || Det[res] <= 0,
48.
49.     auxsystDP = Eigensystem[res];
50.
51.     duDP = DiagonalMatrix[auxsystDP[]];
52.
53.     wDP = ReplacePart[duDP, {i_, i_} /; duDP[[i, i]] <= 0 :> 0];
54.
55.     res = Transpose[auxsystDP[]].wDP.auxsystDP[];
56.     (*Eigensystem dá o eigenvectors em linhas e não colunas*)
57.
58.     lenRes = Length[res];
59.
60.     powerDP = 20;
61.
62.     While[Not[PositiveDefiniteMatrixQ[res]] || Det[res] <= 0,
63.
64.      res = res + IdentityMatrix[lenRes]*\$MachineEpsilon*2^(-powerDP);
65.      (*Adding the smallest nugget possible.
66.      This may not be a good idea. It still creates machine-
67.      precision problems... See below. *)
68.
69.      powerDP--;
70.      ]; ,
71.     res
72.     ];
73.
74.    (*Theoretically the next block of code shouldn't be needed,
75.    but due to numerical machine accuracy, we sometimes need this*)
76.
77.    If[Not[SymmetricMatrixQ[res]],
78.     res = 0.5*(res + Transpose[res]);
79.     ];
80.
81.    res
82.    ];
83.
84. In:= SeedRandom;
85. list = RandomReal[{0, 50}, {10000, 3, 3}];
86.
87. In:= res1 = ParallelMap[positivizeMatrix[#] &, list];
88. res2 = ParallelMap[positivizeMatrix2[#] &, list];
89.
90. In:= ParallelMap[PositiveDefiniteMatrixQ, res1] // Total
91.
92. Out= 10000 True
93.
94. In:= ParallelMap[Det[#] > 0 &, res1] // Total
95.
96. Out= 1749 False + 8251 True
97.
98. In:= ParallelMap[Det[#] > 0 &, res2] // Total
99.
100. Out= 10000 True
101.
102. In:= ParallelMap[Quiet[Inverse[#]] &, res2] // Total;
103.
104. During evaluation of In:= Inverse::sing: Matrix {{39.1928,33.6862,14.4105},{33.6862,<<18>>,<<18>>},{14.4105,32.7791,25.9402}} is singular.
105.
106. During evaluation of In:= Inverse::sing: Matrix {{21.7121,15.7919,28.2016},{15.7919,<<18>>,<<18>>},{28.2016,27.5316,40.2261}} is singular.
107.
108. During evaluation of In:= Inverse::sing: Matrix {{29.7926,31.3398,32.2753},{31.3398,<<18>>,<<18>>},{32.2753,33.9515,34.9649}} is singular.
109.
110. During evaluation of In:= General::stop: Further output of Inverse::sing will be suppressed during this calculation.
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top