Asymmetric Matrix-Valued Covariances for Multivariate Random Fields on Spheres

Asymmetric Matrix-Valued Covariances for Multivariate Random Fields on Spheres

\nameAlfredo Alegría1, Emilio Porcu2,1 and Reinhard Furrer3 Alfredo Alegría. Email: alfredo.alegria.jimenez@gmail.comEmilio Porcu. Email: georgepolya01@gmail.comReinhard Furrer. Email: 1Department of Mathematics, University Federico Santa Maria, Valparaiso, Chile;
2School of Mathematics and Statistics, University of Newcastle, Newcastle upon Tyne, United Kingdom;
3Department of Mathematics and Department of Computer Science, University of Zurich, Zurich, Switzerland.

Matrix-valued covariance functions are crucial to geostatistical modeling of multivariate spatial data. The classical assumption of symmetry of a multivariate covariance function is overlay restrictive and has been considered as unrealistic for most of real data applications. Despite of that, the literature on asymmetric covariance functions has been very sparse. In particular, there is some work related to asymmetric covariances on Euclidean spaces, depending on the Euclidean distance. However, for data collected over large portions of planet Earth, the most natural spatial domain is a sphere, with the corresponding geodesic distance being the natural metric. In this work, we propose a strategy based on spatial rotations to generate asymmetric covariances for multivariate random fields on the -dimensional unit sphere. We illustrate through simulations as well as real data analysis that our proposal allows to achieve improvements in the predictive performance in comparison to the symmetric counterpart.


auchy model; Geodesic distance; Global data; Rotation group; Wendland model.

1 Introduction

In recent years, the use of multivariate geostatistical models has become increasingly popular in disciplines such as climatology, environmental sciences and oceanography. Specifically, the georeferenced observations are assumed as a partial realization of a multivariate Gaussian random field (MGRF) and the spatial dependencies between the observed data are modeled through the matrix-valued covariance function of the field. Covariance functions are the main tool for spatial interpolation methods and the development of accurate predictions. The best linear unbiased predictor, under a mean squared error criterion, is known in multivariate geostatistics as the co-kriging predictor. We refer the reader to [18] for a detailed study on MGRFs and [8] for a review on parametric families of matrix-valued covariance functions.

A typical assumption in geostatistics is that of spatial isotropy, which means that the multivariate covariance depends solely on spatial distance. The practical and theoretical contributions in spatial statistics have been mainly developed on Euclidean spaces with the corresponding Euclidean metric. This choice is supported by the existence of several methods to project the geographical locations onto the plane. However, for data recorded over large portions of planet Earth, such projection generate significant distortions in the distances. The problem of an appropriate metric has been discussed by [4], [9] and [17]. Indeed, if we agree that a sphere is a good approximation of planet Earth, it is more realistic to work under the framework of random fields indexed on a sphere [16]. Moreover, the most natural metric to describe the spatial interactions is the geodesic distance, which corresponds to the shortest path joining two points on the spherical surface.

The construction of flexible families of matrix-valued covariance functions on spheres is a challenging problem, since the corresponding covariance matrices must be positive definite. The literature in this direction has been very sparse with the work of [17], who derive conditions for the validity of the multivariate Matérn and Cauchy models, being a notable exception. On the other hand, [15] shows that compactly supported models, introduced on Euclidean spaces, can be adopted on specific low dimensional spheres. Additional approaches based on latent dimensions are explored by [2].

The parametric covariance families mentioned above depend solely on the geodesic distance (geodesic isotropy). In particular, spatial isotropy implies symmetry of the matrix-valued covariance structure. However, in real data applications, asymmetrical behaviors are usually observed. Indeed, we need to go beyond the class of geodesically isotropic covariances to generate asymmetric models. In this paper, we propose a strategy based on spatial rotations to generate asymmetric covariances from symmetric models. We show that the inclusion of asymmetry produce more flexibility in the statistical modeling, and particularly, that the improvements in terms of predictive performance are significant. Our proposal is the spherical counterpart of the work developed by [14] on Euclidean spaces.

