1
$\begingroup$

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

  1. https://peterroelants.github.io/posts/gaussian-process-tutorial/#Sampling-from-prior
  2. https://stephens999.github.io/stat34800/GP_example.html

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

(* Exponentiated quadratic kernel *)
ExponentiatedQuadratic[xa_, xb_] := 
 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 *)
\[CapitalSigma] = ExponentiatedQuadratic[x, x];

(* 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

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

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

# 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

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

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

# 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).

$\endgroup$

2 Answers 2

2
$\begingroup$

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
$\endgroup$
1
  • $\begingroup$ 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! $\endgroup$
    – ForceBru
    Dec 15, 2022 at 16:12
0
$\begingroup$

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]];
Σ = N[ExponentiatedQuadratic[x, x], 10];

(* 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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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