# Using LinearSolve instead of Inverse does not give a good enough precision

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],
n];
result = ParallelTable[
a[[i]].LinearSolve[b[[i]], Transpose[a[[i]]]], {i, 1, Length[a]}];


And then I get

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

Out= 1000 True

Out= 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.

• What do you get with Total[SymmetricMatrixQ[#, Tolerance -> 1*^-10] & /@ result]? – Michael E2 Feb 13 '19 at 16:17
• @MichaelE2 I get 1000 True. The problem is that other in-built functions will 'complain'. For example, if I use the matrices above as covariance matrices in MultinormalDistribution, I get an error message. And I guess it's not the only one... – An old man in the sea. Feb 13 '19 at 16:21

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.

• I get Total[SymmetricMatrixQ /@ result2] --> 163 False + 837 True, the reverse of what you're showing, if I start with SeedRandom. – Michael E2 Feb 13 '19 at 16:34
• If I use LinearSolve[b[[i]], Transpose[a[[i]]], Method -> "Cholesky"] in the OP's code, then I get 1000 True... – Michael E2 Feb 13 '19 at 16:36
• Huh, that is indeed weird. Version 11.3 for macOS delivers 1000 True also for SeedRandom... – Henrik Schumacher Feb 13 '19 at 16:38
• @MichaelE2, when I use Method -> Cholesky, I don't get that!! I still get many Falses... I'm using Windows 7, with Mathematica Version 11.0.1 – An old man in the sea. Feb 13 '19 at 16:39
• That weird. I'm using V11.3 (but the early one) on MacOS, too: i.stack.imgur.com/yfSnf.png – Michael E2 Feb 13 '19 at 16:40

The problem is rounding error. The option Tolerance is one way Mathematica lets you deal with it in some functions.

SeedRandom;
n = 1000;
m = 10;
a = RandomReal[{-1, 1}, {n, m, m}];
b = RandomVariate[WishartMatrixDistribution[11, IdentityMatrix], 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.