The remainder of the article is as follows. Section 2 contains the background material on MGRFs on spheres as well as a brief summary of some symmetric covariance families, which are the building block to construct asymmetric models. Section 3 introduces a general approach to produce asymmetric covariances. In Section 4, we show through simulation experiments that asymmetric models can provide better predictive results in comparison to the symmetric counterparts. A bivariate data set is analyzed in Section 5. Finally, Section 6 contains some conclusions.

2 Background

We briefly review the literature on MGRFs on spheres and the parametric families of matrix-valued covariance functions introduced in the literature. Although the models considered in this section are symmetric, they are the building block to construct explicit asymmetric models in the subsequent sections.

2.1 MGRFs on Spheres

Throughout the article, we consider MGRFs, , defined on the -dimensional unit sphere , , where denotes the number of components of the field, the transpose operator and the Euclidean norm in . We call a -variate random field. The natural metric on the spherical surface is the geodesic distance defined through the mapping given by .

We assume, without loss of generality, that has zero mean, since we are essentially interested on its covariance structure. Let and consider the matrix-valued mapping , with -th element defined through

The diagonal elements, , are called marginal covariances, whereas the off-diagonal elements, , for , are called cross-covariances. The mapping is (semi) positive definite, i.e., for all positive integer , and , the following inequality holds

Note that, in general, the cross-covariance function is not symmetric in the sense that and . However, when this covariance depends on and only through its geodesic distance (so-called geodesic isotropy), such a covariance function is clearly symmetric. The rest of this section focuses on geodesically isotropic models.

Following [17], we call the class of continuous mappings , being continuous and such that, for a positive definite function , we have . The mapping is called radial part of a geodesically isotropic covariance . Moreover, the inclusion , for all positive integer , is strict.

Throughout, denotes the -ultraspherical polynomial of degree [1], which is defined implicitly as

The following result [11, 19] characterizes completely the class in terms of ultraspherical polynomials. The equalities and summability conditions must be understood in a componentwise sense.

Theorem 2.1.

Let and be positive integers. Let be a continuous mapping with elements , such that , for all . Then, is a member of the class if and only if it admits the representation

where is a sequence of symmetric, positive definite and summable matrices of order .

2.2 Parametric Families in

We now provide a review on parametric families of geodesically isotropic matrix-valued covariance functions. We focus on the Matérn and Cauchy families, as well as on compactly supported models. In particular, the Matérn and Cauchy families are valid on spheres of any dimension, whereas the compactly supported models require additional restrictions on the dimension.

Matérn model

Let be the modified Bessel function of second kind and be the Gamma function [1]. The multivariate Matérn model on the sphere [17] is given by

provided that and are positive, , , with , for all . Here, is the standard deviation of the -th component of the field, are scale parameters, is the collocated correlation coefficient between the components and , and is a smoothness parameter. Some additional conditions for the parameters are required (see [17]). In the following sections, we focus on parsimonious bivariate () models and such conditions are relaxed. For instance, the parameterization has associated condition


For a more detailed study of the parametric conditions, we refer the reader to [10], [3] and [17].

Cauchy model

The conditions for matrix-valued covariance functions of generalized Cauchy type on spheres have been developed by [17]. The generalized Cauchy model is given by

with , and being positive, , and , for all . Here, the interpretation of the parameters is similar to the Matérn model. Additional conditions must be imposed in order to obtain a positive definite mapping (see [17]). For instance, in a bivariate case with , we have the restriction


Moreover, [17] provide explicit conditions on the parameters for the -variate case, with .

Compactly supported models

[15] shows that a radial matrix-valued positive definite function on , with elements having compact support less than , can be coupled with the geodesic distance to obtain a geodesically isotropic positive definite function on .

Formally, consider the mapping , such that , with elements , is positive definite on , with , for all and . Therefore, the mapping defined by belongs to , where denotes the restriction of to a Borel set .

