Efficient Simulation of Stationary Multivariate Gaussian Random Fields with Given Cross-Covariance

The present paper introduces a new approach to simulate any stationary multivariate Gaussian random field whose cross-covariances are predefined continuous and integrable functions. Such a field is given by convolution of a vector of univariate random fields and a functional matrix which is derived by Cholesky decomposition of the Fourier transform of the predefined cross-covariance matrix. In contrast to common methods, no restrictive model for the cross-covariance is needed. It is stationary and can also be reduced to the isotropic case. The computational effort is very low since fast Fourier transform can be used for simulation. As will be shown the algorithm is computationally faster than a recently published spectral turning bands model. The applicability is demonstrated using a common numerical example with varied spatial correlation structure. The model was developed to support simulation algorithms for mineral microstructures in geoscience.


Introduction
The theory of uni-and multivariate random fields has been used extensively in a broad range of science and engineering disciplines such as meteorology, astrophysics and geosciences ( [1] [2]) over the last decades.The multivariate case is still topic of modern spatial statistics research ([3] [4] [5]).
Of particular interest is the general case when the components of the multivariate vector are not independent.Numerous models of the resulting so-called cross-covariance and methods for simulation of such fields have been given in the literature (see [6]).
Another distinction is whether the approach is restricted to the simulation of the fields at regular grid locations (e.g. the ones mentioned for exactly Gaussian) or not (e.g. the ones for approx.Gaussian), the Cholesky decomposition approach ( [15] [16]) and sequential algorithms ( [17]).Unfortunately the later became computationally too costly when the number of locations exceeds a few thousands.
Many models are also limited to certain parametric models for the cross-covariances.
Such as the methods based on convolution, e.g., the so-called kernel convolution and covariance convolution approach ( [6], p. 739).Both are suitable for the isotropic case but the number of parameters is limited to N (or Recently, a very interesting ansatz without this limitation was given in [3].Similar to a kind of spectral turning-bands algorithm it is available for non-grid positions but it creates only approximately Gaussian fields. The spectrum convolution approach proposed in this paper creates exact stationary Gaussian random fields.Such a multivariate field with values in N  is given by con- volution of a N dimensional vector of univariate Gaussian random fields and a functional N N × matrix.Note that this model is completely different from the Cholesky decomposition approach ([15] [16]) and requires only decomposition of N N × matrices.The functional matrix is derived by Cholesky decomposition of the Fourier transform of the predefined cross-covariance matrix (with some normalization).
It will turn out that this approach allows fast generation of large samples since only Fourier transform and matrix multiplication of the initial vector of univariate Gaussian fields is required.Unfortunately, in the present form, the method seems only appropriate for grid locations.Relevant for many research fields, the new model will be used in our research group in a geostatistical framework.It serves as a basis for models of mineral microstructure.For that a very fast algorithm is required since the grids contain several million locations.
In a certain sense it is similar to coregionalization ( [18]), which may also be defined for a functional transformation matrix but unfortunately this approach cannot be made stationary (besides trivial cases).
The outline of the present paper is as follows.Section two introduces notation and provides the theoretical frameworks and results for multivariate random fields as well as a closer look at two other approaches.This is followed by a description of the new model in section three, which contains the proof of validity.Example simulation and a comparison of the computation effort to another model is done in the second last section.The paper is topped off with a conclusion section.

Theory
This section is intended to make the reader familiar with notation, the mathematical objects used in this paper and some theoretical results.

Multivariate Random Fields
A N-variate second-order random field with existing second moments.
Common first and second-order characteristics are the field trend of Z and the matrix ( ) which is called cross-covariance function of Z.This function satisfies positive semidefiniteness in the sense that , 0 If the cross-covariance function depends on x and y only through the distance vector h x y = − , then Z is called stationary.
Without loss of generality, throughout the rest of the paper, the trend m of the random fields is assumed to be constant equal to ( ) 0, , 0 o =  .In contrast to the univariate case, isotropy has different interpretations in the multivariate case (see [4]).In this paper Z is called isotropic if ( ) A multivariate Gaussian field is a random field where all finite dimensional distributions are normal distributions.Similar to the univariate case it is fully described by the aforementioned characteristics m and ρ .
The main theorem of this paper makes use of the famous Kolmogorov-Chentsov theorem.See [19] (thm 2.2.3) or [20] (thm 1) for the proof.It provides a useful cri- terion for establishing the existence of versions of stochastic processes or fields with continuous sample paths.
Theorem 1 (Kolmogorov continuity) Suppose that the random field ( ) on an open domain d T ⊆  satisfies the following condition: There exist positive constants , ,D α β such that E , , .
Then there exists a version of X on the closure of T with a.s.continuous paths.
Immediately one can deduce the following result.
Corollary 1 A stationary random field ( ) T ⊆  with continuous and integrable covariance function ρ has a version with a.s.continuous paths on T.
Proof.On the one hand, by continuity, ρ is bounded on every bounded subset of T.
On the other hand, it holds ( ) From Theorem 1 follows the statement.
The following theorem from [21] is required for the proof of Lemma 2.
Lemma 1 Suppose ( ) X ∈ is a sequence of normal distributed vectors and : lim E and : lim cov exist then X is normal with mean vector b and covariance matrix C.
Also important for our main result is the next lemma.

Previous Simulation Approaches
As introduced in the beginning there are several procedures to construct and sample random fields with predefined cross-covariance ( [5] and [18] section 28.9).Two promising ways are given in detail below.
a) A common method to correlate several random fields is called coregionalization (see [4] section four).For this, let  be a vector of independent stationary and isotropic Gaussian random fields.For a matrix M define the vector of random fields Z by The cross-correlation of Z is now given by the positive semi-definite matrix ( [4]) b) Recently, an improved spectral turning-bands algorithm for simulating stationary multivariate Gaussian random fields was presented in [3].Using this one can simulate any multivariate Gaussian field whose cross-covariance function is continuous and absolutely integrable for each entry.
The approach reformulated and reduced in what follows for the stationary isotropic case in 2   which makes the cross-covariance symmetric such that the imaginary part in the model vanishes.In general it is also available for the anisotropic case.
Extending the very famous result of Bochner for the multivariate case the isotropic cross-covariance at lag : r h = can be written as (see [22]) is called the matrix of spectral densities ( H + denotes the set of symmetric positive semi-definite matrices).It is given by which can be rewritten using Hankel transform (see [23]) as where 0 J is the 0-th order Bessel function of the first kind.
As shown in [3] the random field Z is approximately multivariate Gaussian (for large L) with trend o and cross-covariance ρ .

