Scalable modeling of nonstationary covariance functions with non-folding B-spline deformation

07/01/2020 ∙ by Ronaldo Dias, et al. ∙ University of Campinas University of Washington 0

We propose a method for nonstationary covariance function modeling, based on the spatial deformation method of Sampson and Guttorp [1992], but using a low-rank, scalable deformation function written as a linear combination of the tensor product of B-spline basis. This approach addresses two important weaknesses in current computational aspects. First, it allows one to constrain estimated 2D deformations to be non-folding (bijective) in 2D. This requirement of the model has, up to now,been addressed only by arbitrary levels of spatial smoothing. Second, basis functions with compact support enable the application to large datasets of spatial monitoring sites of environmental data. An application to rainfall data in southeastern Brazil illustrates the method



There are no comments yet.


page 2

page 8

page 9

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Geostatistical methods for spatial and spatio-temporal data are in great demand from fields such as earth and climate sciences, epidemiology and agriculture. A comprehensive overview can be found in cressie2011statistics. Spatially stationary processes are commonly used as models in geostatistical applications, but often the assumption of stationarity and isotropic covariance functions are difficult to hold in real applications; see for example, guttorp1994space, le2006statistical, damian2001bayesian. A recent review of methods that allow nonhomogenous covariance models is schmidt2020flexible.

We propose a semiparametric method of nonstationary spatial covariance function that expands upon the work of sampson1992nonparametric. Data observed in a set

2 Spatial Deformation Model (SDM)

Let be spatial locations with for all in a geostatistical domain and with for all Moreover, let such that for all and are bijective, differentiable functions, for

Let be a Gaussian random field with spatial covariance function For any pair of spatial sites we have with given by

where and is a stationary, isotropic correlation function with parameters We remark that this model may not be identifiable, since is invariant for shifts and rigid rotations in Moreover, scaling by a constant will produce the same if is also scaled by

Nevertheless the model can be used for Kriging and spatial interpolation without issues.

For example, consider and given by the swirl transformation shown in Figure 1. We sampled from a spatio-temporal process on a regular grid of points, with and separable covariance function

where and A realization of a random field on a fine mesh grid is also shown in Figure 1, obtained with the RandomFields package (schlather2015analysis).

Figure 1: Left panel: regular grid in in black, and deformed grid in in red. The contour plots in the margins are the the functions and that provide the deformation. Right panel: simulated random field on using a deformed exponential covariance function that is stationary on

3 Spatial Deformation Estimation


then is the semivariogram of the random field . In particular, for a stationary, isotropic

where is a positive non-decreasing function of . Note valid functions are conditionally negative definite (cressie1993statistics).

If there are ways to obtain a sample covariance matrix, such as temporal replicates, then we can construct a sample variogram as

If the variogram is isotropic on a set of artificial coordinates , then it must be have form

where Therefore, the variogram entries can be seen as a dispersion metric, with

3.1 Non-metric multidimensional scaling

sampson1992nonparametric first step is to consider a non-metric multidimensional scaling approach (hereafter, nMDS; see kruskal1964multidimensional; mardia1979multivariate; cox2000multidimensional). From a set of dispersions we seek a monotone transformation such that


and an artificial set of coordinates with interpoint Euclidean distances that minimizes the stress criterion

Here we need to impose constraints on such that the function is conditionally positive definite. kruskal1964multidimensional only finds an isotonic transformation, so sampson1992nonparametric adapted the optimization algorithm as follows:

  1. Set k = 1. Find an initial configuration using, for example, Kruskal’s isotonic nMDS.

  2. As in Section 3.2 (or Section 4 using our proposed smoother) we smooth the artificial coordinates using the corresponding spline method. Denote them by

  3. Fit a variogram model to the data using the coordinates obtaining

  4. Obtain via MDS performed on

  5. Set and return to step 2 until convergence.

3.2 Spline smoothing

The second step in sampson1992nonparametric approach is to find an approximation for such that This can be done if we approximate the functions with, for example, thin-plate splines (wahba1980some). Thin-plate splines are the minimizers of the variational problem

The solution for a fixed is

where and are estimating by plug-in.

3.3 Tensor product of B-splines

We propose a different approach for the spline smoothing by making use of low rank approximation (finite approximation) for the deformation functions. Specifically, let and assume we have a collection of points be a sequence of points in . Define a map such that


where the functions , form a partition of unity.

Particularly, suppose that are the B-splines basis functions and there is a integer number such that a map can be represented by


where and Let for some index set such that then the functions are well approximated by

where are fixed positive integers and are B-spline basis functions (see, e.g., ramsay2005functional).

To guarantee that the functions do not fold onto themselves, we must guarantee that is locally invertible and differentiable. In fact, a diffeomorphism is desirable (see, e.g., perrin1999modelling). A necessary condition is that everywhere in the domain,


but in practice such constraint is difficult to implement, requiring the evaluation of the Jacobian for every pair We will show in Section 4 that such condition can be translated into a condtion onto the coefficients for Note however that it is enough to ensure since a change in signs would imply a discontinuous Jacobian.

3.4 Simultaneous estimation of covariance function and deformation

The covariance parameters

as well as the mean vector

are estimated by maximizing the profile log-likelihood. Without loss of generality set Thus our proposed B-spline approach can be estimated as

  1. Given solve the optimization problem

    where corresponds to the optimization target seen as a function of only, with parameters fixed.

  2. Given obtain such that


    where are the row vector of basis functions and is the Kronecker product, and for some choice of such that if

  3. Repeat steps 1 and 2 until convergence.

