Long samples from Gaussian Process _prior_

Question

I'm interested in being able to sample a long (N~10^5) sample from a Gaussian process. For a small sample I understand I can quite easily construct an NxN covariance matrix and then choose a random sample from a multivariate gaussian distribution with that covariance. However, when N is around 10^5, the covariance matrix needs to store 10^10 values, which is causing memory issues. I understand that my covariance matrix is quite sparse, in that I only really have values around some small band (I'm using a squared exponential kernel) around the diagonal, but I don't know how to take advantage of that.

Any ideas? I'm coding in Python, so any means to do this in python would be preferable. Maybe there's something from the GPy package?

(This question is essentially a cross post of this)

Thanks.

Answer 1

One way of tackling this problem, is to generate the sample sequentially. This can be done using the formula for the conditioning of Gaussians (See Rasmussen, Williams Formula A.6 in Appendix A.2)

But the choice of squared exponential covariance does not seem ideal to me. It will require some approximations or even fudging. Note that simply thresholding the covariance matrix is probably not what you want to do to make the matrix sparse, as the result might not be positive definite anymore. Preferable were a kernel with compact support. (such as the Wendland Kernels see Rasmussen, Williams Formula 4.21) Only with those Kernels sparsity is ensured.

Say you want to create a sample of the $N$-dimensional multivariate normal random vector $(X_1,\ldots,X_N).$ You can proceed as follows

Algo 1

1. Determine the variance $\sigma^2$ of $X_1$ and create a sample $x_1$ of $X_1$.
2. Determine the covariance and mean vector of $(X_1,X_2)$ and create a sample $x_2$ of $X_2$ conditional on $x_1$.
3. In the general case you have generated $x_1,\ldots,x_k$ and create a sample for $X_{k+1}$ conditional on the values $x_1,\ldots,x_k$.

This Algorithm is exact, i.e. it will produce a multivariate normal distribution with mean and covariance as required. Sadly, it does not solve your problem, since sooner or later the size of your set of conditioning variables $x_1,\ldots,x_k$ will grow beyond a manageable size.

Not always exact is the following approach. Let $n$ be the maximum size of samples you can handle safely in one go.

Algo2

1. Create a sample for $X_1,\ldots,X_{n}$ either by the method above or in one go.
2. Determine the correlations of $X_{n+1}$ with $X_1,\ldots,X_n$. Drop all $X_j$ which have zero correlation. If none have, drop the variable with correlation closest to zero.
3. Determine a sample for $X_{n+1}$ as in Step 3 of Algo 1 but use only the variables not dropped in Step 2 for conditioning.
4. To find $X_k$ with $k>n$, drop all zero correlated variables from your conditioning set. If this is not enough to bring down the size, drop the ones with correlation closest to zero until your conditioning set has size $n$.

This Algorithm is exact as long as you never have more than $n$ variables with non-zero correlation in each step. Once you start dropping variables with low correlations you are only approximating. The result will not be an exact sample from a multivariate normal distribution. How close or far it is is hard to say and will depend on the specific covariance matrix.