Spectrum Convolution Approach
The main part of the paper is given in what follows.The following model can be seen as a spectral variant of coregionalization (Section 2.2) but also bears analogy to the spectral turning bands method.Let ( ) For the spectrum it holds ( ) ( ) ( ) ( ) • Let d  be the d-dimensional hypertorus, the product space of d circles • For n ∈  define the random field It can be written as showing that it is well-defined and Gaussian, since the sum and Riemann integral over a bounded domain of a Gaussian processes with a.s.continuous paths are Gaussian processes again [24] [25].There is a version of j Y with a.s.continuous paths, since its correlation function was assumed to be continuous and integrable such that Corollary 1 can be applied.
• For the trend it follows

Z x Z y A Y x A Y y A x s Y s s A y t Y t t A x s Y s A y t Y t s t A x s Y s Y t A y t s t
using the transformation , u s y v t y = − = − and Equation (2).From that, one can which does exist since all components were assumed to be integrable.According to Lemma 2 the limit ( ) ( ) is well-defined and Gaussian.
• Furthermore, from (3) it follows that Z is stationary and by the convolution theorem the spectrum becomes Modeling and simulation is based on the following corollary.
Corollary 2 If a cross-covariance ρ is previously defined, the matrix A is given by Fourier transform of Â , which can be deduced by ˆˆˆt AA R ρ = using Cholesky decomposition such that Z has cross-covariance ρ .

Simulation
This section is intended to describe and to demonstrate how simulation of the new model works by means of a popular example.The second part evaluates the computational effort in comparison to the spectral turning bands approach.

