I've generated the cov matrix in the following way:

kernel[x1_, x2_] := Exp[-1/2*Norm[x1 - x2]^2];
Xtest = Range[-5, 5, 2];
n = Length[Xtest];
Xtrain = RandomReal[{-5, 5}, 5];
kmat[x1_, x2_] := Module[{mat},
   n = Length[x1];
   n2 = Length[x2];
   mat = ConstantArray[0, {n, n2}];
   For[i = 1, i <= n, i++,
    For[j = 1, j <= n2, j++,
      mat[[i, j]] = kernel[x1[[i]], x2[[j]]];

m = kmat[Xtest, Xtrain] // N[#, {Infinity, 1000}] &;

mean = (m.Inverse[kmat[Xtrain, Xtrain]].fobs) // 
   N[#, {Infinity, 1000}] &;
cov = (kmat[Xtest, Xtest] - 
     m.Inverse[kmat[Xtrain, Xtrain]].Transpose[m]) // 
   N[#, {Infinity, 1000}] &;

Theoretically, cov should be symmetric. However, when I do SymmetricMatrixQ[cov], it returns False. It's the m.Inverse[kmat[Xtrain, Xtrain]].Transpose[m] which returns a non-symmetric matrix when it should not. When I do SymmetricMatrixQ[Inverse[kmat[Xtrain, Xtrain]]] I get True.

My objective is to be able to run RandomVariate[MultinormalDistribution[mean, cov], 4]; which I can't, since Mathematica thinks it's not symmetric or PD...

Any help would be appreciated.

  • 3
    $\begingroup$ Something to note, SymmetricMatrixQ[cov, Tolerance -> 10^-12] returns False, while SymmetricMatrixQ[cov, Tolerance -> 10^-11] returns True. $\endgroup$ – user6014 Aug 3 '17 at 13:20
  • $\begingroup$ Evaluate this to see which values trigger the False: Table[{i, j, cov[[i, j]] == Transpose[cov][[i, j]]}, {i, 1, 6}, {j, 1, 6}]. I'll leave it up to you to sort out how that happened. In playing with this, a handful of times I have gotten a cov that SymmetricMatrixQ returned True for. $\endgroup$ – user6014 Aug 3 '17 at 13:23
  • $\begingroup$ @user6014 Thanks for the comments. I've checked and the values indeed differ, when they should not. Why isn't enough to just use N[,{infinity, big number}]? $\endgroup$ – An old man in the sea. Aug 3 '17 at 14:02
  • $\begingroup$ I don't think this tells the whole story, but applying your N statement retroactively doesn't fix precision discrepancies/losses that happened in the earlier (kmat[Xtest, Xtest] - m.Inverse[kmat[Xtrain, Xtrain]].Transpose[m]) computation. I'd have to really sit down and look at it to see if that's actually what's going on here, though. $\endgroup$ – user6014 Aug 3 '17 at 14:09
  • 3
    $\begingroup$ ...you noticed that Xtrain = RandomReal[{-5, 5}, 5] generates its results in machine precision, didn't you? (Look up the WorkingPrecision option.) Also, consider using DistanceMatrix[] and LinearSolve[]: m = Exp[-DistanceMatrix[Xtest, Xtrain]^2/2]; lsf = LinearSolve[Exp[-DistanceMatrix[Xtrain]^2/2]]; mean = m.lsf[fobs]; cov = Exp[-DistanceMatrix[Xtest]^2/2] - m.lsf[Transpose[m]]; $\endgroup$ – J. M. will be back soon Aug 4 '17 at 0:41

Here's what I do in that situmation (which comes up quite often):

cov = .5 * (cov + Transpose[cov]);

If you want to see if a matrix a is symmetric, subtract it from its transpose, and see if it's zero:

a - Transpose[a] // Norm


a - Transpose[a] // Max

depending on how you want to measure the "distance away from symmetry".

So for a = {{1, 2}, {2, 3}} you will get zero, for a = {{1, 2}, {2.0001, 3}}, you will get a small number. This is what is probably happening in your construction, some element(s) i,j are just slightly different from the corresponding j,i element.


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.