Conditional fractional Gaussian ﬁelds with the package FieldSim

We propose an effective and fast method to simulate multidimensional conditional fractional Gaussian ﬁelds with the package FieldSim . Our method is valid not only for conditional simulations associated to fractional Brownian ﬁelds, but to any Gaussian ﬁeld and on any (non regular) grid of points.


Introduction
Rough phenomena arise in texture simulations for image processing or medical imaging, natural scenes simulations (clouds, mountains) and geophysical morphology modeling, financial mathematics, ethernet traffic, etc. Some are time-indexed, some others, like texture or natural scene simulations, should be indexed by subsets of the Euclidean spaces R 2 or R 3 . Recent data (as the Cosmic Microwave Background or solar data) are even indexed by a manifold.
The fractional Brownian motion (fBm), introduced by Kolmogorov (1940) (and developed by Mandelbrot and Van Ness (1968)) is nowadays widely used to model this roughness. Fractional Brownian motions have been extended in many directions: higher dimensions with fields, anisotropy, multifractionnality, etc.. This paper is devoted to a simulation method for conditional Gaussian fields. It could improve, in the future, natural scene simulations by fixing for instance the valleys.
The simulation of fractional Gaussian processes is not difficult in dimension one (see a review of Coeurjolly (2000)). Let us recall the numerical complexity of some classical methods: the Cholesky method has a complexity of O(N 3 ) where N is the size of the simulated sample path. For specific stationary processes (on a regular grid) the Levinson's algorithm has a complexity of O(N 2 log N) and the Wood and Chan algorithm (see Wood and Chan (1994)) a complexity of O(N log N).
In higher dimension, the Wood and Chan method has been extended to stationary increments fields with the Stein's method (Stein, 2002) ; the fractional Brownian field can therefore be simulated on a regular grid of the plane. For general Gaussian fields on general discrete grid, the Cholesky method is costly and exact simulations are no longer tractable. Approximated methods have been intensively developed (midpoint (Peitgen and Saupe, 1988), turning bands (Yin, 1996), truncated wavelet decomposition) but for specific fields. On manifolds, simulation procedures based on truncated series of eigenfunctions of the Laplace-Beltrami operator are discussed in Gelbaum and Titus (2014).
Our approach, presented in (Brouste et al., 2007(Brouste et al., , 2010, is based on a 2-steps method with an exact simulation step plus a refined fast step, that is an improvement of the midpoint method. It has been implemented in the FieldSim package. The fieldsim simulation method can be applied to general Gaussian processes on general simulation grids (regular and non regular) on Euclidian spaces and even on some manifolds (see Figure 1). It is worth mentioning that the other package RandomFields allows the simulation of a large class of random fields such as Gaussian random fields, Poisson fields, binary fields, Chi2 fields, t fields and max-stable fields (see Schlather et al. (2015)). These simulations are obtained essentially by construction and not from the covariance structure of the process. In RandomFields, conditional random fields (what is the purpose of the present paper) are given for specific spatial and spatio-temporal random fields. It may be noted that the FieldSim package can simulate any conditional Gaussian process as long as you know its covariance function.
We propose here to adapt the FieldSim package to conditional simulations. Definitions and notation will be introduced in the following section with the process class, the setProcess procedure and the fieldsim procedure. The fieldsim procedure adapted to conditional Gaussian fields is described in the next section. Simulations with the package FieldSim are presented in the last section.

Notation and preliminaries Fractional Gaussian fields
Let d be a positive integer and X(·) = X(M), M ∈ R d be a real valued non stationary field with zero mean and second order moments. It is worth emphasizing that we consider in this paper the metric space R d with the Euclidian norm but it can be generalized to smooth and complete Riemannian manifold equipped with its geodesic distance Brouste et al. (2010). The covariance function R(·, ·) is defined by: This function is nonnegative definite (n.n.d.). Conversely, for any n.n.d. function R(·, ·), there exists an unique centered Gaussian field of second order structure given by R(·, ·).
Different classical fractional Gaussian fields have been simulated to illustrate the FieldSim package in Brouste et al. (2007Brouste et al. ( , 2010. In the sequel, M and M are two points of R d and · is the usual norm on R d , d = 1, 2. We can cite 1. The standard fractional Brownian fields are defined through their covariance function (e.g. Samorodnitsky and Taqqu (1994)): where the Hurst parameter H is real in (0, 1).
2. The standard multifractional Brownian fields are defined through their covariance function (see Peltier and Levy-Véhel (1996); Benassi et al. (1997)): and the Hurst parameter is a continuous function H : R d −→ (0, 1), where Γ is the usual Gamma function.

The FieldSim package
In the new version 3.2 of the package FieldSim, new features have been added. The most important add is the process class and the setProcess function.
An object of class process has different slots: • the name of the process. Several names are reserved for classical fractional Gaussian processes: see Table 1 for details. "cond" is used for all kind of conditional simulations (see further).
• the slot values stores the values of the process on the simulation (and visualisation) grid.
• an object of class manifold which is the Riemannian manifold on which the process is lying ; an object of the class manifold has four slots: name which is the name of the manifold we consider. The name "line", "plane", "sphere" and "hyperboloid" are taken for the eponymous manifolds.
atlas which is the union of discretized domains that cover the manifold (must be a matrix where the number of rows is the dimension of the space where the manifold lives); distance which is the distance considered on the manifold ; 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. This class is already described in Brouste et al. (2010).
• the slot covf which contains the covariance function of the Gaussian process ; • the slot parameter which contains all the parameter s associated to the covariance function of the process. Here are the classical parameters associated to classical process.
All the examples presented can be defined with the setProcess command (see Table 1). With the following command, the user can set a fBm with Hurst parameter 0.7 on a regular grid of the interval It is worth mentioning that the slot values is empty since there is no simulation done. Then as usual, the user can use the fieldsim function in order to simulate the Gaussian process associated to covf on the manifold grid defined in manifold.

R> fieldsim(linefBm)
In the fieldsim function, we can add the quantity Ne, the number of points of the grid to be simulated in the exact step, and nbNeighbor, the number of neighbors used in the refined step. By default, Ne is equal to the size of the grid given in atlas. The slot values are now set with the simulated values. There exist different visualisation procedures to draw the results, for instance R> plot(linefbm,"default") We recall that the discretization grids can be modified with the setAtlas command. Depending on the manifold, there is several type of grid: "regular", "random" and "visualization". For instance, R> setAtlas(linefBm,"regular",1 ) R> fieldsim(linefBm) R> plot(linefBm,"default")

The fieldsim procedure for conditional Gaussian fields
In order to build conditional fractional Gaussian fields, we consider a conditioning set N = {N 1 , . . . , N k } , N i ∈ R d , i = 1, . . . , k, and the conditioning values x = (x 1 , . . . , x k ) T ∈ R k . Then we will say that X(·) = { X(M), M ∈ R d } is the conditional Gaussian field associated to the field X(·) (of covariance function R) and to the conditioning pair (N , x) if the finite dimensional laws of X(·) is the same as the finite dimensional laws of X(·) given the event {(X (N 1 ), . . . , X(N k )) T =: X N = x}. We denote by m(·) (resp. R(·, ·)) the mean (resp. covariance) function of the process X(·). The following lemma allows us to determine m(·) and R(·, ·) according to R(·, ·) (sketch of proof is given in (Piterbag, 1996, Section A.1)).

Lemma 1. Let us consider the centered Gaussian vector
and covariance matrix In the Gaussian field context, Lemma 1 allows us to write down an explicit expression of the mean function and the autocovariance function of the conditional Gaussian field associated to R(·, ·) and to (N , x). Let us put Y 1 = X(M 1 ) and Y 2 = X(M 2 ) the values of the field X(·) at points M 1 ∈ R d and M 2 ∈ R d respectively, and Z = X N ∈ R k . Therefore, all quantities in (1) and (2) can be expressed in term of the autocovariance function R. Precisely, Consequently, the mean function of the conditional Gaussian field is given by Then the autocovariance function of a conditional Gaussian field (using (1,2)-coordinate of Equation (2)) is given by For instance, for k = 1, we get Let us recall that the goal of this paper is to give a procedure that yields discretization of sample path of the conditional Gaussian field over a space discretization {S e , S r } of R d associated to the n.n.d. autocovariance function R and the conditioning set and values (N , x) . In the sequel, we denote by X(·) this sample path. Since the mean function (3) is known, we can consider the centered field X(·) = X(·) − m(·).The fieldsim procedure for conditional Gaussian fields takes place as follows.
Exact simulation step. Given a space discretization S e , a sample of a centered Gaussian vector (X(M)) M∈S e with covariance matrix R given by { R} i,j = R(M i , M j ), M i , M j ∈ S e , is simulated. Here R is defined by (4). 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 O M be a neighbors set of M (for the Euclidean distance) and X O M be the space generated by the variables X(M ), M ∈ O 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 O M (M) be the best linear combination of variables of X O 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 For storage and computing time, the accurate simulation step must concern only a small number of variables whereas the second step can relate a larger number of variables. That leads to an effective and fast method to simulate any Gaussian field.
It is worth mentioning that the setProcess command will check if {E(X N X T N )} −1 exists for common conditional simulations.

Some examples of conditional fractional Gaussian fields
We focus, in this paper, on the conditional Gaussian fields associated to the previously mentioned fields but every other classical Gaussian fields can be also simulated: standard bifractional Brownian motion, space-time deformed fractional Brownian motion, etc. (see Brouste et al. (2007)). We also consider conditional simulations associated to fractional Gaussian fields on manifolds (hyperboloid and sphere) (see Brouste et al. (2010) for covariance function definition).
The procedure fieldsim is extended to the conditional Gaussian fields. We can find in the Table 1, the setProcess reference short-card.

On the line
The fractional Gaussian processes on the line are fast to simulate.
They can be obtained with the fieldsim procedure. For instance, the mBm at Figure 2 is obtained with: R> funcH<-function(x){ .3+x* .6} R> cond.mBm<-setProcess("cond-mBm-line",list(Gamma=matrix(c(1/2,1,3/4, .5,1, ),2,3), par=funcH)) R> fieldsim(cond.mBm) R> plot(cond.mBm) In the simulation below, the points of the set N belong to the visualization grid. When this is not the case, the plot could show a failure for the conditioning in the region of high variability. To avoid this, it is possible to add the points of the set N to the visualization grid. For instance, in the previous exemple, to add the point 1/6 to the visualization grid, we can use the following lines of code: R> atlas.cond.mBm <-sort(c(cond.mBm@manifold@atlas[1,],1/6)) R> cond.mBm@manifold@atlas<-matrix(atlas.mBm,nrow=1) An other solution is to use finer grids which contains the points of the set N .

Figure 3:
Conditional simulations associated to fractional Brownian field (on the left) and multifractional Brownian field (on the right). The real time (resp. CPU time) in seconds is equal to 1270.911 (resp. 6.769) for the fractional Brownian field and 1782.533 (resp. 9.953) for the multifractional Brownian field.

On the hyperboloid and on the sphere
Conditional simulations can be extended to fractional Gaussian fields on manifold associated to the fractional Brownian field on the hyperboloid with H = 0.7, The two processes are simulated on a regular grid of 5400 points of R 3 with 100 points for exact simulation step and 5300 for refined step (with 4 neighbors).

Conclusion and perpsectives
We propose a generic method to simulate multidimensional conditional fractional Gaussian fields.
Our method is valid for any Gaussian field and on any (non regular) grid of points as soon as the covariance function is available. This method is constructed to be universal (conditional simulation, simulation on a manifold) and is consecutively not as fast as other methods defined for specific fields. In the near future, the Fieldsim should also possess such specific methods.
Our method is adapted to conditional simulations and, consequently, permits now to simulate easily several natural scenes (clouds, mountains) with valleys and fixed topographic points. Such a simulation is presented in Figure 6.  The process class. "fBm" for fractional Brownian motion, "mBm" for multifractional Brownian motion, "2pfBm" for the standard bi-fractional Brownian motion, "stdfBm" for the space-time deformed fractional Brownian motion, "AfBf" for anisotropic fractional Brownian field and "fBs" for fractional Brownian sheet.