Simulation Procedure and Example
Based on Theorem 2 and Corollary 2 one can describe a very simple simulation algorithm.Let functions , i j A and R be given.For simulation in a bounded domain do the following.
• Discretize the domain by a grid  .• Apply the inverse discrete Fourier to Ẑ to obtain the field Z.
The algorithm uses the same steps as the usual approach based on Fourier transform for an arbitrary Gaussian field ( [7] [10] [26]) with additional multiplication of ( ) There are numerous models to obtain valid cross-covariance functions.The popular multivariate Matérn model ( [3] [27]) extending the univariate case has only a finite number of real parameters for scale and shape of the field.Define the isotropic Matérn covariance function, that is , .

Let us consider the quadratic domain [ ] [ ]
0,500 0,500 × discretized by a grid of 500 500 × nodes.The model and simulation procedure was implemented using the commercial software Mathematica ® .
The components of ρ are shown in Figure 1.

A realization of the initial univariate random fields 1
Y and 2 Y such as a sample of the resulting cross-correlated bivariate Gaussian field is illustrated in Figure 2.

Computational Effort
This short paragraph compares theoretical and practical computation times of spectral turning bands and spectrum convolution on the same machine for varying numbers of sample locations.
For the spectral turning bands method it can be deduced that the computational cost are of ( ) O LM , where M is the number of locations x to sample and L the number of

O M
. An advantage over circulant-embedding techniques is the possibility to split the area of locations into small subsets allowing for a considerable reduction of memory storage requirements and for parallel computations.
Considering the computational costs for our spectrum convolution approach, since d and N are small numbers, matrix multiplication and Fourier transform are the main operations to be considered.Both can be performed in  and process these subsets consecutively and/or independently.
Considering a practical example performed on a single core 2 ghz laptop cpu.Recall that all code was written in Mathematica software.In Figure 3 the computational time for 10 realizations and different numbers of sample locations M for each approach is shown.Spectrum convolution (blue) is much faster than spectral turning bands, which was computed for L M = (red) and 10 log L M = (green).If L scales logarithmic both implementations have the same asymptotic, which is also verified theoretically, but spectrum convolution is still more than 40 times faster.

Conclusion
As a contribution to the topic of stochastic processes in general and to random fields in particular, a new multivariate modeling approach was presented.It allows modeling and simulation of exact stationary multivariate Gaussian random fields where also the case of isotropy is covered.It is remarkable that any such a random field can be obtained, provided the components of its cross-covariance are continuous and integrable functions.The model is easy to implement since only simulation of normal variables, matrix multiplication, Cholesky decomposition and Fourier transform are required tools.It is shown numerically that the algorithm is more than 40 times faster than spectral turning bands and it can also be modified for parallel computing but, unfortunately, it is limited to sampling in locations on a regular grid.
on the lag r x y = − , where ⋅ denotes the euclidean norm.

Fix
non-negative random reals with density h and , , , l n l n θ φ are mutually independent random variables uniformly distributed over the interval [ [ 0, 2π and independent of the , position n.The deterministic "lower triangular matrix"-valued function A is uniquely given by Cholesky decomposition such that d-torus is obtained the d-dimensional hypercube by gluing the opposite faces together giving a cubic domain with periodic boundaries.Because of periodicity for a continuous function :

Y
≡ for 1 i d ≤ ≤ .• For the cross-covariance function of n Z one can obtain the formula

Ẑ
g by multiplication of the matrix ( ) Â g and the vector  ( )  ( )


Here K ν is a modified Bessel function of the second kind and 0 a > is a spatial scale parameter referred to as a correlation length.The smoothness parameter 0 ν > defines the Hausdorff dimension and the differentiability of the sample paths, since for a positive integer k, the sample paths are k times differentiable if and only if k

Figure 2 .
Figure 2. Top, samples of the univariate Gaussian random fields.Bottom, resulting sample of the bivariate Gaussian field, x-component (left), y (right).
. Thus, the method is of the same scale than circulant-embedding (one of the fastest alternatives to date with ( )log O MM ).Again one can split the area of locations into smaller subsets

Figure 3 .
Figure 3. Computation time as a function of the number of sample locations for spectral turning bands with L M = (red), 10log L M =(green) and spectrum convolution approach (blue). d For each grid point g ∈  deduce the vector