In particular, the inclusion between the classes imply that, under this construction principle, we can obtain models on spheres of dimension 1, 2 and 3. Examples can be generated from Table 1 in [6]. For instance, the Wendland model is defined as

where , , , and , for all . Here and the interpretation of the parameters is similar to the previous examples. In addition, we have the condition (see [6])


In the subsequent sections, we consider a bivariate Wendland model with the parsimonious choice .

The characterization of the cross scale parameter as a function of the marginal ones, and , is a very useful strategy to avoid an excessive number of parameters in the estimation procedure. In some cases, this parsimonious strategy can provide better predictive results than the full model (see [10]).

Additional parametric models can be constructed using latent dimensions [2], or by extending the spherical convolution approach [21] as well as covariance approximation methods [12] to the multivariate case.

3 Asymmetric Covariances on

A strategy to construct asymmetric models is now provided. We start with a -variate zero mean Gaussian field, , with geodesically isotropic, and thus symmetric, covariance structure .

Let be a collection of rotation matrices, i.e., is an orthogonal matrix with determinant one, for each . Such matrices satisfy the relation . Consider a new -variate random field defined through


The field introduced in Equation (4) is Gaussian and it has zero mean. Furthermore, the covariance function associated to , denoted as , is given by

for each . Obviously, preserves the marginal structure associated to , since , for each . However, the cross-covariances , for , are not necessarily geodesically isotropic. Thus, does not belong to the class . Moreover, we have that , for .

Here, the rotation matrices are handled as additional parameters of the model. In order to avoid identifiability problems in the estimation of such parameters, we need to add some constraints in a similar fashion to the work of [14]. Since only the product is relevant, we propose the condition , where denotes the identity matrix of order . We illustrate these considerations for the cases and .

Example 3.1.

In the unit circle , we have the general rotation matrix

where denotes the rotation angle. For each , we define . Thus, the constraint is equivalent to the condition . In particular, we use the following identifiability condition


Direct inspection shows that and , where and .

For instance, consider a bivariate scenario with and , for . In this case, we have the additional asymmetry parameter . The general expression for the bivariate asymmetric covariance functions on , generated from our approach, is given by


with denoting the Hadamard product. Recall that is a summable sequence of symmetric positive definite matrices. In this case, the ultraspherical polynomials reduces to the Tchebyshev polynomials .

Example 3.2.

Consider now the unit sphere . The Rodrigues’s rotation formula provides the following representation for a rotation matrix, with respect to an axis by an angle ,


Here, we take the collection of matrices . Note that we have the relation . For the bivariate case (), considering the identifiability condition (5), and , we have the additional parameters . Since the pairs and define the same rotation, we must to restrict the range of to the interval . The axis of rotation can be parameterized in terms of spherical coordinates: , with and . In conclusion, we have the vector of additional parameters .

4 Simulation Study

We illustrate simulation experiments related to bivariate Gaussian fields on with asymmetric covariances. We use initial covariance functions either of Matérn, Cauchy and Wendland type, as introduced in the previous sections. In particular, we consider the mentioned models with the following parameterizations:

  • Matérn model with , which reduces to the exponential model

    where .

  • Cauchy model with ,

    where .

  • Wendland model with ,

    where .

The parameterization of (M1) and (M2) ensures that for . For model (M3), is exactly equal to 0 beyond the cut-off distance . The large number of parameters considered in our simulation experiments supports the parsimonious parameterization used for . In particular, for the generalized Cauchy model (M2) we consider as the average of the marginal scale parameters, since under this assumption we have the condition (2), which can be checked easily. The same happens for the other models (see conditions (1) and (3)). Figure 1 shows these covariances under a specific parameter setting: , , , .

