4
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ What do you get with Total[SymmetricMatrixQ[#, Tolerance -> 1*^-10] & /@ result]? $\endgroup$ – Michael E2 Feb 13 at 16:17
  • 2
    $\begingroup$ @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... $\endgroup$ – An old man in the sea. Feb 13 at 16:21
5
$\begingroup$

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"]]
 ]

enter image description here

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.

$\endgroup$
5
$\begingroup$

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.