# How to sample from a MultinormalDistribution with a positive semidefinite covariance matrix?

I'm trying to sample from a Gaussian process (GP) following these tutorials:

I translated Python code from (1) into Wolfram/Mathematica:

(* Exponentiated quadratic kernel *)
Outer[{a, b} |-> Exp[-0.5 (a - b)^2], xa, xb]

(* Want to evaluate GP at these points *)
x = With[{nsamples = 41},
Range[-4, 4, Abs[-4 - 4]/nsamples // N]]

(* Covariance matrix *)

(* Sample from multivariate normal *)
MultinormalDistribution[0*x, \[CapitalSigma]] // RandomVariate
(* This is equivalent to sampling from the GP *)


The last line produces this error:

MultinormalDistribution::posdefprm: The value {{1.,0.981144,0.926682,0.842546,0.737433,0.621321,0.503936,0.393459,0.295726,0.213966,<<32>>},<<9>>,<<32>>} at position 2 in MultinormalDistribution[{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,<<32>>},{{1.,0.981144,0.926682,0.842546,0.737433,0.621321,0.503936,0.393459,0.295726,0.213966,<<32>>},<<9>>,<<32>>}] is expected to be a symmetric positive definite matrix.


Basically, it says that \[CapitalSigma] is not a symmetric positive definite matrix.

Indeed, PositiveDefiniteMatrixQ[\[CapitalSigma]] == False, so it's not positive definite.

However, NumPy has no problem sampling from a multivariate normal with this covariance matrix, as shown in (1) and (2). I can reproduce this too:

import numpy as np

xa, xb = xa[:, None], xb[:, None]
return np.exp(-0.5 * (xa - xb.T)**2)

nsamples = 41
x = np.linspace(-4, 4, nsamples)

# No errors, samples fine
ys = np.random.RandomState(5).multivariate_normal(mean=np.zeros_like(x), cov=Sigma)


Wikipedia says that the covariance matrix can be positive semi-definite (PSD). Indeed, \[CapitalSigma] is positive semidefinite:

In[296]:= \[CapitalSigma] = ExponentiatedQuadratic[x, x];

In[297]:= PositiveSemidefiniteMatrixQ[\[CapitalSigma]]

Out[297]= True


Julia's Distributions.jl can sample from this multivariate normal too, although I have to jump through some hoops to achieve this:

using Distributions
import PDMatsExtras

[exp(-0.5 * (a-b)^2) for a in xa, b in xb]

x = range(-4, 4, step=abs(-4-4) / 41) |> collect

# Tell MultivariateNormal that Sigma is PSD
MultivariateNormal(0 * x, PDMatsExtras.PSDMat(Σ)) |> rand
# 42-element Vector{Float64}:
#   0.8988201432760597
#  -0.4700210875624489
#  -0.6040028058243867
#  -0.5706889628081354
#  -0.9903264880614615
#  ...


How do I sample from MultinormalDistribution with a positive semidefinite covariance matrix in Mathematica? The end goal here is to sample from a Gaussian process that evolves really slowly since it has a high lengthscale (it's equal to 1 in this question).

From the help: "MultinormalDistribution allows [Mu] to be any vector of real numbers, and [CapitalSigma] any symmetric positive definite p*p matrix of real numbers with p=Length[[Mu]]"

A simple hack consist in adding to the covariance matrix a small multiple of the identity matrix. This does not change much, but makes the original matrix positive definite.

\[CapitalSigma] =  ExponentiatedQuadratic[x, x] + 10^-15  IdentityMatrix[42];


Now it is positive definite:

PositiveDefiniteMatrixQ[\[CapitalSigma]]
(* True *)


And we can now get samples:

MultinormalDistribution[0*x, \[CapitalSigma]] // RandomVariate

• Wow, I didn't know the "add a positive number to the diagonal to nudge a matrix closer to being positive definite" trick worked with such small numbers too. Everything works fine now, thanks! Dec 15, 2022 at 16:12

The issue is that with nsamples = 41 machine precision isn't enough to guarantee that RandomVariate sees a positive definite covariance matrix. I've changed the 0.5 in the ExponentiatedQuadratic definition to /2 and used the Cholesky decomposition.

(* Generate covariance matrix with more than machine precision numbers *)
ExponentiatedQuadratic[xa_, xb_] := Outer[{a, b} |-> Exp[-(a - b)^2]/2, xa, xb]
n = 41;
k = 4;
x = With[{nsamples = n}, Range[-k, k, 2 k/nsamples]];

(* Determine Cholesky decomposition *)
L = CholeskyDecomposition[Σ];

(* Generate nSim sets of n+1 independent standard normals *)
nSim = 10000;
z = RandomVariate[NormalDistribution[0, 1], {nSim, n + 1}];

(* Construct random samples from the desired multivariate normal *)
y = Table[x + (# . z[[i]] & /@ Transpose[L]), {i, 1, nSim}];

(* Check on means and variances (which should be x and 0.5, respectively *)
m = Mean[y[[All, #]]] & /@ Range[n + 1]
(* {-4.00881, -3.81503, -3.62098, -3.42581, -3.22903, -3.03097, -2.83265,
-2.63519, -2.43913, -2.24421, -2.04963, -1.85462, -1.65893, -1.46296,
-1.26743, -1.07284, -0.879112, -0.685695, -0.491939, -0.297479,  -0.102286,
0.0935449, 0.290046, 0.487398, 0.685661, 0.884427, 1.08277, 1.27965,
1.47447, 1.66738, 1.85909, 2.05059, 2.24285, 2.43666, 2.6324, 2.82984,
3.02813, 3.22611, 3.42283, 3.61788, 3.81144, 4.00405} *)

Variance[y[[All, #]]] & /@ Range[n + 1]
(* {0.502011, 0.499354, 0.495732, 0.492683, 0.49092, 0.490278, 0.490374, 0.49104,
0.492478, 0.494772, 0.497239, 0.498768, 0.498626, 0.496772, 0.494011,
0.491736, 0.491059, 0.492416, 0.49564, 0.499567, 0.50207, 0.501536,
0.498285, 0.494063, 0.490457, 0.48828, 0.488018, 0.489974, 0.493747,
0.497996, 0.501095, 0.502165, 0.501512, 0.50013, 0.499214, 0.500204,
0.504047, 0.509509, 0.513275, 0.512931, 0.509431, 0.505923} *)


Note that 10 might not be large enough in

Σ = N[ExponentiatedQuadratic[x, x], 10];


if n is much larger than 41.