Failure in MultinormalDistribution

My RandomVariate fails to evaluate when sampling from a multinormal distribution, any idea what's happening here?

This is Mathematica 12.1 on Mac:

d = 6;
sigma = {{0.2456589529728737, 0.06467112508291382,
0.07381608109139866, 0.01761648518373636, 0.05925859721052876,
0.045898960400932734}, {0.06467112508291387,
0.35501955566599797,
0.05242856376555134, -0.009970312195220573,
0.015359942596352098,
0.02448314785173377}, {0.07381608109139871,
0.05242856376555134, 0.3026715722222353, -0.0819524814567561,
0.07919970253622352, -0.033627335100296744}, \
{0.01761648518373638, -0.009970312195220557, -0.08195248145675606,
0.5519153257149109, -0.1068423071796184,
0.35028088018797104}, {0.05925859721052879,
0.015359942596352072,
0.07919970253622348, -0.10684230717961839, 0.29853333010040306,
0.002543954098042138}, {0.0458989604009327,
0.02448314785173373, -0.033627335100296814,
0.35028088018797104, 0.002543954098042013, 0.6962012633235793}};
ii = IdentityMatrix[d];
RandomVariate@MultinormalDistribution[sigma] (* fails *)
RandomVariate@MultinormalDistribution[sigma + ii] (* works *)


The example above was generated by code below, it's a random covariance matrix with geometrically decaying eigenvalues, so not even that close to singular --

randomRotation[n_] := Module[{M, z, q, r, d, ph},
z = RandomVariate[NormalDistribution[0, 1], {n, n}];
{q, r} = QRDecomposition[z];
d = Diagonal[r];
ph = d/Abs[d];
M = q*ph; (* Note, determinant may be -1 giving "pseudo-rotation",
switch 2 rows of the matrix to guarantee true rotation *)
indices = If[Det[M] > 0, Range[n], {2, 1}~Join~Range[3, n]];
M[[indices]]
];

d = 20;
rot = randomRotation[d];
ii = IdentityMatrix[d];

sigma = rot.DiagonalMatrix[Table[1/i, {i, 1, d}]].Inverse[rot];
RandomVariate@MultinormalDistribution[sigma] (* fails *)
RandomVariate@MultinormalDistribution[sigma+ii] (* works *)
$$$$

• "...symmetric positive definite..." for the covariance. Sigma does not meet the condition. That said, an error should be for that, instead it seems to fail in the definition of the distribution itself (that is, MultinormalDistribution[sigma] should note the condition failure). You could try something like RandomVariate@MultinormalDistribution[N[Rationalize[sigma, 0], 8]] – ciao Oct 27 '20 at 22:43
• ah, thanks...it looks like symmetry failure from numerical noise. On Mac the error was not very informative, it just said MultinormalDistribution called with 1 argument; 2 arguments expected` – Yaroslav Bulatov Oct 27 '20 at 22:58
• For whatever it's worth on Windows one gets the same 1 vs 2 arguments error. Also, see the answer to mathematica.stackexchange.com/questions/152987/… which gives a function to force the approximately symmetric matrix to be symmetric. – JimB Oct 28 '20 at 0:27