Using LinearSolve instead of Inverse does not give a good enough precision
The problem is rounding error. The option Tolerance
is one way Mathematica lets you deal with it in some functions.
SeedRandom[0];
n = 1000;
m = 10;
a = RandomReal[{-1, 1}, {n, m, m}];
b = RandomVariate[WishartMatrixDistribution[11, IdentityMatrix[10]], n];
result = Table[a[[i]].LinearSolve[b[[i]], Transpose[a[[i]]]], {i, 1, Length[a]}];
Total[PositiveDefiniteMatrixQ /@ result]
Total[SymmetricMatrixQ[#, Tolerance -> 1*^-10] & /@ result]
(*
1000 True
1000 True
*)
Here's why 1*^-10
is about the right tolerance:
ListPlot[#, PlotLabel -> Max[#]] &[
$MachineEpsilon*Through[(LinearSolve /@ b)["ConditionNumber"]]
]
The condition number estimates the relative loss of precision due to LinearSolve
. The maximum absolute value of the entries of the matrices b
are near 1
, so $MachineEpsilon
is the right number to use to multiply the condition number. There is further rounding error in the matrix multiplication, so the tolerance needs to be a little greater than the estimated maximum error shown in the plot label above.
Using the CholeskyDecomposition
explicitly not only seems to remove the problem, it is also faster: Moreover, this gets rid of one of the matrix-matrix multiplications and, probably more important, it requires onle one triangular matrix solve instead of two.
Notice also that I cast the matrix-matrix products in a way that Transpose[#].# &
is applied last. This ensures that the resulting matrix is numerically symmetric and positive (semi-)definite.
n = 1000;
m = 50;
a = RandomReal[{-1, 1}, {n, m, m}];
b = RandomVariate[WishartMatrixDistribution[m + 1, IdentityMatrix[m]],
n];
result = Table[
a[[i]].LinearSolve[b[[i]], Transpose[a[[i]]]], {i, 1,
Length[a]}]; // AbsoluteTiming // First
result2 = Table[
With[{L = Transpose[CholeskyDecomposition[b[[i]]]]},
Transpose[#].# &[LinearSolve[L, Transpose[a[[i]]]]]
], {i, 1, Length[a]}]; // AbsoluteTiming // First
Total[SymmetricMatrixQ /@ result]
Total[SymmetricMatrixQ /@ result2]
0.437883
0.19313
529 False + 471 True
1000 True
Beware that CholeskyDecomposition
does not apply pivoting which may also cause accuracy problems when small entries appear on the main diagonal.