In the special case where , for all , we say that the covariances are separable, otherwise they are non-separable. We make simulation experiments for both cases. For each non-separable model, after adding asymmetry, the vector of parameters is given by . In the separable case, on the other hand, the vector of parameters reduces to . In our studies, we consider 225 spatial sites in a grid generated combining 15 equispaced azimuthal and polar angles. Specifically, the polar angles are taken as , whereas the azimuthal angles as , for .

Figure 1: Parametric families of matrix-valued covariances (M1)–(M3) considered in our simulation experiments. Here , , , .

We assess the typical variance of the likelihood-based estimates when the asymmetric effect is added. In particular, we consider the pairwise composite likelihood (CL) approach developed by [5], which is a fast likelihood approximation method. We set , , , , and . Under the separable scenario, we set . Figure 2 reports the bias of the CL estimates for each model on the basis of 500 independent simulated data sets. We have used a special case of CL method, by considering only pairs of data whose spatial distance is less than 1 radian, although in our experience, smaller cut-off distances can also ensure stability of the method. Even though we are dealing with a large number of parameters, the boxplots show that the direct maximization of the CL objective function provides quite reasonable results. The empirical confidence interval for is in the non-separable model (M1), and similar intervals are obtained for the other models. Another appealing estimation procedure is to consider a two steps maximization routine (see [14]), which is useful in terms of computational complexity.

Figure 2: Bias of the pairwise CL estimates, on the basis of 500 independent repetitions. We consider models (M1)–(M3) (from top to bottom), with separable (left) and non-separable (right) structures.
Separable Non-separable
(M1) MSPE 0.983 0.968 0.964 0.935 0.964 0.953 0.954 0.924
LSCORE 1.410 1.402 1.398 1.383 1.396 1.386 1.387 1.373
(M2) MSPE 0.969 0.956 0.961 0.930 0.924 0.910 0.929 0.893
LSCORE 1.398 1.393 1.395 1.381 1.357 1.348 1.352 1.327
(M3) MSPE 0.981 0.966 0.980 0.959 0.993 0.988 0.972 0.945
LSCORE 1.409 1.399 1.407 1.396 1.414 1.409 1.405 1.388
Table 1: Cross-validation scores for models (M1)–(M3), separable and non-separable, using symmetric (S) and asymmetric (A) versions, with .
Separable Non-separable
(M1) MSPE 0.969 0.945 0.979 0.919 0.940 0.919 0.950 0.882
LSCORE 1.401 1.389 1.405 1.369 1.383 1.372 1.387 1.345
(M2) MSPE 0.971 0.950 0.975 0.908 0.938 0.904 0.924 0.856
LSCORE 1.400 1.390 1.406 1.366 1.356 1.339 1.358 1.326
(M3) MSPE 0.992 0.969 0.993 0.911 0.981 0.951 0.970 0.897
LSCORE 1.415 1.405 1.414 1.364 1.408 1.393 1.403 1.366
Table 2: Cross-validation scores for models (M1)–(M3), separable and non-separable, using symmetric (S) and asymmetric (A) versions, with .

Now, we study the predictive performance of asymmetric models through a cross-validation analysis. We simulate bivariate fields from the asymmetric versions of models (M1)–(M3), with both separable and non-separable parameterizations. We set , , , , as in the previous example, and consider two scenarios for the collocated correlation coefficient: , and two scenarios for the asymmetry parameter: . Again, in the separable case, we consider a similar experiment with . We quantify the accuracy of the co-kriging predictor, with both the symmetric and asymmetric covariances, in terms of the mean squared prediction error (MSPE) and the log-score (LSCORE) (see [20]). These indicators are evaluated using a drop-one prediction strategy. Tables 1 and 2 report the results for and , respectively. It is clearly visible that the asymmetric models produce better results. Note that if the asymmetry and the correlation between the fields are strong, then the improvements are more significant.

5 A Real Data Example

This section assesses the statistical performance of asymmetric models on a real data example. We analyze a bivariate data set of Temperatures (T) and Precipitations (P). These data outputs have been obtained from the Community Climate System Model (CCSM4.0, [7]) provided by National Center for Atmospheric Research (NCAR), CO, USA. The joint modeling of these variables is of major interest in climatological disciplines.

