The "SPDE approach", as popularized by lindgren2011explicit, consists in characterizing stationary continuous Markov random fields on (
) as solutions of stochastic partial differential equations (SPDE). This approach has two benefits:
Discrete approximations of the solutions of these SPDEs, obtained using Galerkin methods such as the Finite Element method, are used in place of the original field in numerical computations. lindgren2011explicit actually derived the expression of the precision matrix of the weights of the discrete representation of the solution, thus facilitating the use of this approach.
By tinkering with these same SPDEs, generalizations of stationary continuous Markov random fields on can be defined on manifolds, and oscillating and even non-stationary random fields can be produced (lindgren2011explicit; fuglstad2015exploring).
This work aims at generalizing the SPDE approach to fields that are not continuously Markovian while still keeping the benefits mentioned above. First, the motivations for this work are detailed in order to point out the type of random fields that will be used throughout the developments. Then, an explicit formula for the covariance matrix of the weights of the finite element representation of such fields is provided, and an error analysis is carried out. Finally, an algorithm based on a Chebyshev polynomial approximation and allowing to compute simulations of these weights with a linear computational complexity is introduced.
Note : This paper is a stub presenting the main results obtained by the authors. It is planned to be expanded.
the spatial Gaussian white noise on
defined on a complete probability space. It can be seen as a Gaussian random measure satisfying:
where is the collection of bounded Borel sets of and Leb denotes the Lebesgue measure.
2.1 Solutions of stochastic partial differential equations
Let be a continuous, polynomially bounded function such that:
And let be the pseudo-differential operator defined over sufficiently regular functions of by :
Consider then the stochastic partial differential equation (SPDE) defined over by (vergara2018general):
where is a spatial Gaussian white noise and the equality is understood in the second order sense, i.e. both sides have same (generalized) covariance. The existence and uniqueness of a stationary solution of (1) are guaranteed if is inferiorly bounded by the inverse of a strictly positive polynomial (vergara2018general).
where is a family of deterministic basis functions and
is a vector of Gaussian weights.lindgren2011explicit provided an expression for the precision matrix of these weights in the special case where is a real polynomial taking strictly positive values on , which corresponds to the case where is a continuous Markov random field.
2.2 Generalized random fields
Let be an isotropic stationary real Gaussian random field on with spectral density . In particular, is a positive radial function of . lang2011fast showed that then, can be seen as a Generalized random field defined by:
where is once again a spatial Gaussian white noise. They used this characterization of Gaussian fields with spectral density to derive algorithms for the fast generation of samples of such fields over rectangular lattices of
using Fast Fourier transform.
A second motivation for this work is to combine this characterization of Gaussian fields with a given spectral density and the SPDE approach to come up with a way to generate samples of these fields on domains more complex than lattices, namely irregular grids, general bounded domains of and even Riemannian manifolds. This type of generalization was in particular exploited by lindgren2011explicit in the particular case of continuous Markov random fields.
In the next section, the approximation of such Generalized random fields using the Finite element method is presented, and an error analysis on this approximation is carried out.
3 Finite element approximation of Generalized random fields
Le and let be either a bounded (convex and polygonal) domain of , or a compact -dimensional smooth Riemannian manifold. Denote , the separable Hilbert space of (real) square-integrable functions on . Denote the inner product of .
Let denote a densely defined, self-adjoint, positive semi-definite linear differential operator of second order, defined in a domain with Dirichlet boundary conditions on . is diagonalizable on a orthonormal basis of, denoted , are arranged so that the eigenvalues satisfy (courant1966methods):
3.1 Theoretical framework
Differential operator on
For , denote . Then we define the action of the differential operator on by:
Remark that the subspace is itself a Hilbert space with respect to the inner product and corresponding norm defined by:
The following lemma gives a sufficient condition so that .
If satisfies then .
, , so in particular . Therefore, such that , and so which allows to conclude that the series is convergent given that the series is convergent. ∎
In the particular case where , where denotes the Laplacian (or the Laplace-Beltrami operator) on , the operator satisfies the following property:
where denotes the extension of the Fourier transform over .
Generalized random fields of
Denote the spatial Gaussian white noise on defined on a complete probability space . A characterization of based on the Hilbert space is given by the following lemma.
Denote . For an integer , take a set of linearly independent elements of and denote . Let’s show that
is a Gaussian vector. Indeed, the characteristic function of this vector is:
Using the fact that the are independent standard Gaussian variables yields:
where is the positive definite matrix whose entries are , . Therefore, we can conclude that is a Gaussian vector, and in particular, by definition of , (10) is satisfied by .
Hence, can be identified to the spatial Gaussian white noise on by defining for the measure where is the indicator function of the set . ∎
Denote the Hilbert space of -valued random variables satisfying and equipped with the scalar product . The next result introduces a class of Generalized random fields defined through the white noise that can be identified with elements of .
Let such that :
Then, denotes the Generalized random field defined by:
can be identified to an element defined by:
for a sequence of independent standard normally-distributed random variables, through the linear functional of : , .
Clearly, is an element of given that
Let’s now show that the linear functional is equal to . Indeed, ,
So, using Lemma 3.3,
From now on, Generalized random fields of the form will be directly identified with their representation, and we will write:
where is a sequence of independent, standard normally-distributed random variables. In the next section, the simulation of such random fields through a finite element scheme is presented.
3.2 Finite element approximation of a Generealized random field
Let denote a triangulation of with mesh size and a family of linearly independent basis functions associated with such that .
Denote the linear span of , which is a -dimensional space. Denote an orthonormal basis of with respect to the scalar product . The discretization of the operator on is denoted and is defined by:
Let and be the matrices defined by:
is a symmetric positive definite matrix called Mass matrix and is a symmetric positive semi-definite matrix called stiffness matrix. Denote the symmetric positive definite square root of , and its inverse.
is diagonalizable on and its eigenvalues are those of the matrix defined by:
In particular, the application:
is an isometric isomorphism that maps the eigenvectors of
is an isometric isomorphism that maps the eigenvectors ofto the eigenfunctions of .
Take an eigenvalue of and denote an associated eigenvector. Then, and so, where . Hence, , , which, using the fact that is self-adjoint, gives
In particular, given that is also a basis of , we denote the invertible change-of-basis matrix between and . Then (8) can be written,
And therefore, , . And so, given that ,
Therefore is an eigenvalue of and maps the eigenvectors of to the eigenfunctions of .
is an isometry: Indeed, , . Consequently, given that it is also linear, is injective. And finally, using the rank–nullity theorem, is bijective (as an injective application between two vector spaces with same dimension). ∎
Let’s denote the eigenvalues of (or equivalently ) and the associated orthonormal eigenfunctions, obtained by mapping by the orthonormal eigenbasis of . In particular,
Take . The discretization of the operator on , denoted , can now be defined as:
Let be a -valued random variable defined by :
where is a vector whose entries are independent zero-mean (standard) normally-distributed random variables. Then, is called white noise on .
Let be white noise in . Then can be written where .
Using the linearity of , can be written where . But also, in the basis , we get . In particular, using the fact that is bijective, which proves the result. ∎
The -valued random variable defined by:
is called finite element approximation of the generalized random field defined by (6).
The finite element approximation defined by (10) can be decomposed as:
for some multi-normally distributed weights with mean and covariance matrix:
Notice that . But also,
Therefore, given that is bijective, where , which proves the result. ∎
Theorem 3.1 provides an explicit expression for the the covariance matrix of the weights of the finite element approximation of a field defined by (6). This expression agrees in particular with the expression of the precision matrix of those same weights exposed in (lindgren2011explicit) for the particular case of continuous Markovian fields on . Given the generality of domains (convex bounded or Riemannian manifold), differential operators and functions , a wide class of fields are now open to study.
Note in particular that Riemannian manifolds provide a generic framework for the study of non-stationary fields on a manifold. Indeed, stationary fields on a Riemannian manifold equipped with a metric can be seen as non-stationary random fields on the manifold with locally-varying anisotropies defined by .
3.3 Error analysis of the finite element approximation
First, we recall the notations used in this paper.
Let be a family of finite element spaces indexed by a mesh size over a domain . Let’s denote be the number of basis functions associated with the triangulation of with mesh size .
Let denote a densely defined, self-adjoint, positive semi-definite linear differential operator of second order defined on a subset of , and let denote its discretization over . Let and be the eigenvalues of and , listed in non-decreasing order.
Assumption 3.1 (Growth of the eigenvalues of ).
There exists tree constants , and such that:
Assumption 3.2 (Derivable of ).
is derivable on , and there exist and such that:
Assumption 3.3 (Asymptotic behavior of ).
There exists a constant such that satisfies , i.e.
Assumption 3.4 (Dimension of the finite element space).
There exists two constants , such that:
Assumption 3.5 (Mesh size).
Assumption 3.6 (Eigenvalues and eigenvectors of ).
There exists constants , and exponents and such that the eigenvalues and the eigenvectors of the discretisation of the operator onto the finite element space associated with a triangulation of mesh size satisfy:
for all and for all
Following the notations of the previous sections, let and be the random fields defined by:
The expected error defined by :
is bounded using the following result.
If , , and satisfy Assumptions (1-6), and if the growth of eigenvalues defined in Assumption 3.1 satisfies:
where the constants are defined in Assumptions (1-6), then, for sufficiently small (cf. Assumption 3.5), the error between the generalized random field defined by (6) and its finite element approximation defined by (10) is bounded by:
where is a constant independent of .
The proof of this theorem is provided in Appendix B and is an adaptation of proof of Theorem 2.10 of (bolin2017numerical), which provide a bound for the same approximation error for the particular case where , , and is a positive definite.
In the particular case where , is a strongly elliptic differential operator of order . Weyl’s law (weyl1911asymptotische) gives an exponent for which Assumption 3.1 is satisfied:
Moreover, if we assume that the finite element spaces are quasi-uniform and are composed of continuous piecewise polynomial functions of degree , then Assuptions 3.4 and 3.6 are satisfied for the exponents (bolin2017numerical; strang1973analysis):
4 Finite element simulations
Simulation of the finite element approximation (10) of fields satisfying (6) comes down to the simulation of the weights of the decomposition onto the basis functions of the finite element space. These weights are Gaussian, with covariance matrix (11).
A straightforward method to generate admissible samples would consist in computing:
where is a vector with independent standard Gaussian components. But this method supposes that the matrix has been diagonalized and that its eigenvalues and its eigenvectors have been stored. Both operations being tremendously costly ( for the diagonalization and for the storage), another method that doesn’t involve them is highly preferable.
We propose to rather simulate the weights using the Chebyshev simulation algorithm presented in (pereira2018efficient). Indeed, the algorithm can provides vector samples with covariance matrix (11) with a computational complexity that is linear with the number of non-zeros (which is small as is a sparse matrix) and an order of approximation which is fixed using a criterion that ensures that the statistical properties of the output are valid. Concerning the storage needs, only a matrix as sparse as needs to be stored as the algorithm relies on matrix-vector multications. This algorithm is reminded below.
In this work we provided an explicit expression for weights of the finite element approximation of a Generalized random field defined over of domain consisting of a bounded domain of or a -dimensional smooth Riemannian manifold, by where , is a second-order self-adjoint positive definite differential operator, and is a Gaussian white noise. An error bound on this approximation was derived and an algorithm for fast and efficient sampling of this field was exposed.
Appendix A Laplacian and Fourier transform
a.1 Multivariate Fourier series and transform
Let be a -periodic function of i.e. is -periodic with respect to each variable and suppose that . Then can be represented as the limit on of its Fourier series defined by (OSBORNE2010115):
where the coefficients are given by:
The Fourier transform of is then defined from its Fourier series as a train of impulses (corresponding to the Fourier transform of each term of Fourier series):
where and is the Dirac impulse at .
Let’s first consider the case and denote . The eigenvalues and eigenfunctions of on with Dirichlet boundary conditions are given for by (grebenkov2013geometrical):
Let and define by:
Then the coefficients of the Fourier series of satisfy :
where , and is an eigenfunction of on with Dirichlet boundary conditions, as defined in (17). (Note : if , otherwise)
In particular, the Fourier series of (restricted to ) coincides (up to a normalization constant) with the development of in the eigenbasis of the Laplacian.
By induction, the same process yields,