On Simulation of Manifold Indexed Fractional Gaussian Fields

To simulate fractional Brownian motion indexed by a manifold poses serious numerical problems: storage, computing time and choice of an appropriate grid. We propose an eﬀective and fast method, valid not only for fractional Brownian ﬁelds indexed by a manifold, but for any Gaussian ﬁelds indexed by a manifold. The performance of our method is illustrated with diﬀerent manifolds (sphere, hyperboloid).


Introduction
Rough phenomena arise in various fields (Frisch and Parisi 1985;Leland et al. 1994;Mandelbrot 1975;Peitgen and Saupe 1988;Pentland 1984): texture simulations and image processing, natural scenes (clouds, mountains) simulations, fluid mechanics, financial mathematics, ethernet traffic. . . Some phenomena, like ethernet traffic or financial data, are time-indexed. Other phenomena, like image processing, should be indexed by subsets of the Euclidean spaces R 2 or R 3 . But, some other rough phenomena are not indexed by an Euclidean space, but by a manifold. Let us mention for instance the cosmic microwave background (Marinucci et al. 2007; National Aeronautics and Space Administration 2010) or solar data (Koenig and Chainais 2008), that lead to spherical data. Real-world textures are usually not supported by an Euclidean space, but by a surface. Applying an Euclidean indexed texture on, say, a sphere, always generates artefacts. Generating manifold indexed textures is therefore a challenge.
Fractional Gaussian fields, like fractional Brownian fields, are good candidates for modeling rough phenomena. Fractional Gaussian fields may be indexed by an Euclidean space or a manifold. The simulation of Euclidean indexed Gaussian fields poses a lot of problem, in particular storage and complexity. Recently, Brouste et al. (2007) propose a new and fast algorithm for simulating Euclidean indexed Gaussian fields. The aim of this paper is to extend Conversely, for any n.n.d. function R(·, ·), there exists an unique centered Gaussian field of second order structure given by R(·, ·).

The function fieldsim
We give here a summary of the function fieldsim that allows to simulate an Euclidean random field. In this case M is a subspace like [0, 1] 2 for instance with the usual Euclidean norm. This function yields discretization of sample path of the Gaussian field over a space discretization {S e , S r } of M, associated with any n.n.d. function R(·, ·). It is consisted of the two following steps.
Accurate simulation step. Given a space discretization S e , a sample of a centered Gaussian vector: (X(M )) M ∈Se with covariance matrix R given by R i,j = R(M, M ), M, M ∈ S e , is simulated. This simulation is obtained by an algorithm based on Cholesky decomposition of the matrix R.
Refined simulation step. Let S r be the remaining space discretization. For each new point M ∈ S r at which we want to simulate the field, X(M ) is generated by using only a set of neighbors instead of all the simulated components (as in the accurate simulation step). Precisely, let N M be a neighbors set of M (for the Euclidean distance in R 2 ) and X N M be the space generated by the variables X(M ), M ∈ N M . Let us remark that the neighbors set is defined with all the already simulated variables (in the accurate and refined simulation step). Let X X N M (M ) be the best linear combination of variables of X N M approximating X(M ) in the sense that the variance of the innovation ε where U is a centered and reduced Gaussian variable independent of the already simulated components.
Note that the variable X X N M (M ) and the variance V ar(ε X N M (M )) are completely determined by the covariance structure of the sequence X(M ), X(M ), M ∈ N M . For storage and computing time, the accurate simulation step must concern only a small variables number whereas the second step can relate a larger variables number. That leads to an effective and fast method to simulate any Gaussian field.
In Brouste et al. (2007), the procedure is implemented in the R package FieldSim. A natural discretization space is of the form (k2 −J , l2 −J ), k, l = 0, . . . , 2 J for J a positive integer. The accurate step is applied for points (k2 −Ja , l2 −Ja ), where J a is a level to be chosen (in general J a = 1 or 2). The refined step is applied for the remainder.