Figure 3: Residuals for temperature (left) and cubic root of precipitation (right) of model output of December 2015.

The spatial resolution is 2.5 degrees in latitude and longitude. We assume planet Earth as a sphere of radius 6378 kilometers. We pay attention to the geographical region delimited by the longitudes and latitudes degrees (see Figure 3) and consider December 2015. The maximum geodesic distance between the points in this region is 70% of the maximum distance on planet Earth. We work with the cubic root of precipitations in order to attenuate the skewness of this variable. The units are Kelvin degrees for temperatures and centimeters for transformed precipitations. In order to achieve geodesic isotropy, we remove through splines the trend and cyclic patterns of the variables, using longitude and latitude as covariates. The resulting residuals are approximately Gaussian with zero mean.

We consider the bivariate exponential model (M1) with the following features:

  • Model 1. Separable exponential model without asymmetry.

  • Model 2. Non-separable exponential model without asymmetry.

  • Model 3. Separable exponential model after adding asymmetry.

  • Model 4. Non-separable exponential model after adding asymmetry.

Table 3 illustrates the pairwise CL estimates for each model. We have considered CL method with pairs of data whose geodesic distance is less than 1000 kilometers. The estimates of the scale parameters are given in kilometers. We compare the models in terms of the Log-Composite Likelihood (Log-CL) value at the optimum, as well as through their predictive performance. The accuracy of the prediction, for each model, is evaluated in terms of MSPE and LSCORE, using a drop-one prediction strategy (see Table 4). It is clear that the asymmetric versions provide better results with respect to the symmetric ones. Note that for this specific data set, Model 3 outperforms Model 2, which means that the inclusion of asymmetry can produce more flexibility in comparison to the inclusion of non-separability.

Model 1 34.643 0.342 5102 5102 5102 - - -
Model 2 34.275 0.342 4789 5414 5414 - - -
Model 3 34.524 0.342 5089 5089 5089 0.043 1.512 2.678
Model 4 34.098 0.342 4815 5427 5427 0.047 1.504 2.236
Table 3: Composite Likelihood estimates for Models 1-4.
parameters Log-CL MSPE LSCORE
Model 1 4 20460 1.309
Model 2 5 20485 1.306
Model 3 7 20500 1.157
Model 4 8 20518 1.139
Table 4: Comparison of Log-CL values and cross-validation scores for Models 1-4.

6 Conclusions

We have provided a strategy to construct asymmetric matrix-valued covariance functions on spheres. We show that our proposal can produce significant improvements in the predictive performance of the model. Our developments have been illustrated through simulated data, under several models and parametric settings. We have also analyzed a real data set of temperatures and precipitations on a large portion of planet Earth, where the asymmetric models produce better results.

An interesting research direction is the construction of parametric tests, based on the asymmetry parameter , to contrast the hypothesis of symmetry in a given data set. In particular, a challenging problem is the adaptation of the methodology used by [13], where the authors evaluate several types of common assumptions on multivariate covariance functions, to the spherical context.


