Asymmetric MatrixValued Covariances for Multivariate Random Fields on Spheres
Abstract
Matrixvalued 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 matrixvalued 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 cokriging predictor. We refer the reader to [18] for a detailed study on MGRFs and [8] for a review on parametric families of matrixvalued 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 matrixvalued 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 matrixvalued 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 matrixvalued 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 matrixvalued mapping , with th element defined through
The diagonal elements, , are called marginal covariances, whereas the offdiagonal elements, , for , are called crosscovariances. The mapping is (semi) positive definite, i.e., for all positive integer , and , the following inequality holds
Note that, in general, the crosscovariance function is not symmetric in the sense that and . However, when this covariance depends on and only through its geodesic distance (socalled 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 matrixvalued 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
(1) For a more detailed study of the parametric conditions, we refer the reader to [10], [3] and [17].
 Cauchy model

The conditions for matrixvalued 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
(2) Moreover, [17] provide explicit conditions on the parameters for the variate case, with .
 Compactly supported models

[15] shows that a radial matrixvalued 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])
(3) 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]).
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
(4) 
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 crosscovariances , 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
(5) 
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
where
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 ,
where
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 cutoff 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 nonseparable. We make simulation experiments for both cases. For each nonseparable 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 .
We assess the typical variance of the likelihoodbased 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 cutoff 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 nonseparable 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.
Separable  Nonseparable  
S  A  S  A  S  A  S  A  
(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 
Separable  Nonseparable  
S  A  S  A  S  A  S  A  
(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 
Now, we study the predictive performance of asymmetric models through a crossvalidation analysis. We simulate bivariate fields from the asymmetric versions of models (M1)–(M3), with both separable and nonseparable 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 cokriging predictor, with both the symmetric and asymmetric covariances, in terms of the mean squared prediction error (MSPE) and the logscore (LSCORE) (see [20]). These indicators are evaluated using a dropone 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.
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. Nonseparable exponential model without asymmetry.

Model 3. Separable exponential model after adding asymmetry.

Model 4. Nonseparable 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 LogComposite Likelihood (LogCL) 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 dropone 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 nonseparability.
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 
parameters  LogCL  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 
6 Conclusions
We have provided a strategy to construct asymmetric matrixvalued 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.
Acknowledgments
Alfredo Alegría is supported by Beca CONICYTPCHA/Doctorado Nacional/201621160371. Emilio Porcu is supported by Proyecto Fondecyt Regular number 1130647. Reinhard Furrer acknowledges support of the Swiss National Science Foundation SNSF144973. We additionally acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for Coupled Model Intercomparison Project (CMIP).
References
 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 crosscovariance 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). CrossCovariance Functions for Multivariate Geostatistics. Statistical Science, 30(2):147–163.
 Gneiting [2013] Gneiting, T. (2013). Strictly and nonstrictly positive definite functions on spheres. Bernoulli, 19(4):1327–1349.
 Gneiting et al. [2010] Gneiting, T., Kleiber, W., and Schlather, M. (2010). Matérn crosscovariance 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érnlike 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). SpatioTemporal Covariance and CrossCovariance 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 crossvalidation for massive spatial data. Environmetrics, 21(34):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.