The function fieldsim adapted to manifolds
The previous procedure can be adapted to the case of fields indexed by a manifold. The main problem stands in the discretization grid choice. There is, in the case of the sphere for instance, no equidistributed grid and it is difficult to define a concept of finer grid such as in the case of field indexed by [0, 1] 2 . Moreover, this choice can be related to the software that one wishes to use to represent the manifold.
We propose to develop a function for which the user enters a grid of his/her choice for each step (accurate and refined) and a distance to be used in order to determine the neighbors.
Thus one can use our function for any manifold when having a covariance function, a grid and a distance. We give here two examples of manifolds: the sphere and the hyperboloid.
Let us precise how we use it for this two manifolds. We denote by S g the set of the point of M at which we want to generate the process (the visualization grid). This set choice can be induced for instance by the software that one wishes to use to represent the manifold (see Appendix C). Let us recall that we want to generate a path of field indexed by the manifold M (discretized at S g ) and with covariance function R(·, ·). We need also to specify the concept of neighbors. A natural choice is to use the geodesic distance between two points. So the closest neighbors of some point M are the points closest to M according to this distance. In general, the cardinal of S g is too large to use only the accurate simulation step. Moreover contrary to the grid chosen in the case of field indexed by [0, 1] 2 , one is not able in general to use any more a concept of finer grid. Indeed, for our sphere visualization grid for instance, using overlapping sub-grids (an atlas of 6 maps), each sub-grid have some parts of the sphere with very few points and there are some others with points accumulation (in particular on the six poles). That is why in this case, we propose to first simulate the process at uniform random points and next to simulate the process at the S g points. The generation of uniform random points on the manifold as sphere or hyperboloid are given in Appendix B. Precisely, let S u the uniform random points set. We propose to run the function fieldsim with the set

How to chose the parameters?
This function requires in addition to R(·, ·), the integers N e , the points number for the accurate step, and nbNeighbor, the neighbors number. Choosing these parameters is a delicate task. This problem can be related to the question of evaluating the simulation accuracy. One needs to define an appropriate accuracy criterion. One would wish to evaluate the "fractionality" of the simulated field. That is what was proposed in Brouste et al. (2007) by estimating the fractional index. This is only partially satisfactory. A second, and more classical, way is to estimate the mean square error between the simulated sample path and the true one. A third one consists in evaluating a distance between the covariance of the simulated field and the true covariance. Second, one needs the effective computation of the error. We have investigated the mean square error in a theoretical way. Apart obvious results (i.e. the number of neighbors used in the refined grid tends to infinity), the computations are so intricate that we did not succeed. The same occurs for studying theoretically a distance between the covariance of the simulated field and the true covariance.
In this subsection, we propose two protocols to help the user in the error control and the parameters choice. To do that, we propose to estimate the distance between the covariance of the simulated field and the true covariance since in our context this function is known. Firstly, one can simulate several paths and estimate this error on the part of the process simulated in an exact way and that of the process simulated in an refined way. The first error only corresponds to the sampling error. The second one contains in more the approximation error due to the use of some neighbors instead of all the already simulated variables. Thus for a covariance function and a given grid, one compares these two kinds of error. In particular one uses this protocol to chose the neighbors number that results from a compromise between computational times and error loss control. Table 1 gives the results for different values of number of neighbors obtained for n s = 10, 000 sample paths simulated over the same  100 uniform random points. The points number for the accurate step is N e = 50, so the points number for the refined step is N r = 50. The covariance function is given by R S 1 (·, ·) defined at the Section 3.1 (with H = 0.25). The distance between the covariance of the simulated field and the true covariance is estimated by the mean 2 -distance between the estimated covariance function coefficients and the true ones. Firstly we remark that the ratio of the errors (accurate and refined) goes to 1 as the neighbors number increases. For this example, one also observes that there is an important gain while passing from 8 to 12 neighbors. But it seems unnecessary to increase this number of neighbors to some value greater than 15.
Secondly, concerning the choice of the size of the random grid with respect to the size of S g , one uses the following protocol. For several grid size values (here N e fixed to 50), one estimates the distance between the covariance of the simulated field and the true covariance (over several paths). Next one retains the value of the grid size for which this error is smallest. Table 2 gives the results for different sizes of random grid S u obtained for n s = 10, 000 sample paths simulated. The size of the visualization grid is equal to 192. For this example, we observe that we do not need to increase the size or the random grid beyond 81.