We remark that the matrices are sparse when using B-splines, therefore solutions to step 2 are scalable. To evaluate the inequality consider the following: write

or simply where is a vector of B-spline basis functions evaluated at , similarly is a vector of B-spline basis evaluated at , and is a matrix of spline coefficients for the tensor product approximation of the -th function . Similarly,

is the partial derivative of with respect to and if we write then This allows us to see that

so the determinant of the Jacobian as a function of is an inner product of

, weighted by the skew-Symmetric matrix

Ensuring the inequality for all values of remains a difficult task, but we can chose a set of basis for which does not depend on

4 Constrained spatial deformation estimation

We will use an approach similar to musse2001topology. Consider B-splines of degree 1 on . Assume that there are equally spaced inner knots, where . Since the knots are equally spaced, they can be written as where In this case, the B-spline bases are given by

with derivatives

Consider where and is between and (where and ). Then there are only 4 bases that evaluate to non-zero values, indexed by , , and , so

and therefore, looking only at the non-zero pairs , we have


note is therefore proportional to

The above equations describe a plane in with coefficients depending on . Now, since and , where (and similarly for ), we have four restrictions to consider:

  • When and ,

  • When and ,

  • When and ,

  • When and ,

This collection of constraints, for and together imply in a non-folding deformation map, and can be enforced with constrained optimization routines. We have employed svanberg2002class constrained optimization algorithm, avaliable in the nloptr package (johnson2020nloptr). The code is available as an R package in

5 Simulation study

To evaluate the performance of the algorithm, we conducted a simulation study based on the swirl function shown in Figure 1. Each sample has spatial points on a regular grid in the geographical domain We simulated from a Gaussian random field with mean function and covariance function with parameters and First, consider a single realization of the spatial random field. The estimated deformation maps are shown in Figure 2. We have obtained the constrained B-spline deformations (hereafter, bdef) with and basis functions. The estimated deformation function does not fold, even though the true deformation function is difficult to be recovered. The case when performs better than or

indicating that there are features in the deformation map that need a large number of degrees of freedom to be estimated. On the other hand, the functions using

sampson1992nonparametric (SG) method have smoothing parameters and set to match the bdef approach. They start showing folding at and cannot recover deformation maps that require a number of degrees of freedom larger than Note that the SG maps were stretched or shrunk to fit the plot area, but the bdef maps did not require this step.

Figure 2: Comparison of estimated deformation functions for simulated data. The upper row corresponds to the proposed regularized B-spline approach with and respectively. The bottom row corresponds to the sampson1992nonparametric method with and which are equivalent in degrees of freedom to the upper cases.

The comparison of estimated covariance matrices allow us to overlook the identifiability issues with rotations, shifts and scaling of the estimated maps. In Figure 3 we show a scatterplot of the upper-diagonal entries of the estimated covariance matrices for the data, versus the true covariance matrices. We remark that the bdef method shows no apparent bias and becomes more accurate as the number of degrees of freedom increase. On the other hand, SG has good performance at a small number of degrees of freedom (more smoothing), but tends to overestimate matrix entries as the number of degrees of freedom for the thin-plate spline increases.

Figure 3: Comparison of upper-diagonal entries for the estimated covariance matrix of the simulated data, versus the true covariance matrix. The upper row corresponds to the proposed regularized B-spline approach with and respectively. The bottom row corresponds to the sampson1992nonparametric method with and which are equivalent in degrees of freedom to the upper cases.

6 Case study: Rainfall data in southeastern Brazil

The dataset we use to illustrate our method comes from meteorological surveys conducted by INMET – Instituto Nacional de Meteorologia, Brazil. The measurements are made at every 15 minutes, and daily accumulated values are made available. Following rozante2010combining, we grouped the rainfall data in periods of 10 days.

Since we seek temporally stationary data, we decided to restrain the data collection to 2018-01-01 to 2018-03-30 (rainfall season), for a total of 9 time periods (period 1: 2018-01-01 to 2018-01-10, …, period 9: 2018-03-22 to 2018-03-31).

We selected the 50 stations of the southeastern region of Brazil that had complete observations available during the period. The stations are shown in the left panel of Figure 4. Data can be obtained at

In the right panel of Figure 4, we show the estimated deformation map, using B-spline basis functions. The estimated deformation map reveals topographical features shown in the left panel of the same Figure. Weather stations located in northwestern flat lands of Minas Gerais are treated as closer to each other than stations near the rough coast of Rio de Janeiro and São Paulo states, where elevation changes are abrupt.

Figure 4: Left panel: map of meteorological stations located in southeastern Brazil. Right panel: estimated deformation functions using B-spline basis.

In Figure 5 we perform conditional simulation of the Gaussian random field (Kriging) for the 10-day periods starting in 2018–01–01 and 2018–01–11, showing the resulting prediction maps.

Figure 5: Conditional simulation of the Gaussian random field (Kriging) for the 10-day periods starting in 2018–01–01 and 2018–01–11.


Ronaldo Dias was supported by the FAPESP foundation, grant 2018/04654-9. Guilherme Ludwig was supported by FAPESP grants 2018/04654-9 and 2019/03517-0, as well as FAEPEX grant 2253/19.


Appendix A Derivation of the non-folding constraint

We can derive the constraint presented in Section 4 as follows: Observe that


Further simplifying leads to