Tucker Tensor Analysis of Matérn Functions in Spatial Statistics

2.1 Problem Settings in Spatial Statistics

Below, we formulate five tasks. These computational tasks are very common and important in statistics. Fast and efficient solution of these tasks will help to solve many real-world problems, such as the weather prediction, moisture modeling, and optimal design in geostatistics.

2.2 Matérn Covariance and Its Fourier Transform

A low-rank approximation of the covariance function is a key component of the tasks formulated above. Among of the many covariance models available, the Matérn family [45, 25] is widely used in spatial statistics, geostatistics [7], machine learning [4], image analysis, weather forecast, moisture modeling, and as the correlation for temperature fields [49]. The work [25] introduced the Matérn form of spatial correlations into statistics as a flexible parametric class with one parameter determining the smoothness of the underlying spatial random field.

The main idea of this low-rank approximation is shown in Figure 1 and explained in details in Section 3.3. Figure 1 demonstrates two possible ways to find a low-rank tensor approximation of the Matérn covariance function. The first way (marked with “?”) is not so trivial and the second via the Fast Fourier Transform (FFT), low-rank and the inverse FFT (IFFT) is more trivial. We use here the fact that the FT of the Matérn covariance is analytically known and has a known low-rank approximation. The IFFT can be computed numerically and does not change the tensor ranks.

Figure 1

Two possible ways to find a low rank tensor approximation of the Matérn covariance matrix Cν,ℓ⁢(r).

The Matérn covariance function is defined as

where distance r:=∥x-y∥, x,y two points in ℝd, ν>0 defines the smoothness of the random field, and 𝒦ν denotes the modified Bessel function of order ν. The larger ν, the smoother the random field. The parameter ℓ>0 is a spatial range parameter that measures how quickly the correlation of the random field decays with distance, with larger ℓ corresponding to a slower decay (keeping ν fixed). When ν=12 (see [55]), the Matérn covariance function reduces to the exponential covariance model and describes a rough field. The value ν=∞ corresponds to a Gaussian covariance model, which describes a very smooth field, that is infinitely differentiable. Random fields with a Matérn covariance function are ⌊ν-1⌋ times mean square differentiable. Thus, the hyperparameter ν controls the degree of smoothness.