Numerical results
In this section, we illustrate the method proposed here through simulations for two manifolds: sphere and hyperboloid.

Some examples of sphere indexed fractional fields
Let S be the unit sphere of R 3 , and let d S be its geodesic distance. We give hereafter several examples of fields indexed by S. The spherical fractional Brownian fields were introduced by Istas (2005) in the following way. There exists a centered Gaussian field called spherical fractional Brownian field (sfBf) whose covariance function is given by where O is any given point of S, if and only if H ∈ (0, 1/2]. We can show (see Istas (2009)) that the following functions  Finally we can restrict the fractional Brownian fields indexed by R 3 to the sphere. We obtain the following covariance function: where H ∈ (0, 1).
Figures 1 and 2 gives some representations of such fields. For each sphere, we have chosen N e = 100, N r = 900 and nbNeighbor = 15. The discretized grid S g is given as in Appendix C with N g = 100.

Some examples of hyperboloid indexed fractional fields
Let H be the hyperboloid: Since H is unbounded, we will simulate the field on a hyperbolic cap. We give hereafter several examples of fields indexed by H.

The hyperbolic fractional Brownian fields (hfBf) introduced by Istas (2005) have covariance function given by
where O is any given point of H and H ∈ (0, 1/2]. We can show (see Istas (2009)) that the following functions , are covariance functions for H ∈ (0, 1/2]. Figure 3 gives some representations of such fields. For each example, we have chosen N e = 100, N r = 900 and nbNeighbor = 15. The discretized grid S g is given as in Appendix C with N g = 150.

A. Using FieldSim
The new version of FieldSim is a set of R functions that allows performing simulations of manifold indexed fractional Gaussian fields with known covariance function. Three classes of functions are implemented: Three building functions: setManifold permits the user to create an object of class manifold; constructcovf that produces covariance functions; constructgrid that produces several grids for different manifolds.
Simulation function fieldsim that performs simulations of the path of the manifold indexed Gaussian fields.
Print function visualize that gives graphical representation of path of some manifold indexed Gaussian fields.
The R environment is the only user interface returned to R. In order to make it easier for the reader not familiar to R language, we detail the call to functions and the command used to produce graphical outputs.

A.1. Setting manifold
A manifold is an object of the S4 class manifold with four slots which are: name which is the name of the manifold we consider (character); atlas which the union of discretized domains that covers the manifold (must be a matrix where the number of line is the dimension of the space where the manifold lives); distance which is the distance considered on the manifold (function); origin which is the origin considered on the manifold (must be a point on the manifold).
The setter setManifold permits the user to create an object of class manifold with all its slots. We give below the example of the field indexed by some interval.

A.2. Building grids
Some grids can be constructed with the function constructgrid(manifold, typegrid, Ng). It has, three argument, the manifold on which the grid is settled, the type of grid the user want to construct and the number of points of this grid. There are four types of grid: regular, finer, random and visualization. Some informations on the random grid and on the visualization grid for the sphere and the hyperboloid manifolds are postponed in Appendixs B and C.

