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[41]:= 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[[1]]];
  11.  
  12.     wDP = ReplacePart[duDP, {i_, i_} /; duDP[[i, i]] <= 0 :> 0];
  13.  
  14.     res = Transpose[auxsystDP[[2]]].wDP.auxsystDP[[2]];
  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[38]:= 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[[1]]];
  52.  
  53.     wDP = ReplacePart[duDP, {i_, i_} /; duDP[[i, i]] <= 0 :> 0];
  54.  
  55.     res = Transpose[auxsystDP[[2]]].wDP.auxsystDP[[2]];
  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[48]:= SeedRandom[1234];
  85. list = RandomReal[{0, 50}, {10000, 3, 3}];
  86.  
  87. In[53]:= res1 = ParallelMap[positivizeMatrix[#] &, list];
  88. res2 = ParallelMap[positivizeMatrix2[#] &, list];
  89.  
  90. In[55]:= ParallelMap[PositiveDefiniteMatrixQ, res1] // Total
  91.  
  92. Out[55]= 10000 True
  93.  
  94. In[58]:= ParallelMap[Det[#] > 0 &, res1] // Total
  95.  
  96. Out[58]= 1749 False + 8251 True
  97.  
  98. In[59]:= ParallelMap[Det[#] > 0 &, res2] // Total
  99.  
  100. Out[59]= 10000 True
  101.      
  102. In[63]:= ParallelMap[Quiet[Inverse[#]] &, res2] // Total;
  103.  
  104. 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.
  105.  
  106. 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.
  107.  
  108. 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.
  109.  
  110. During evaluation of In[63]:= 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. OK, I Understand
 
Top