# Is this a Symmetric Matrix or not?

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]]];
];
];
mat
];

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.

• Something to note, SymmetricMatrixQ[cov, Tolerance -> 10^-12] returns False, while SymmetricMatrixQ[cov, Tolerance -> 10^-11] returns True. – user6014 Aug 3 '17 at 13:20
• 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. – user6014 Aug 3 '17 at 13:23
• @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}]? – An old man in the sea. Aug 3 '17 at 14:02
• 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. – user6014 Aug 3 '17 at 14:09
• ...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]]; – 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


or

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.