A.3. Building covariance function
In order to simulate further Gaussian processes or fields, we have to define a covariance function. In fact the user can construct himself this covariance function (semi-definite positive on the manifold). In order to treat more easily some particular case, a function have been implemented that constructs for known manifolds some particular covariance functions. Let us examine the call of this function constructcovf(manifold, "fBm", H = 0.6). The first argument is a manifold object, the second parameter is the type of covariance (in this case, fBm covariance function) and the third argument (depending of the type) is parameter of this covariance (for the fBm, it is the Hurst parameter H). For the moment, the user can choose between two types, namely: "fBm" (fractional Brownian field) and "mBm" (multifractional Brownian field).

A.4. Sphere indexed fractional fields
The function fieldsim takes as arguments a manifold, and a covariance function. It returns the values of the field on the atlas of the manifold computed by our method. In this case, the slots atlas, distance and origin are fixed to default values, namely a visualization grid, the usual geodesic distance on the sphere R> dS <-function(xi, xj){ + u <-sum(xi * xj) + if (u < -1) u <--1 + if (u > 1) u <-1 + return(acos(u)) + } and an origin fixed to rbind(1, 0, 0). As explained in Appendix C, it is natural that NaN values exist in the sphere default visualization atlas.

Defining the covariance function
In order to compute sample paths of the spherical fractional Brownian field, it is possible to use the function R> H <-0.4 R> R.S.1 <-constructcovf(sphere, "fBm", H = H) to compute its covariance function. Other type of covariance functions have been defined in Section 3.1 which use the manifold distance. For instance, let us construct the second example R> dS <-sphere@distance R> R.S.2 <-function(xi, xj) exp(-dS(xi, xj)^(2 * H))

Simulating via the FieldSim method
The commands (the first one to integrate the simulation grid to the manifold object and the second to compute the sample path) R> sphere@atlas <-simulationgrid R> resS <-fieldsim(sphere, R.S.1, Ne = 100, nbNeighbor = 15) simulate a spherical fBm on the mesh given by the concatenation of the random grid and the visualization grid. The quantity Ne is the number of points of the grid to be simulated in the accurate step and nbNeighbor the number of neighbors used in the refined step.

Visualizing the sphere indexed fractional field
Note that the function visualize needs the packages rgl (Adler and Murdoch 2010) and RColorBrewer (Neuwirth 2007) that load automatically with the FieldSim package. Before using the function visualize, the user has to defined the visualization atlas of the manifold object and the proper sample path associated with this atlas R> sphere@atlas <-S.g R> res <-resS [(dim(S.u)[2] + 1):length(resS)] After that, the command R> visualize(sphere, res) allows us to visualize the sphere indexed fractional field.

A.5. Hyperboloid indexed fractional fields
To obtain hyperboloid indexed fractional fields, similar functions can be used by starting with R> hyper <-setManifold("hyperboloid") B. Uniform random generator on manifold B.1. Spherical uniform random generator We parameterize the unit sphere by the polar coordinate system x = cos θ sin φ, y = sin θ sin φ, z = cos φ, where φ ∈ [0, π] and θ ∈ [0, 2π]. To generate uniform random points on the unit sphere, we can use the following way. Firstly we generate θ with an uniform distribution on [0, 2π]. Secondly (and independently of the first sample), we generate φ as a random variable on [0, π] with density sin(φ)/2.

B.2. Hyperboloid uniform random generator
We parameterize the hyperboloid by the hyperbolic coordinate system x = cos θ sinh φ, y = sin θ sinh φ, z = cosh φ, where φ ≥ 0 and θ ∈ [0, 2π]. Let us remark that there does not exist uniform random generator on hyperboloid since this space is not bounded. So one considers the hyperbolic cap defined by {x 2 + y 2 − z 2 = −1, 1 ≤ z ≤ Z}, where Z is some real in (1, ∞). To generate uniform random points on the hyperbolic cap, we can use the following way. Firstly we generate θ with an uniform distribution on [0, 2π]. Secondly (and independently of the first sample), we generate φ as a random variable on [0, arccoshZ] with density sinh(φ/(Z − 1)).