Alfredo Alegría is supported by Beca CONICYT-PCHA/Doctorado Nacional/2016-21160371. Emilio Porcu is supported by Proyecto Fondecyt Regular number 1130647. Reinhard Furrer acknowledges support of the Swiss National Science Foundation SNSF-144973. We additionally acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for Coupled Model Intercomparison Project (CMIP).


  • Abramowitz and Stegun [1965] Abramowitz, M. and Stegun, I. A. (1965). Handbook of mathematical function: with formulas, graphs and mathematical tables. Dover Publications.
  • Alegría et al. [2017] Alegría, A., Porcu, E., Furrer, R., and Mateu, J. (2017). Covariance Functions for Multivariate Gaussian Fields Evolving Temporally Over Planet Earth. arXiv preprint, arXiv:1701.06010.
  • Apanasovich et al. [2012] Apanasovich, T. V., Genton, M. G., and Sun, Y. (2012). A valid Matérn class of cross-covariance functions for multivariate random fields with any number of components. Journal of the American Statistical Association, 107(497):180–193.
  • Banerjee [2005] Banerjee, S. (2005). On geodetic distance computations in spatial modeling. Biometrics, 61(2):617–625.
  • Bevilacqua et al. [2016] Bevilacqua, M., Alegria, A., Velandia, D., and Porcu, E. (2016). Composite Likelihood Inference for Multivariate Gaussian Random Fields. Journal of Agricultural, Biological, and Environmental Statistics, 21(3):448–469.
  • Daley et al. [2015] Daley, D., Porcu, E., and Bevilacqua, M. (2015). Classes of compactly supported covariance functions for multivariate random fields. Stochastic Environmental Research and Risk Assessment, 29(4):1249–1263.
  • Gent et al. [2011] Gent, P. R., Danabasoglu, G., Donner, L. J., Holland, M. M., Hunke, E. C., Jayne, S. R., Lawrence, D. M., Neale, R. B., Rasch, P. J., Vertenstein, M., Worley, P. H., Yang, Z.-L., and Zhang, M. (2011). The community climate system model version 4. Journal of Climate, 24(19):4973–4991.
  • Genton and Kleiber [2015] Genton, M. G. and Kleiber, W. (2015). Cross-Covariance Functions for Multivariate Geostatistics. Statistical Science, 30(2):147–163.
  • Gneiting [2013] Gneiting, T. (2013). Strictly and non-strictly positive definite functions on spheres. Bernoulli, 19(4):1327–1349.
  • Gneiting et al. [2010] Gneiting, T., Kleiber, W., and Schlather, M. (2010). Matérn cross-covariance functions for multivariate random fields. Journal of the American Statistical Association, 105(491):1167–1177.
  • Hannan [2009] Hannan, E. (2009). Multiple Time Series. Wiley Series in Probability and Statistics. Wiley.
  • Jeong and Jun [2015] Jeong, J. and Jun, M. (2015). A class of matérn-like covariance functions for smooth processes on a sphere. Spatial Statistics, 11:1–18.
  • Li et al. [2008] Li, B., Genton, M. G., and Sherman, M. (2008). Testing the covariance structure of multivariate random fields. Biometrika, pages 813–829.
  • Li and Zhang [2011] Li, B. and Zhang, H. (2011). An approach to modeling asymmetric multivariate spatial covariance structures. Journal of Multivariate Analysis, 102(10):1445–1453.
  • Ma [2016] Ma, C. (2016). Stochastic representations of isotropic vector random fields on spheres. Stochastic Analysis and Applications, 34(3):389–403.
  • Marinucci and Peccati [2011] Marinucci, D. and Peccati, G. (2011). Random fields on the sphere: representation, limit theorems and cosmological applications, volume 389. Cambridge University Press.
  • Porcu et al. [2016] Porcu, E., Bevilacqua, M., and Genton, M. G. (2016). Spatio-Temporal Covariance and Cross-Covariance Functions of the Great Circle Distance on a Sphere. Journal of the American Statistical Association, 111(514):888–898.
  • Wackernagel [2003] Wackernagel, H. (2003). Multivariate Geostatistics: An Introduction with Applications. Springer, New York, 3rd edition.
  • Yaglom [1987] Yaglom, A. M. (1987). Correlation Theory of Stationary and Related Random Functions. Volume I: Basic Results. Springer, New York.
  • Zhang and Wang [2010] Zhang, H. and Wang, Y. (2010). Kriging and cross-validation for massive spatial data. Environmetrics, 21(3-4):290–304.
  • Ziegel [2014] Ziegel, J. (2014). Convolution roots and differentiability of isotropic positive definite functions on spheres. Proceedings of the American Mathematical Society, 142(6):2063–2077.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description