The efficient simulation of isotropic Gaussian random fields on the unit sphere is a task encountered frequently in numerical applications. A fast algorithm based on Markov properties and fast Fourier transforms in 1d is presented that generates samples on an grid in . Furthermore, an efficient method to set up the necessary conditional covariance matrices is derived and simulations demonstrate the performance of the algorithm. An open source implementation of the code has been made available at https://github.com/pec27/smerfs.
Methods for modelling spatially distributed data with isotropic random fields occur in many areas, including astrophysics [2, 7], geophysics , optics , image processing  and computer graphics . The building block for almost all of which is the Gaussian random field in , although many applications require realisations in other geometries, the most important of which may be the unit sphere , especially for geophysics and astrophysics, which often use spherical data.
The computational complexity of realising an lattice of points of a Gaussian random field in depends considerably upon the structure of the covariance function. In the worst (non-isotropic) case, we must calculate the covariance between each point and perform a Cholesky decomposition of the resulting matrix, at a cost of operations . Fortunately, the majority of Gaussian random fields we are interested in are isotropic (i.e., the covariance depends only on the geodesic distance between points), allowing us to perform a spectral decomposition into harmonic (Fourier) components. In the periodic case on (i.e., a torus) with a regular lattice, this provides a dramatic computational improvement due to the existence of the fast Fourier transform (FFT) algorithm, which allows the realisation to be computed in operations.
On , similar results for isotropic Gaussian random fields also apply, i.e., we can perform a spectral decomposition into the spherical harmonic functions. Unfortunately, the corresponding spherical transform algorithm has fewer symmetries to exploit, and for a general grid on the sphere, steps are required for an isotropic covariance function. If one chooses an iso-latitude system such as the equidistant cylindrical projection (i.e., a regular grid in θ and ϕ, the inclination and azimuthal angles) or the HEALPix discretisation (, as used for the cosmic microwave background data of WMAP), one can exploit the Fourier element of the azimuthal angle to reduce this to (see, e.g., ).
As one moves to higher and higher resolutions however, it becomes even more attractive to find algorithms with superior scaling. One such possibility is to consider Gaussian Markov random fields (hereafter GMRFs) (e.g., [23, 10, 16, 21, 13, 22]) which extend the Markov property to higher dimensions and orders. The essential property of these fields is that the probability measure for a conditional realisation of these fields depends upon only the boundary of the constraining set (and its derivatives) rather than the full set. Such a property allows one to construct a realisation recursively, i.e., to iteratively expand the realisation around its boundary, a process sometimes described as a telescoping representation .
The requirement for an isotropic GMRF to obey this property is that the covariance function is the Green’s function of an operator that is a polynomial in the Laplacian (with appropriate coefficients to ensure strong ellipticity; see, e.g., ). This turns out to describe rather a large class of covariance functions .
The contribution of this paper is an algorithm that generates samples on an grid of the sphere in once the conditional covariance matrices are computed. This is achieved by decomposing the field into 1d Gaussian Markov random fields. These can be sampled together with the derivatives point by point and then transformed to an isotropic Gaussian random field on the unit sphere by FFT.
The paper is structured as follows: In Section 2, we derive the decomposition of an isotropic Gaussian random field into 1d GMRFs via Fourier transforms, and compute the conditional covariance matrices. A computationally efficient method to set up the conditional covariance matrices is presented in Section 3. We collect the results of the previous sections in Section 4, where we present the algorithms explicitly. Finally, the performance and convergence of the introduced algorithm is demonstrated in a simulation example in Section 5. Our implementation is available online at https://github.com/pec27/smerfs.
2 Decomposition of isotropic Gaussian random fields into 1d Gaussian Markov random fields
Let us assume that T is a zero mean 2-weakly isotropic Gaussian random field (GRF for short) on the unit sphere in , i.e., on
where denotes the Euclidean norm. Then T admits an expansion with respect to the surface spherical harmonic functions as mappings , which are given by
for , and , and by
for and . Here denote the associated Legendre polynomials which are given by
for , and , where are the Legendre polynomials given by Rodrigues’ formula (see, e.g., )
for all and . Associated Legendre polynomials with negative satisfy
This expansion of T is given by (see, e.g., [8, Corollary 2.5])
where is a sequence of complex-valued, centred, Gaussian random variables with the following properties:
is a sequence of independent, complex-valued Gaussian random variables.
The elements of with satisfy and independent and are distributed.
The elements of with are real-valued and distributed.
The elements of with are deduced from those of by the formulae
Here is called the angular power spectrum.
In what follows, let us consider the case that there exist such that
for all , i.e., is an eigenvalue of
with corresponding eigenfunctions , where denotes the spherical Laplacian (also known as the Laplace–Beltrami operator). In the spirit of , T is the solution of the stochastic partial differential equation
where W is white noise which admits the formal Karhunen–Loève expansion
Here is a set of independent, complex-valued standard normally distributed random variables independent of the real-valued standard normally distributed random variables . For , the same relations as in condition (iv) hold. For fixed , define the random field
for . We observe that
and obtain for any ,
where we set
Moreover, we obtain for the relation
Since is generated by a sum of centred Gaussian random variables, it is clear that is centred Gaussian. In the following lemma, we show that they are independent for positive m and compute the covariance.
The sequence consists of pairwise independent centred Gaussian random fields on with covariance given by
while the functions with negative index are determined by the relation
We have already seen that is centred Gaussian and that . It remains to show independence as well as to compute the covariance. Therefore, let us fix and . By the previous computations it holds that
using the properties of the random variables in , where for and 0 else. Similarly, we obtain
which finishes the proof since uncorrelated Gaussian random variables are independent. ∎
We define the operators
Let us recall the associated Legendre differential equations for and (see, e.g., ):
which is solved by the associated Legendre polynomials. Therefore, we obtain that
i.e., solves the stochastic differential equation
where is white noise on and is a well-defined operator of order M since is of order with positive spectrum.
For , let us first consider the open sets with complements . Since white noise satisfies the strong Markov property, we obtain for all that
has the property. In other words, we obtain for all ,
Therefore, it is sufficient to know for when generating a sample value for any where the random field is already constructed for all . Since we are considering Gaussian fields which have Gaussian derivatives, it is sufficient to know the mean and the covariance between any two points and derivatives up to order .
Since all derivatives are well-defined, we obtain the covariance functions, which we will refer to as cross covariances,
These derivatives can analytically be calculated with the identities of the associated Legendre polynomials. We will discuss a computational efficient method in Section 3. Let us set
such that is the matrix consisting of all covariances between different derivatives at and . The conditional distribution of given is known to be Gaussian with conditional mean
and conditional covariance matrix
By using the values at for non-negative m, it is straightforward to sample the vector from the conditional distribution by
where the random variables are independent complex-valued standard normally distributed for , and independent standard normally distributed for . This construction is sometimes referred to as a state-space model .
Note that we will also need an initial sampling, and for this we choose the equator (), i.e.,
3 Efficient covariance computation
We have shown in Lemma 2.1 that
i.e., we obtained a representation in term of the associated Legendre polynomials. This allows us to compute it but from an algorithmic perspective; the computation is expensive for large m and as the recurrence relations are unsuitable for numerical computation with limited precision. As such we find an alternative formulation in terms of Legendre functions that is numerically superior at even moderate-m values.
Let us assume that the are distinct in equation (2.1) and have the representation with to avoid singularities. More specifically, for , either or is in , and therefore is not well-defined. A typical example is and with , which satisfies the desired partial fraction decomposition with and . Then the coefficients in equation (2.1) can be decomposed (via partial fractions) into
for some finite constants , provided the are distinct.
We note that the partial sums in equation (3.1) can be written in the form
which is nothing more than the bilinear expansion for the Green’s function [17, (17.72)] for the operator . This operator can be alternatively formulated as the solution to the differential equation
One eigenfunction of this operator with eigenvalue zero is the associated Legendre function , where denotes the generalisation of to non-integers .
Associated Legendre polynomials can be generalised to complex degree λ and order μ as
via the (Gauss) hypergeometric function
where denotes the rising Pochhammer symbol. The relevant case of μ being an integer (i.e., ) is only defined in the limit, which we find by (see [1, 15.3.3])
and (see [1, 15.1.2])
This function is well-defined at with , and singular at for all if .
For a second linearly independent solution to (3.3) we use , which is not linearly dependent since is not an integer; see, e.g., [1, 8.2.3]. We can then construct the Green’s function by using the Wronskian (e.g., [15, 14.2.3])
which satisfies continuity, regularity at and , and by substitution,
In conclusion, we obtain that
Substituting the associated Legendre functions by the hypergeometric functions as in (3.4), i.e., setting
and similarly for q and , we derive our efficient computation method for the cross-covariance matrices .
We observe that the above formula requires many re-computations of associated Legendre functions. To be more efficient, we transform the vector of derivatives up to order to a vector of linear combinations of derivatives of increasing order. This does not change anything since we are in the end only interested in the random field itself but not in its derivatives. Therefore, define for the family of operators by
which satisfies . Thus setting
is much easier to compute while the properties of the random field itself are not changed. Let us keep the notation even if we use the modified matrices .
Below we outline an algorithm for the construction of the isotropic GRF on . Since the sampling of the GRF is merely filtered white noise, in our tests the time is dominated by the computation of the (inhomogeneous) filter coefficients, and as such we have described the pre-computation in Algorithm 1, which just has to be done once, and the generation of samples in Algorithm 2.
For , let us fix a grid on and for some . One possible symmetric example is to set such that
and in . Furthermore, define the covariance operator by setting and in (3.2).
Algorithm 1 (Pre-computation of covariance matrices.).
Algorithm 1 is devoted to the pre-computations that set up the covariance matrices that just have to be computed once and can be reused for sampling on the same grid. We remark that the algorithm can be optimised if the discrete grid z is symmetric around zero, i.e., if it satisfies . Then and just have to be computed for all and the negative values follow by symmetry. This is the case for the example grid z above.
Samples of the GRF on the sphere are generated with Algorithm 2 based on the pre-computed families of matrices A, B and . Since the computational cost to generate white noise is negligible, the complexity of the algorithm reduces to FFT in 1d, which makes it computationally fast. Therefore, the algorithm is attractive if many samples have to be computed. An example are ice crystals, which can be modelled by lognormal random fields . For more information on the relation of lognormal and Gaussian random fields, the reader is referred to .
Algorithm 2 (Gaussian random field generation.).
To illustrate the algorithm and show the performance of the presented algorithm, we show simulation results in this section which include the generation of a sample and the convergence of the covariance function with respect to the discretisation.
Let us first show random field samples that the algorithm generates. Therefore, we fix the angular power spectrum given by
The algorithm in Section 4 generates the one-dimensional random fields on a discrete grid first, which are then transformed to random fields on the sphere by FFT. To give the reader an idea of , we show in Figure 1 samples of and the real part of along with its standard deviation or RMS expectation .
The standard deviation is dependent on θ since the random field is of zero mean but not translation invariant.
Note that tends to zero at . Since T has a continuous first derivative (i.e., , as it is in this case), then for all the standard deviation is zero for , i.e., the nonzero iso-latitude Fourier transforms must tend to zero near the poles to meet the smoothness criterion.
The process of generating a sample is illustrated in Figure 2. We start with generating the random fields , which become after FFT the discrete sample of the random field at the equator. Conditionally, and are sampled and added to the picture after FFT. This process is continued until all of the sphere is covered with random numbers and the random field on the sphere is complete.
A method to validate that a method for sampling a GRF is appropriate and that it converges is to compute the covariance function via the examination of many samples and compare it to the analytic one inferred from the angular power spectrum. We have used the GRF with angular power spectrum from equation (5.1) and corresponding covariance function
which is plotted in the left two panels of Figure 3 as solid blue line.
In order to analyse the errors of our method, we have constructed filters at multiple resolutions of , keeping and (i.e., so the maximum angular separation between adjacent grid points is ). We generated samples on all resolutions to estimate the covariance by
and to compare it to the theoretical one. The dotted and dashed lines in the left two panels of Figure 3 show the results for the different . The left figure is based on equatorial evaluations, while the middle one estimates the covariance along meridional lines.
The error was computed by taking on each grid the maximum over all grid points in a set , which was once the equator and once the meridian, of the difference of the theoretical and statistical covariance, i.e., the error was computed by
The results are shown in the right panel of Figure 3. The maximum error on the covariance function falls as (where ) as can be seen in the error plot.
Funding source: Knut och Alice Wallenbergs Stiftelse
Award Identifier / Grant number: 621-2014-3995
Funding source: Vetenskapsrådet
Award Identifier / Grant number: 621-2014-3995
Funding statement: The work of the second author was supported in part by the Knut and Alice Wallenberg foundation and the Swedish Research Council under Reg. No. 621-2014-3995.
The first author would like to thank Ben Lowing, and the second author would like to thank the members of the KAW project “Stochastics for big data and big systems” for fruitful discussions.
 M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, Dover Publications, Mineola, 1965. Search in Google Scholar
 W. L. Brogan, Modern Control Theory, 3rd ed., Prentice-Hall, Upper Saddle River, 1991. Search in Google Scholar
 F. Cohen, Z. Fan and M. Patel, Classification of rotated and scaled textured images using Gaussian Markov random field models, IEEE Trans. Pattern Anal. Mach. Intelligence 13 (1991), no. 2, 192–202. 10.1109/34.67648 Search in Google Scholar
 K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke and M. Bartelmann, HEALPix: A framework for high-resolution discretization and fast analysis of data distributed on the sphere, Astrophys. J. 622 (2005), 759–771. 10.1086/427976 Search in Google Scholar
 Y. Hoffman and E. Ribak, Constrained realizations of Gaussian fields – A simple algorithm, Astrophys. J. Letters 380 (1991), L5–L8. Search in Google Scholar
 A. Lang and C. Schwab, Isotropic Gaussian random fields on the sphere: Regularity, fast simulation and stochastic partial differential equations, Ann. Appl. Probab. 25 (2015), no. 6, 3047–3094. 10.1214/14-AAP1067 Search in Google Scholar
 F. Lindgren, H. Rue and J. Lindström, An explicit link between Gaussian fields and Gaussian Markov random fields: The stochastic partial differential equation approach, J. R. Stat. Soc. Ser. B. Stat. Methodol. 73 (2011), no. 4, 423–498. 10.1111/j.1467-9868.2011.00777.x Search in Google Scholar
 E. Mendez and K. O’Donnell, Observation of depolarization and backscattering enhancement in light scattering from Gaussian random surfaces, Optics Commun. 61 (1987), no. 2, 91–95. 10.1016/0030-4018(87)90225-2 Search in Google Scholar
 T. Nousiainen and G. M. McFarquhar, Light scattering by quasi-spherical ice crystals, J. Atmospheric Sci. 61 (2004), no. 18, 2229–2248. 10.1175/1520-0469(2004)061<2229:LSBQIC>2.0.CO;2 Search in Google Scholar
 F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller and B. V. Saunders, NIST Digital Library of Mathematical Functions, http://dlmf.nist.gov/, Release 1.0.15 of 2017-06-01. Search in Google Scholar
 K. F. Riley, M. P. Hobson and S. J. Bence, Mathematical Methods for Physics and Engineering: A Comprehensive Guide, 2nd ed., Cambridge University Press, Cambridge, 2002. Search in Google Scholar
 Y. A. Rozanov, Markov Random Fields, Springer, New York, 1982. Search in Google Scholar
 H. Rue and L. Held, Gaussian Markov Random Fields: Theory and Applications, Monogr. Statist. Appl. Probab. 104, Chapman & Hall, London, 2005. Search in Google Scholar
 G. Szegő, Orthogonal Polynomials, 4th ed., Amer. Math. Soc. Colloq. Publ. 23, American Mathematical Society, Providence, 1975. Search in Google Scholar
 A. Tewfik, B. Levy and A. Willsky, Internal models and recursive estimation for 2-d isotropic random fields, IEEE Trans. Inform. Theory 37 (1991), no. 4, 1055–1066. 10.1109/18.86997 Search in Google Scholar
 D. Vats and J. Moura, Telescoping recursive representations and estimation of Gauss–Markov random fields, IEEE Trans. Inform. Theory 57 (2011), no. 3, 1645–1663. 10.1109/TIT.2011.2104612 Search in Google Scholar
© 2018 Walter de Gruyter GmbH, Berlin/Boston
This work is licensed under the Creative Commons Attribution 3.0 Public License.