Cfrgtkky

If I want to calculate $B^{-1}A$, then instead of using Inverse, I should in theory just be able to use LinearSolve[B,A].

Now if I have $B$ is positive definite and symmetric, then $A^TB^{-1}A$ should also be PD and Symmetric. But it's not...

Here's an example of what I mean:

n = 1000;

m = 10;

a = RandomReal[{-1, 1}, {n, m, m}];

b = RandomVariate[WishartMatrixDistribution[11, IdentityMatrix[10]], 

   n];

result = ParallelTable[

   a[[i]].LinearSolve[b[[i]], Transpose[a[[i]]]], {i, 1, Length[a]}];

And then I get

In[323]:= Total[PositiveDefiniteMatrixQ /@ result]

Total[SymmetricMatrixQ /@ result]



Out[323]= 1000 True



Out[324]= 388 False + 612 True

This is a precision problem. If I use Inverse the problem worsens.

However, this seems like something that could have been corrected... Or is it something very difficult to do? I don't get any warning about badly-conditioned matrices.

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.