# Gaussian sampling in high dimension

I have a covariance function $f(x)$, where $x = (x_1, x_2, x_3)$ is a point in three-dimensional space. I need to generate a Gaussian field with given covariance function on a 3D grid of points, that is very large (the storage of full covariance matrix is impossible, so a Cholesky decomposition is fail). What is the best approach to deal with this kind of problem?

Thanks,

Does the covariance between two points $x$ and $y$ depend only on $x-y$? If so, then a spectral method is one way to approach the problem. This is discussed in many textbooks on geostatistics.
• I'm curious what you mean by "fails." Obviously, to get an adequate approximation to your spectral density you may have to increase $N$ (and then subsample the resulting vectors), but I don't see any fundamental reason why this shouldn't work. May 28, 2014 at 20:07
• You're right. I have to increase $N$ to obtain Gaussian distribution according to the central limit theorem. But the problem is, the coefficient $A_i$ of the series $A_i\cos(\omega_i t+\phi_i)$ tends to zero quite fast. So the fact increasing N is not change anything... May 29, 2014 at 11:57
Note, that there exists a large class of methods which allows large matrix multiplications for the data of specific structure. As you have a 3D grid of points, you have a separate fixed grid for each dimension of sizes $n_i, i = \overline{1, 2, 3}$ . So, if you can transform your covariance function $f(x, y)$, both $x, y \in \mathbb{R}^3$ into product of three covariance functions: $f(x, y) = \prod_{i = 1}^3 f_i(x_i, y_i)$ then there exist fast methods for Gaussian processes regression.