# Long samples from Gaussian Process _prior_

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.

• Perhaps sparse data dimensionality reduction would be of interest. – Carl Nov 13 '18 at 22:35
• You can use the sparse gaussian approximations in GPy. See here. – InfProbSciX Nov 14 '18 at 14:11
• @InfProbSciX: This sounds quite like what I want, but I don't fully follow what has been done in what you have linked to. If I understand correctly, they are using the sparsity to more efficiently generate the posterior distribution. Could you tell me how to do that easily for the prior? I'm quite new to GPy so I apologize for the naivety. – SarthakC Nov 14 '18 at 19:53
• To clarify, I don't have any data points, I just want a long stochastic curve that is a Gaussian process with mean zero and has a squared exponential kernel – SarthakC Nov 14 '18 at 19:55

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.

• This looks extremely helpful! I was also not aware of Kernels with compact support with are positive definite, thanks for that. I think algorithm 2 should work for my purpose. I will have a better look at this when I get to work today and will mark your answer as the accepted one then. – SarthakC Nov 15 '18 at 15:24