The d-dimensional Fourier transform 𝑭d(C(r,ν) of the Matérn kernel, defined in equation (2.2), in ℝd is given by [45]

where β=β⁢(ν,ℓ,n) is a constant and |ξ| is the Euclidean distance in ℝd. The following tensor approach also applies to the case of anisotropic distance, where r2=〈A⁢(x-y),(x-y)〉, and A is a positive diagonal d×d matrix.

3.1 General Definitions

CP and Tucker rank-structured tensor formats have been applied for the quantitative analysis of correlation in multidimensional experimental data for a long time in chemometrics and signal processing [59, 8]. The Tucker tensor format was introduced in 1966 for tensor decomposition of multidimensional arrays in chemometrics [65]. Though the canonical representation of multivariate functions was introduced as early as in 1927 [29], only the Tucker tensor format provides a stable algorithm for decomposition of full-size tensors. A mathematical approval of the Tucker decomposition algorithm was presented in papers on higher-order singular value decomposition (HOSVD) and the Tucker ALS algorithm for orthogonal Tucker approximation of higher-order tensors [10]. For higher dimensions, the so-called Matrix Product States (MPS) (see the survey paper [57]) or the Tensor Train (TT) [53] decompositions can be applied. However, for three-dimensional applications, the Tucker and CP tensor formats remain the best choices. The fast convergence of the Tucker decomposition was proved and demonstrated numerically for higher-order tensors that arise from the discretization of linear operators and functions in ℝd for a class of function-related tensors and Green’s kernels in particular, it was found that the approximation error of the Tucker decomposition decayed exponentially in the Tucker rank [33, 36].

These results inspired the canonical-to-Tucker (C2T) and Tucker-to-canonical (T2C) decompositions for function-related tensors in the case of large input ranks, as well as the multigrid Tucker approximation [37].

A tensor of order d in a full format is defined as a multidimensional array over a d-tuple index set:

Here, 𝐀 is an element of the linear space

equipped with the Euclidean scalar product 〈⋅,⋅〉:𝕍n×𝕍n→ℝ, defined as

Tensors with all dimensions having equal size nℓ=n, ℓ=1,…,d, are called n⊗d tensors. The required storage size scales exponentially with the dimension, nd, which results in the so-called “curse of dimensionality”.

To avoid exponential scaling in the dimension, the rank-structured separable representations (approximations) of the multidimensional tensors can be used. The simplest separable element is given by the rank-1 tensor

with entries ui1,…,id=ui1(1)⁢⋯⁢uid(d), which requires only n1+…+nd numbers for storage.

The rank-1 canonical tensor is a discrete counterpart of the separable d-variate function, which can be represented as the product of univariate functions

An example of the separable d-variate function is f⁢(x1,x2,x3)=e(x1+x2+x3). Then, by discretization of this multivariate function on a tensor grid in a computational box, we obtain a canonical rank-1 tensor.

A tensor in the R-term canonical format is defined by a finite sum of rank-1 tensors (Figure 2, left)

where 𝐮k(ℓ)∈ℝnℓ are normalized vectors, and R is the canonical rank. The storage cost of this parametrization is bounded by dRn. An element a⁢(i1,…,id) of the tensor 𝐀=∑i=1R⊗ν=1dui⁢ν can be computed as

An alternative (contracted product) notation is used in computer science community:

where 𝐂=diag⁡{c1,…,cd}∈ℝR⊗d, and U(ℓ)=[𝐮1(ℓ)⁢…⁢𝐮R(ℓ)]∈ℝnℓ×R. An analogous multivariate function can be represented by a sum of univariate functions

For d≥3, there are no stable algorithms to compute the canonical rank of a tensor 𝐀, that is, the minimal number R in representation (3.2), and the respective decomposition with the polynomial cost in d, i.e., the computation of the canonical decomposition is an N-P hard problem [27].

Figure 2

Canonical (left) and Tucker (right) decompositions of three-dimensional tensors.

The Tucker tensor format (Figure 2, right) is suitable for stable numerical decompositions with a fixed truncation threshold. We say that the tensor 𝐀 is represented in the rank-𝐫 orthogonal Tucker format with the rank parameter 𝐫=(r1,…,rd) if

where {𝐯νℓ(ℓ)}νℓ=1rℓ∈ℝnℓ represents a set of orthonormal vectors for ℓ=1,…,d, and 𝐁=[𝐁ν1,…,νd]∈ℝr1×⋯×rd is the Tucker core tensor. The storage cost for the Tucker tensor is bounded by d⁢r⁢n+rd, with r=|𝐫|:=maxℓ⁡rℓ. Using the orthogonal side matrices V(ℓ)=[v1(ℓ)⁢…⁢vrℓ(ℓ)] and contracted products, the Tucker tensor decomposition may be presented in the alternative notation

In the case d=2, the orthogonal Tucker decomposition is equivalent to the singular value decomposition (SVD) of a rectangular matrix.

3.2 Tucker Decomposition of Full Format Tensors

We use the following algorithm to compute the Tucker decomposition of the full format tensor. The most time-consuming part of the Tucker algorithm is higher-order singular value decomposition (HOSVD), the computation of the initial guess for matrices V(ℓ) using the SVD of the matrix unfolding A(ℓ), ℓ=1,2,3 (Figure 3), of the original tensor along each mode of a tensor [10]. Figure 3 illustrates the matrix unfolding of the full format tensor 𝐀 (see (3.1)) along the index set I1={1,…,n1}.

Figure 3

Unfolding of a three-dimensional tensor along the mode Iℓ with ℓ=1.

The second part of the algorithm is the ALS procedure. For every tensor mode, a “single-hole” tensor of reduced size is constructed by the mapping all of the modes of the original tensor except one into the subspaces V(ℓ). Then the subspace V(ℓ) for the current mode is updated by SVD of the unfolding of the “single hole” tensor for this mode. This alternates over all modes of the tensor, which are updated at the current iteration of ALS. The final step of the algorithm is computation of the core tensor by using the ultimate mapping matrices from ALS.

The numerical cost of Tucker decomposition for full size tensors is dominated by the initial guess, which is estimated as O⁢(nd+1) when all nℓ=n are equal , or O⁢(n4) for our three-dimensional case. This step restricts the available size of the tensor to be decomposed since, for conventional computers, the three-dimensional case nℓ>102 is the limiting case for SVD.

The multigrid Tucker algorithm for full size tensors allows the computational complexity to be linear in the full size of the tensor, O⁢(nd), i.e., O⁢(n3) for three-dimensional tensors [37]. It is computed on a sequence of diadically refined grids and is based on implementing the HOSVD only at the coarsest grid level, n0≪n. The initial guess for the ALS procedure is computed at each refined level by the interpolation of the dominating Tucker subspaces obtained from the previous coarser grid. In this way, at fine three-dimensional Cartesian grids, we need 29 only O⁢(n3) storage (to represent the initial full format tensor) to contract with the Tucker side matrices, obtained by the Tucker approximation via ALS on the previous grids.

3.3 Illustration of the Low-Rank Approximation Idea

In this subsection we describe a possible ways to find a low-rank tensor approximation of the Matérn covariance matrix C⁢(r,ν) (Figure 1). Let 𝑭d=⊗ν=1d𝑭ν be the d-dimensional Fourier transform, where 𝑭-d=⊗ν=1d𝑭ν-1 is its inverse and ⊗ denotes the Kronecker product. We assume that 𝑼⁢(ξ)=𝑭d⁢(C⁢(r,ν)) is known analytically and has a low-rank tensor approximation 𝐔=∑j=1r⊗ν=1d𝐮j⁢ν. Since the Fourier and inverse Fourier transformations do not change the Kronecker tensor rank of the argument [51], by applying the inverse Fourier, we obtain a low-rank representation of the covariance function by applying the inverse Fourier:

where 𝐮~ν⁢i:=𝑭ν-1⁢(𝐮ν⁢i).

3.4 Sinc Approximation of the Matérn Function

The Sinc method provides a constructive approximation of the multivariate functions in the form of a low-rank canonical representation. It can be also used for the theoretical proof and for the rank estimation. Methods for the separable approximation of the three-dimensional Newton kernel and many other spherically symmetric functions that use the Gaussian sums have been developed since the initial studies in chemical [5] and mathematical literature [63, 6, 23, 17]. Here, we use a tensor-decomposition approach for lattice-structured interaction potentials [32]. We also recall the grid-based method for a low-rank canonical representation of the spherically symmetric kernel function q⁢(∥x∥), where x∈ℝd, d=2,3,…, by its projection onto the set of piecewise constant basis functions; see [3] for the case of Newton and Yukawa kernels for x∈ℝ3.

Following the standard schemes, we introduce the uniform n×n×n rectangular Cartesian grid Ωn with mesh size h=2⁢bn (we assume that n is even) in the computational domain Ω=[-b,b]3. Let {ψ𝐢} be a set of tensor-product piecewise constant basis functions

for the 3-tuple index 𝐢=(i1,i2,i3)∈ℐ, ℐ=I1×I2×I3, with iℓ∈Iℓ={1,…,n}, where ℓ=1, 2, 3. The generating kernel q⁢(∥x∥) is discretized by its projection onto the basis set {ψ𝐢} in the form of a third-order tensor of size n×n×n, which is defined entry-wise as

The low-rank canonical decomposition of the third-order tensor 𝐐 is based on applying exponentially convergent sinc-quadratures to the integral representation of the function q⁢(p), p∈ℝ, in the form

specified by the weights a1⁢(t),a2⁢(t)>0. Figure 4 illustrates a scheme of the proof of existence of the canonical low-rank tensor approximation. It could be easier to apply the Laplace transform to the Fourier transform of a Matérn covariance matrix than to the Matérn covariance. To approximate the resulting Laplace integral we apply the sinc quadrature. The number (2⁢M+1) of terms in the approximate sum (3.5) is the canonical tensor rank.

Figure 4

Scheme of the proof of existence of low-rank tensor approximation, r=2⁢M+1.

In particular, the sinc-quadrature for the Laplace–Gauss transform

can be applied, where the quadrature points (tk) and weights (ak) are given by

Under the assumption 0<a0≤|p|<∞, this quadrature can be proven to provide an exponential convergence rate in M (uniformly in p) for a class of functions a⁢(z) that are analytic in a certain strip |z|≤D of the complex plane such that the functions a1⁢(t)⁢e-p2⁢a2⁢(t) decay polynomially or exponentially on the real axis. The exponential convergence of the sinc-approximation in the number of terms (i.e., the canonical rank R=2⁢M+1) was analyzed elsewhere [63, 6, 23].

We assume that a representation similar to (3.5) exists for any fixed x=(x1,x2,x3)∈ℝ3 such that ∥x∥>a0>0. Then we apply the sinc-quadrature approximation (3.5) and (3.6) to obtain the separable expansion

providing an exponential convergence rate in M:

By combining (3.4) and (3.7), and taking into account the separability of the Gaussian basis functions, we arrive at the low-rank approximation of each entry of the tensor 𝐐:

Recalling that ak>0, we define the vector 𝐪 as

Then the third-order tensor 𝐐 can be approximated by the R-term (R=2⁢M+1) canonical representation

where 𝐪k(ℓ)∈ℝn. Given a threshold ε>0, M can be chosen as the minimal number such that in the max-norm

The skeleton vectors can be reindex by k↦k′=k+M+1, 𝐪k(ℓ)↦𝐪k′(ℓ) (k′=1,…,R=2⁢M+1), ℓ=1,2,3. The symmetric canonical tensor 𝐐R∈ℝn×n×n in (3.8) approximates the three-dimensional symmetric kernel function q⁢(∥x∥) (x∈Ω), centered at the origin, such that 𝐪k′(1)=𝐪k′(2)=𝐪k′(3) (k′=1,…,R).

In some applications, the tensor can be given in the canonical tensor format, but with large rank R and discretized on large grids n×n×…×n; thus, computation of the initial of guess in the Tucker-ALS decomposition algorithm becomes intractable. This situation may arise when composing the tensor approximation of complicated kernel functions from simple radial functions that can be represented in the low-rank CP format.

For such cases, the canonical-to-Tucker decomposition algorithm was introduced [37]. It is based on the minimization ALS procedure, similar to the Tucker algorithm, described in Section 3.2 for full-size tensors, but the initial guess is computed by just the SVD of the side matrices U(ℓ)=[𝐮1(ℓ)⁢…⁢𝐮R(ℓ)]∈ℝnℓ×R, ℓ=1,2,3; see (3.3). This schema is the Reduced HOSVD (RHOSVD), which does not require unfolding of the full tensor.

Another efficient rank-structured representation of the multidimensional tensors is the mixed-tensor format [31], which combines either the canonical-to-Tucker decomposition with the Tucker-to-canonical decomposition, or standard Tucker decomposition with the canonical-to-Tucker decomposition, in order to produce a canonical tensor from a full-size tensor.

3.5 Laplace Transform of the Covariance Matrix

The integral representations like (3.5) can be derived by the Laplace transform either directly to the Matérn covariance function or to its spectral density (2.3). For example, in the case of the Newton kernel, q⁢(p)=1p, and the Laplace–Gauss transform representation takes the form

In this case, 𝐪k(ℓ)=𝐪-k(ℓ), and the sum of (3.8) reduces to k=0,1,…,M, implying that R=M+1. Therefore, the Laplace transform representation of the Slater function q⁢(p)=e-2⁢α⁢p (i.e., exponential covariance) with p=∥x∥2 can be written as

When the Matérn spectral density in (2.3) has an even dimension parameter d=2⁢d1, d1=1,2,… and ν=0,1,2,…, the Laplace transform

can be applied after substituting p=∥ξ∥2 in

If -η=ν+d1=12,32,52,…,2⁢k+12,…, then the Laplace transform is

3.6 Covariance Matrix in Rank-Structured Tensor Format

In what follows, we consider the CP approximation of the radial function q⁢(r) in the positive sector, i.e., on the domain [0,b]3 (by symmetry, the canonical tensor can be extended to the whole computational domain [-b,b]3.) Let the covariance function q⁢(r)=Cν,ℓ⁢(r) in (2.2) be represented by the rank-R symmetric CP tensor on an n×n×n tensor grid, denoted by Ωn⊂Ω=[0,b]3 as described in the previous sections,

with the same skeleton vectors 𝐪k(ℓ)∈ℝn for ℓ=1,2,3.

We define the covariance matrix 𝐂=[c𝐢,𝐣]∈ℝ𝐧×𝐧 entry-wise by

Using the tensor representation (3.9), we represent the large n3×n3 matrix 𝐂 in the rank-R Kronecker (tensor) format as

where the symmetric Toeplitz matrix 𝐐k(ℓ)=Toepl⁡[𝐪k(ℓ)]∈ℝn×n, ℓ=1,2,3, is defined by its first column, which is specified by the skeleton vectors 𝐪k(1)=𝐪k(2)=𝐪k(3) in the decomposition (3.9).

Figure 5 illustrates eight selected canonical generating vectors from 𝐪k(1) for k=1,…,R, R=34, on a grid of size n=2,049 for the Slater function e-∥x∥p, which defines the corresponding Toeplitz matrices.

Figure 5

Selected eight canonical vectors from the full set 𝐪k(1), k=1,…,R, see (3.9).

A Toeplitz matrix can be multiplied by a vector in O⁢(n⁢log⁡n) operations via complementing it to the circulant matrix. In general, the inverse of a Toeplitz matrix cannot be calculated in the closed form; unlike circulant matrices, which can be diagonalized by the Fourier transform.

In the rest of this subsection, we introduce the numerical scheme, based on certain specific properties of the skeleton matrices in the symmetric rank-R decomposition (3.10) that is capable of rank-structured calculations of the analytic matrix functions ℱ⁢(𝐂R). We discuss the most interesting examples of the functions ℱ1⁢(𝐂)=𝐂-1 and ℱ2⁢(𝐂)=𝐂12, where 𝐂=𝐂R.

Given an symmetric positive definite matrix such that ∥𝐀∥=q<1, the matrix-valued function is given as the exponentially fast converging series

where the matrix 𝐀, acting in the multidimensional index set, allows the low-rank Kronecker tensor decomposition. Then the low-rank tensor approximation of ℱ⁢(𝐀) can be computed by the “add-and-compress” scheme, where each term in the series above (3.11) is summed using the rank-truncation algorithm in the corresponding format.

To limit the rank-structured evaluation of ℱ1⁢(𝐂) to the described framework, we propose the special rank-structured additive splitting of the covariance matrix 𝐂 with the easily invertible dominating part. To that end, we construct the diagonal matrix 𝐐0(1) by assembling all of the diagonal sub-matrices in 𝐐k(1), k=1,…,R (in the following, we simplify the notation by omitting the upper index (1)):

and modify each Toeplitz matrix 𝐐k by subtracting its diagonal part,

Using the matrices defined above, we introduce the rank R+1 additive splitting of 𝐂 , which is defined by the skeleton matrices 𝐐0 and 𝐐^k since 𝐂 is Kronecker symmetrical,

Hence, we have

Likewise, since 𝐐0⊗𝐐0⊗𝐐0 is a scaled identity, we obtain

We assume that ∥𝐐0-1⁢𝐐^k(1)∥<1 for k=1,…,R in some norm; thus, we can apply the “add-and-compress” scheme described above.

We illustrate the “add-and-compress” computational scheme in the following example. We consider the covariance matrix 𝐂R, obtained by a rank-Rsinc approximation of the Slater function e-∥x∥p with R=40 on a grid with n=1,025 sampling points. Figure 6 demonstrates the decay in both the matrix norms 𝐐k(1) (left), and the scaled, preconditioned matrices 𝐐0-1⁢𝐐^k(1), k=1,…,R (right). We use the scaling factor of 1n. The right figure indicates that the analytic matrix functions ℱ⁢(𝐂R) can be evaluated by using an exponentially fast convergent power series supported by the “add-and-compress” strategy to control the tensor rank.

Figure 6

Scaled norms ∥𝐐k(1)∥ (left) and ∥𝐐0-1⁢𝐐^k(1)∥ (right) vs. k=1,…,R.

In the Kriging calculations (Task 3 above), the low-rank tensor structure in the covariance matrix 𝐂R can be directly utilized if the sampling points in the Kriging algorithm form a smaller m1×m2×m3 tensor sub-grid of the initial n×n×n tensor grid Ωn with mℓ<n. The same argument also applies to the evaluation of conditional covariance. In the general case of “non-tensor” locations of the sampling points, some mixed tensor factorizations could be applied.

3.7 Numerical Illustrations

In what follows, we check some examples of the low-rank Tucker tensor approximation of the p-Slater function C⁢(x)=e-∥x∥p, and Matérn kernels with full-grid tensor representation. We demonstrate the fast exponential convergence of the tensor approximation in the Tucker rank. The functions were sampled on the n1×n2×n3 three-dimensional Cartesian grid with nℓ=100, ℓ=1,2,3.

For a continuous function q:Ω→ℝ, where Ω:=∏ℓ=1d[-bℓ,bℓ]⊂ℝd, and 0<bℓ<∞, we introduce the collocation-type function-related tensor of order d:

where (xi1(1),…,xid(d))∈ℝd are grid collocation points, indexed by ℐ=I1×…×Id,

which are the nodes of equally spaced subintervals with a mesh size of hℓ=2⁢bℓnℓ-1.

We test the convergence of the error in the relative Frobenius norm with respect to the Tucker rank for p-Slater functions with p=0.1,0.2,1.9,2.0. The Frobenius norm is computed as

where 𝐐(r) is the tensor reconstructed from the Tucker rank-r decomposition of 𝐐.

Figure 7

Convergence in the Frobenius error (3.13) with respect to the Tucker rank for the function (3.14) withp=0.1,0.2,1.9,2.0 (left); Decay of singular values of the weighted Slater function (right).

Figure 8

Cross section of the three-dimensional radial function (3.14) with p=0.1 (left) and p=1.9 (right) at level z=0.

Figure 9

Multigrid Tucker: convergence with respect to Tucker ranks of a Slater function with p=1 on a sequence of grids.

Figure 10

Convergence with respect to the Tucker rank of three-dimensional spectral density of Matérn covariance (3.15)with α=0.1 (left) and α=100 (right).

Figure 7 shows convergence of the Frobenius error with respect to the Tucker rank for the p-Slater function discretized on the three-dimensional Cartesian grid

for different values of the parameter p. These functions are illustrated in Figure 8 for p=0.1 and p=1.9. Figure 9 shows for a Slater function with p=1 the dependence of the Tucker decomposition error in the Frobenius norm (3.13) on the Tucker rank, for the increasing grid parameter n.

Figure 10 shows the convergence with respect to the Tucker rank for the spectral density of Matérn covariance

where α∈(0.1,100) and d=1,2,3. The Tucker decomposition rank is strongly dependent on the parameter α and weakly depend on the parameter ν. The three-dimensional Matérn functions with the parameters ν=0.4, α=0.1 (left) and ν=0.4, α=100 (right) are presented in Figure 11, showing the function at z=0.

These numerical experiments demonstrate the good algebraic separability of the typical multidimensional functions used in spatial statistics, that lead us to apply the low-rank tensor decomposition methods to the multidimensional problems of statistical data analysis.

Figure 11

The shape of three-dimensional spectral density of Matérn covariance (3.15) with α=0.1 (left) and α=100 (right).

Computing Matrix-Vector Product.

Let 𝐂=∑i=1r⊗μ=1d𝐂i⁢μ. If 𝐳 is separable, i.e., ∥𝐳-∑j=1rb⊗ν=1d𝐳j⁢ν∥≤ε, then

If 𝐳 is non-separable, then low-rank tensor properties cannot be employed and either the FFT idea [51] or the hierarchical matrix technique should be applied instead [38, 21, 19, 22].

Trace and Diagonal of 𝐂.

Let 𝐂≈𝐂~=∑i=1r⊗μ=1d𝐂i⁢μ. Then