{'fraction_non_alphanumeric': 0.059814179580438866, 'fraction_numerical': 0.013255649595780256, 'mean_word_length': 3.604285714285714, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 38, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High Dimensional Multi-Level Covariance Estimation and Kriging

With the advent of big data sets much of the computational science and engineering communities have been moving toward data-driven approaches to regression and classification. However, they present a significant challenge due to the increasing size, complexity and dimensionality of the problems. In this paper a multi-level kriging method that scales well with dimensions is developed. A multi-level basis is constructed that is adapted to a random projection tree (or kD-tree) partitioning of the observations and a sparse grid approximation. This approach identifies the high dimensional underlying phenomena from the noise in an accurate and numerically stable manner. Furthermore, numerically unstable covariance matrices are transformed into well conditioned multi-level matrices without compromising accuracy. A-posteriori error estimates are derived, such as the sub-exponential decay of the coefficients of the multi-level covariance matrix. The multi-level method is tested on numerically unstable problems of up to 50 dimensions. Accurate solutions with feasible computational cost are obtained.


Introduction
Big data sets arise from many fields, including, but not limited to commerce, astrophysical sky-surveys, enviromental data, and tsunami warning systems. With the advent of big data sets much of the computational science and engineering communities have been moving toward data-driven approaches to regression and classification. These approaches are effective since the underlying data is incorporated. However, they present a numerical challenge due to increasing size, complexity and dimensionality.
Due to the high dimensionality of the underlying data many modern machine learning methods, such as classification and regression algorithms, seek a balance between accuracy and computational complexity. How efficient this balance is depends on the approach. Linear methods are fast, but only work well when there is linear separation of the data.
For non-linear description of the data, kernel approaches have been effective under certain circumstances. These methods rely on Tikhonov regularization of the data to obtain a functional representation, where it is assumed that the noise model of the phenomena is known. However, this assumption is not necessarily satisfied in practice and can lead to significant errors as the algorithm cannot distinguish between noise and the underlying phenomena.
Sparse grid polynomial approximation [3,4,6,7,10,[26][27][28]31] has emerged as one of the most appealing methods for approximating high-dimensional functions. It is simple to use and converges significantly faster if the function at hand presents some degree of regularity. In [30] the authors develop a sparse grid approach to high dimension classification. However, the noise is assumed to be known, which can still lead to significant error.
To incorporate the variability of the noise model a class of machine learning algorithms based on Bayes method have been developed. In this approach the noise model is assumed to be known up to a class of probability distributions and an optimal choice is made that fits the training data and noise. For example, from the Geo-statistics community a well known approach to identifying the underlying data and noise model is known as kriging [25]. The noise model parameters are estimated from the Maximum Likelihood Estimation (MLE) of the likelihood function.
Kriging methods are effective in separating the underlying phenomena from the noise model. In practice the covariance matrices tend to be ill-conditioned with increasing number of observations making kriging methods numerically fragile. Moreover, applications are limited to 2 or 3 dimensions as the computational cost increases significantly with dimension.
In [8] a novel algorithm to solve kriging problems is proposed. The method is fast and robust. In particular, it can solve kriging problems that where not tractable with previous methods. However, this approach is limited to 2 or 3 dimensions and the computational cost scales very fast with spatial dimension, thus making it impractical for high dimensional problems. Other kriging methods have been developed using skeletonization factorizations [24], low-rank [29] and Hierarchical Matrices [19,21] approaches.
In this paper we extend the kriging approach in [8] to high dimensions. This novel approach effectively identifies the high dimensional underlying phenomena from the noise in a numerically stable manner. Our approach transforms high dimensional ill-conditioned covariance matrices to numerically stable multi-level covariance matrices without compromising accuracy. In Section 2 the problem formulation is introduced. In Section 3 Multivariate polynomial approximations and sparse grids are discussed. In section 4 it is shown how to construct a multi-level basis with a random projection that scales well with dimension. In section 5 the construction of the multi-level covariance matrix is discussed. In section 6 the multilevel estimator and predictor are formulated. In section 7 error estimates for the decay of the multi-level covariance matrix and the inverse solution are derived. These estimates are based on the Smolyak sparse grid formulas.
Note that for the less mathematically inclined reader this section can be skipped. In section 8 the multi-level kriging method is tested on numerically unstable high dimensional problems. The results are obtained with good accuracy and feasable computational cost.

Problem setup
Consider the following model for a Gaussian random field Z: where d is the number of spatial dimensions, k : R d → R p is a functional vector of the spatial location x, β ∈ R p is an unknown vector of coefficients, and ε is a stationary mean zero Gaussian random field with parametric covariance function C(x, x ′ ; θ) = cov{ε(x), ε(x ′ )} with an unknown vector of positive parameters θ ∈ R d . Suppose that we obtain N observations and stack them in the data vector Z = (Z(x 1 ), . . . , Z(x N )) T from locations S := {x 1 , . . . , x N }, where x 1 = x 2 = x 3 = · · · = x N −1 = x N . Let C(θ) = cov(Z, Z T ) ∈ R N ×N be the covariance matrix of Z and assume it is positive definite for all θ ∈ R w . Define M = k(x 1 ) . . . k(x N ) T ∈ R n×p and assume it is of full rank, p. The model (1) leads to the vectorial formulation where ε is a Gaussian random vector, ε ∼ N (0, C(θ)). The aim now is to: • Estimate the unknown vectors β and θ; • Predict Z(x 0 ), where x 0 is a new spatial location. These two tasks are particularly computationally challenging when the sample size N and number of dimensions d are large.
The unknown vectors β and θ are estimated with the log-likelihood function which can be profiled by generalized least squares witĥ A consequence of profiling is that the maximum likelihood estimator (MLE) of θ then tends to be biased. To address the prediction part, consider the best unbiased predictorẐ(x 0 ) = λ 0 + λ T Z where λ = (λ 1 , . . . , λ n ) T . The unbiased constraint implies λ 0 = 0 and M T λ = k(x 0 ). The minimization of the mean squared prediction error E[{Z(x 0 )−λ T Z} 2 ] under the constraint M T λ = k(x 0 ) yieldŝ where c(θ) = cov{Z, Z(x 0 )} ∈ R n andβ is defined in (4). Now, let α := (α 1 , . . . , α d ) ∈ Z d , |α| := α 1 + · · · + α d , To mitigate this problem we need to control the number of polynomial terms along each dimension while still retaining good accuracy by selecting a judicious set of polynomial functions k(x).
The second challenge is that the covariance matrix C(θ) in many practical cases is illconditioned, leading to slow and inaccurate estimates of θ. In this paper, we propose a new transformation of the data vector Z leading to a decoupled multi-level description of the model (1) without any loss of structure. This multi-level representation leads to significant computational benefits when computing the kriging predictorẐ(x 0 ) in (5) for large sample size N and high dimensions d.

Multivariate polynomial approximation
We now have to choose an appropriate model k : R d → R p of random field Z. The choice of k (and thus P p (S)) will determine the accuracy of the representation of the underlying deterministic function. Indeed, there are two main factors that we should consider. First, the size of p and second, the accuracy of the approximation in a suitable functional space. The balance between these two factors should be carefully considered.
ii) Restrict the number of polynomials along each dimension by using the set of monomials contained in i.e. k(x) is built from the monomials in Q Λ(w) .
The choice of index set Λ(w) will determine the accuracy of the polynomial representation and the size of p. We first consider the index set for a full Tensor Product (TP) grid This is a poor choice as p increases d i=1 (p i +1) and quickly becomes intractable. Another choice is the Total Degree (TD) index set, which grows as d + w w , but still suffers from the curse of dimensionality, thus impractical for high dimensions. Smolyak (SM) and Hyperbolic Cross (HC) index sets are a popular choice for large dimensions. SM have been used extensively in sparse grid approximations of high dimensional analytic functions [27,28]. HC have been used in the context of Fourier approximations of high dimensional functions. In Table 1 the different choices for Λ(w) are summarized.
In Figure 1 a comparison between TP,TD,SM and HC is shown. As observed, the SM and HC grows much slower with respect to the higher polynomial order.
Approx. space

Sparse grids
Consider the problem of approximating a function u : Γ → R on the domain Γ d . A good choice for approximating such functions are Smolyak sparse grids with an appropriate set of abscissas. Note that in this paper a sparse grid representation is not explicitly used, unless the observation nodes correspond exactly to the abscissa grid. Rather, it is used to estimate the accuracy of a sparse multi-level covariance matrix. For q ∈ N define the following spaces be the span of tensor product polynomials of degree at most p = (p 1 , . . . , p N ); i.e., P p (Γ d ) = d n=1 P pn (Γ n ) with P pn (Γ n ) := span(y m n , m = 0, . . . , p n ), n = 1, . . . , d. Consider the univariate Lagrange interpolant along the n th dimension of Γ d where I m(i) n k=1 is a Lagrange basis of the space P pn (Γ n ), i ≥ 0 is the level of approximation and m(i) ∈ N + is the number of collocation nodes at level i ∈ N + where m(0) = 0, m(1) = 1 and We can now construct an interpolant by taking tensor products of I m(i) n along each dimension n. However, as d increases the dimension of P p increases as d n=1 (p n + 1). Thus even for moderate d dimensions the computational cost of the Lagrange approximation becomes intractable. However, in the case of sufficient complex analytic regularity of the QoI with respect to the random variables defined on Γ d , the application of Smolyak sparse grids are better suited. In the rest of this section the construction of the classical Smolyak sparse grid (see e.g. [6,31]) is summarized. More details can be found in [4]. Given an integer w ≥ 0, called the approximation level, and a multi-index i = (i 1 , . . . , i d ) ∈ N d + , let g : N d + → N be a strictly increasing function in each argument and define a sparse grid approximation of function v(y) ∈ C 0 (Γ d ) or equivalently written as From the previous expression, the sparse grid approximation is obtained as a linear combination of full tensor product interpolations. However, the constraint g(i) ≤ w in (7) is typically chosen so as to forbid the use of tensor grids of high degree in all directions at the same time.
The following result proved in [4] states that the sparse approximation formula S m,g w,d is exact in P Λ m,g (w) (Γ d ): Here V denotes a Banach space defined on Γ d and The most typical choice of m and g is given by (see [6,31]) Other choices are shown in Table 2. This choice of m and g combined with the choice of Clenshaw-Curtis (extrema of Chebyshev polynomials) abscissas leads to nested sequences of one dimensional interpolation formulas and a sparse grid with a highly reduced number of nodes compared to the corresponding tensor grid. Another good choice includes the Chebyshev Gauss-Lobatto abscissas. For any choice of m(i) > 1 the Clenshaw-Curtis abscissas are given by Approx. space sparse grid: m, g polynomial space: Λ(w) Total Suppose u ∈ C k mix (Γ; R), which is the space of continuous mixed derivatives of degree k and is defined as The authors in [6] show that The weakness of this bound is that the coefficient C(k, d) is in general not known, thus making it difficult to estimate a-posteriori bounds on the solution. A better choice is the L ∞ (Γ d ) bounds derived in [28]. These bounds are explicit, and the coefficients are directly based on analytic extensions of v(y) on a well defined region Θ ⊂ C d . In [27,28] the authors derive error estimates for isotropic and anisotropic Smolyak sparse grids with Clenshaw-Curtis and Gaussian abscissas, where the error v − S m,g w,d [v] L ∞ (Γ) exhibits algebraic or sub-exponential convergence with respect to the number of collocation knots η (See Theorems 3.10, 3.11, 3.18 and 3.19 in [28] for more details). However, for these estimates to be valid v ∈ C 0 (Γ, R) has to admit an extension in the region defined by the polyellipse in C d , Eσ 1 ,...,σ d := Π d n=1 E n,σn ⊂ Θ, where E n,σn = z ∈ C,σ n ≥ δ n ≥ 0; Re z = e δn + e −δn 2 cos(θ), Im z = e δn − e −δn 2 sin(θ), θ ∈ [0, 2π) , andσ n > 0 (see Figure 2).

Remark 1.
For simplicity, in [28], the error bounds assume that the term M(v) is bounded by  one. We remove this assumption and modify Theorem 3.10 and 3.11 in [28] to show the error bounds with the termM (v), where We now introduce Chebyshev polynomials. These will be useful in deriving error estimates of the sparse grid. Let T k : Γ 1 → R, k = 0, 1, . . . , be a Chebyshev polynomial over [−1, 1], which are defined recursively as: The following theorem describe the approximation of analytic functions with Chebyshev polynomials.
Theorem 1. Let q be analytic and bounded by M on E log ρ , ρ > 1, then the expansion holds for all y ∈ E log ρ and |α k | ≤ M/ρ k . Furthermore if y ∈ [−1, 1] then Since the interpolation operator I m(i) n is exact on the space P pn−1 , then if v is analytic on Γ n we have from Theorem 1 that where Λ m(i) is the Lebesgue constant and is bounded by 2π −1 (log (m − 1) + 1) (see [3]). Thus, for n = 1, . . . , d, with ρ n = logσ n we have where σ n =σ n 2 > 0 and C(σ n ) := 2 (e σn −1) . We have that for all n = 1, . . . , d and set σ ≡ min n=1,...,d σ n , i.e. for an isotropic sparse grid the overall asymptotic sub-exponential decay rateσ will be dominated by the smallestσ n .

Remark 2.
As the reader might have realized, it is not necessary to show that there is a bounded analytic extension on Eσ 1 ,...,σ d . Since the convergence rate of the isotropic sparse grid is controlled byσ, it is sufficient to show that v is bounded in Eσ ,...,σ := d n=1 E n,σ and the boundM (v) will be smaller. However, an analytic extension on Eσ 1 ,...,σ d will be relevant for anisotropic sparse grids [27], which leads to faster convergence rates. We leave this as a future avenue to explore. i n and Λ m.g (w) := {p ∈ N d : n (p n + 1) ≤ w}. The cardinality of Λ m.g (w) for the HC grows much slower than SM or TD and works very well in practice. However, to my knowledge no L ∞ (Γ d ) error bounds exist.
Remark 5. For many practical cases all the dimensions of Γ d are not equally important. In this case the dimensionality of the sparse grid can be significantly reduced by means of a dimension adaptive [17] or anisotropic sparse grid. It is straightforward to build related anisotropic sparse approximation formulas by making the function g act differently on the input random variables y n . Anisotropic sparse stochastic collocation [27] combines the advantages of isotropic sparse collocation with those of anisotropic full tensor product collocation.
Another choice is to build quasi-optimal isotropic/anisotropic sparse grids [7,26]. These grids have provable exponential convergence rates. However, in general, it is not shown how to synthesize quasi-optimal grids.

Multi-level approach
The general approach of this paper and multi-level basis construction are now presented. Denote by P p (S) the span of the columns of the design matrix M. Let L ∈ R p×N be an orthogonal projection from R n to P p (S) and W ∈ R (N −p)×N be an orthogonal projection from R n to P p (S) ⊥ (the orthogonal complement of P p (S)). Moreover assume that the operator W L is orthonormal.
By applying the operator W to (2) we obtain Z W = WZ = W(Mβ + ε) = Wε. Our first observation is that the trend contribution Mβ is filtered out from the data Z. We can now formulate the estimation of the covariance parameters θ without the trend. The new loglikelihood function becomes where C W (θ) = WC(θ)W T and Z W ∼ N N −p (0, WC(θ)W T ). A consequence of the filtering is that we obtain an unbiased estimator. The decoupling of the likelihood function is not the only advantage of using C W (θ). The following theorem also shows that C W (θ) is more numerically stable than C(θ).

Proof. To see this let
Now, it is not hard to see that Theorem 3 states that the condition number of C W (θ) is less or equal to the condition number of C(θ). Thus computing the inverse of C W (θ) (using a direct or iterative method) will generally be more stable.
In practice, computing the inverse of C W (θ) will be much more stable than C(θ) depending on the choice of the index set Λ(w). This has many significant implications as it will now be possible to solve numerically unstable problems.
There are other advantages to the structure of the matrix C W (θ). In section 7 we show that for a good choice of the P(S) the entries of C W (θ) decay rapidly, and most of the entries can be safely eliminated without loosing accuracy. A level dependent criterion approach is shown in section 5 that indicate which entries are computed and which ones are not. With this approach a sparse covariance matrixC W is constructed such that it is close to C W in a matrix norm sense, even if the observations are highly correlated with distance.

Random projection multi-Level basis
In this section the construction of Multi-Level Basis (MB) is shown. The MB can then be used to: (i) form the multi-level likelihood (12), (ii) sparsify the covariance matrix C W (θ), and (iii) improve the conditioning over the covariance matrix C(θ). But first, let us establish notations and definitions: • Given Q Λ m,g (w) and the locations S construct the design matrix M. Furthermore, form a second set of monomialsQ a Λ m,g (w) := Q Λ m,g (w+a) for a = 0, 1, . . . , i.e. Q Λ m,g (w) ⊂Q a Λ m,g (w) . Denote the accuracy parameterp ∈ N as the cardinality ofQ a Λ m,g (w) . From the set of monomialsQ a Λ m,g (w) , for some user given parameter a ∈ N 0 , and the set of observations S generate the design matrixM a . Denote also the space Pp(S) as the span of the columns ofM a .
where Γ is the gamma function, 0 < ν, 0 < ρ < ∞, and K ν is the modified Bessel function of the second kind.
Remark 6. The Matérn covariance function is a good choice for the random field model. The parameter ρ controls the length correlation and the parameter ν changes the shape. For example, if ν = 1/2 + n, where n ∈ N + , then (see [2]) The first step is to decompose the domain Γ d into a multi-level domain decomposition. A first choice that comes to mind is a kD-tree decomposition of the space R d . This is a good choice for lower dimensions, however, as the number of dimensions d become larger a better approach is to use a Random Projection (RP) tree [12].
First start with the root node B 0 0 at level 0 that contains all the observation nodes in S. Now, split these nodes into two children cells B 1 1 and B 1 2 at level 1 according to the following rule: Project all the nodes x ∈ S in the cell onto the unit vector v.

Split the cell with respect to the median of the projections.
For each cell B 1 1 and B 1 2 repeat the procedure until there is at most p nodes at the leaf nodes. Thus a binary tree is obtained, which is of the form where t is the maximal depth (level) of the tree. Now, let B be the set of all the cells in the tree and B n be the set of all the cells at level 0 ≤ n ≤ t. The construction of the RP tree is described in Algorithms 1 and 2. In addition, for each cell a unique node number, current tree depth, threshold level and projection vector are also assigned. This will be useful for searching the tree.

Remark 7.
A kD-tree can also be constructed with Algorithms 1 and 3. Now, suppose there is a one-to-one mapping between the set of unit vectors E := {e 1 , . . . , e N }, which is denoted as leaf unit vectors, and the set of locations {x 1 , . . . , x N }, i.e. x n ←→ e n for all n = 1, . . . , N. It is clear that the span of the vectors {e 1 , . . . , e N } is R N . The next step is to construct a new basis of R n that is multi-level and orthonormal.
(a) Start at the maximum level of the random projection tree, i.e. q = t.
(b) For each leaf cell B q k ∈ B q assume without loss of generality that there are s observations i,j , d q i,j ∈ R and for some a ∈ N + . It is desired that the new discrete MB vector ψ q j be orthogonal to Pp(S), i.e., for all g ∈ Pp(S): Figure 3: Random projection tree construction.
Input: S, node, currentdepth, n 0 Output: Tree, node begin if Tree = root then node ← 0, currentdepth ← 0 Tree ← MakeTree(S, node, currentdepth + 1, n 0 ) else Tree.node = node Tree.currentdepth = currentdepth -1 for all e i ∈ C q k . Now, suppose that the matrix M q has rank a and then perform the Singular Value Decomposition (SVD). Denote by UDV the SVD of M q , where U ∈ Rp ×s , D ∈ R s×s , and V ∈ R s×s .
iii) To satisfy equation (13) make the following choice where the columns a + 1, . . . , s form an orthonormal basis of the nullspace N 0 (M). Similarly, the columns 1, . . . , a form an orthonormal basis of R s \N 0 (M). Since the vectors in C q k are orthonormal then φ q 1 , . . . , φ q a , ψ q a+1 , . . . , ψ q form an orthonormal basis of S. Moreover ψ q a+1 , . . . , ψ q satisfy equation (13), i.e., are orthogonal to Pp(S) and are locally adapted to the locations contained in the cell B q k .
Input: S Output: Rule, threshold, v begin choose a coordinate direction Rule(x) := x · v ≤ threshold = median end Algorithm 3: ChooseRule(S) function for kD-tree iv) Denote by D q k the collection of all the vectors ψ q a+1 , . . . , ψ q . Notice that the vectors φ q 1 , . . . , φ q a , which are denoted with a slight abuse of notation as the scaling vectors, are not orthogonal to Pp(S). They need to be further processed.
v) Let D q be the union of the vectors in D q k for all the cells B q k ∈ B q k . Denote by W q (S) as the span of all the vectors in D q .
(c) The next step is to go to level q − 1. For any two sibling cells denote B q lef t and B q right at level q denote C q−1 k as the collection of the scaling functions from both cells, for some indexk.
is not a leaf cell, then repeat steps (b) to (d), but replace the leaf unit vectors with the scaling vectors contained in The ML basis vectors will span the space When the algorithm stops a series orthogonal subspaces W −1 (S), . . . , W t (S) (and their corresponding basis vectors) are obtained. These subspaces are orthogonal to V −1 (S) := span{φ −1 1 , . . . , φ −1 p }. Note that the orthonormal basis vectors of V −1 (S) also span the space P p (S).
Remark 8. Following Lemma 2 in [9] it can be shown that where W −1 (S) := Pp(S)\P p (S). Also, it can then be shown that at most O(Nt) computational steps are needed to construct the multi-level basis of R N .

From the basis vectors of the subspaces (P
since the total number of orthonormal vectors that span P p (S) is p. Conversely, the total number of orthonormal vectors that span (P p (S)) ⊥ is N − p.
Let L be a matrix where each row is an orthonormal basis vector of P p (S). For i = −1, . . . , t let W i be a matrix where each row is a basis vector of the space T , therefore a) The matrix P := W L is orthonormal, i.e., PP T = I; b) Any vector v ∈ R n can be written as The matrix W contains at most O(Nt) non-zero entries and L contains at most O(Np) non-zero entries. This implies that for any vector v ∈ R n the computational cost of applying Wv is at most O(Nt) and Lv is at most O(Np). Lemma 1. Assuming that n 0 = 2p, for any level q = 0, . . . , t there is at mostp2 q multi-level basis vectors. For level q = −1 there is at most p −p multi-level vectors.
Proof. By construction any leaf cell has at mostp multi-level vectors andp scaling vectors that are to be used for the next level.
Starting at the finest level t, for each cell B t k ∈ B t there is at mostp multi-level vectors and p scaling vectors that are to be used for the next level. Since there is at most 2 t cells then there is at most 2 tp multi-level vectors. Now, for each pair of left and right (siblings) cells at level t the parent cell at level t − 1 will have at most 2p scaling functions. Thus at mostp multi-level vectors andp scaling vectors are obtained that are to be used for the next level. Now, the rest of the cells at level t are leafs and will have at mostp multi-level vectors andp scaling vectors that are to be used for the next level.
Since there is at most 2 t−1 cells at level t − 1, there is at most 2 t−1p multi-level vectors. Now, follow an inductive argument until q = 0. Finally at level q = −1 there will be at most p −p Lemma 2. For any level q = 0, . . . , t any multi-level vector ψ q m associated with a cell B q k ∈ B q has at most 2 t−q+1p non zero entries.
Proof. For any leaf cell at the bottom of the tree (level t) there is at most 2p observations, thus the number of non zero entries of level q multi-level vectors is 2p. Combining the left and right cells, the parent cell has at most 4p observations, thus the associated multi-level vectors has 4p non zero entries. By induction at any level l the number of nonzero entries is at most 2 t−q+1p . Now for any leaf cell at any other level l < t the number of nonzero entries is at most 2p. Following an inductive argument the result is obtained.
Remark 9. Extended multi-level basis: A variation to the multi-level basis is obtained by replacingQ a Λ m,g (w) , a ∈ N + , with an extended set of monomials. LetΛ m,g (w), w ∈ N + , be an index set such thatΛ m,g (w) := {p ∈ N d + , 2p ∈ N d + : p ∈ Λ m,g (w)}, and any repeated index inΛ m,g (w) is removed. From this index set form the set of monomials Q ã Λ m,g (w) , for all a ∈ N + . As defined before, the accuracy parameterp ∈ N is the cardinality of Q ã Λ m,g (w) and the design matrixM a is formed fromQ ã Λ m,g (w) and S.

Multi-level covariance matrix
The multi-level covariance matrix C W (θ) and sparse versionC W (θ) can be now constructed.
Recall from the discussion in Section 4 that C W (θ) := WC(θ)W T . From the multi-level basis construct in Section 4.1 the following operator is built: j are formed from all the interactions of the MB vectors between levels i and j, for all i, j = −1, . . . , t. The structure of C W (θ) is shown in Figure 4. Thus for any ψ ĩ l and ψ j k vectors there is a unique entry of C i,j W of the form (ψ ĩ k ) T C(θ)ψ j l . The blocks C i,j W , where i = −1 or j = −1, correspond to the case where the accuracy termp > p.

Sparsification of multi-level covariance matrix
A sparse version of the covariance matrix C W (θ) can be built by using a level and distance dependent strategy: 1. Given a cell B l k at level l ≥ −1 identify the corresponding tree node value Tree.node and the tree depth Tree.currentdepth. Note that the Tree.currentdepth and the MB level q are the same for q = 0, . . . , t. However, for q = −1 the MB is associated to the Tree.currentdepth = 0.
2. Let K ⊂ S be all the observations nodes contained in the cell B l k . 3. Let τ ≥ 0 be the distance parameter given by the user. 4. Let the Targetdepth be equal to the desired level of the tree. In the case that it is −1 then the Targetdepth is zero.
The objective now is to find all the cells at the Targetdepth that overlap a ball with radius τ from the projection of each of the nodes x ∈ K onto the vector Tree.v. This is done by searching the tree from the root node. At each traversed node check that all the nodes x ∈ K satisfy the following rule: x · Tree.v + τ ≤ T ree.threshold.
If this is true then the search continues down the left tree, otherwise both the left and right tree are searched. The full search algorithm is described in Algorithms 4, 5, and 6. The sparse matrix blocks C i,j W (θ) can be built with a simple rule from Algorithm 6. Compute all the entries of C i,j W (θ) that correspond to the interactions between any two cells B i k ∈ B i and B j l ∈ B j such that the maximum Euclidean distance between any two observations x * ∈ B i k and y * ∈ B j l is less or equal to τ i,j ≥ 0 i.e. x * − y * ≤ τ i,j . Thus compute the entries (ψ ĩ k ) T C(θ)ψ j l for all ψ ĩ k ∈ D i (i.e. columns of D i ) and ψ j l ∈ D j (see Algorithm 7). In Figure 5 an example for searching local neighborhood cells of randomly placed observations in R 2 is shown. The orange nodes correspond to the source cell. By choosing a suitable value for τ the blue nodes in the immediate cell neighborhood are found by using Algorithms 4, 5 and 6.
Remark 10. Since the matrixC W (θ) is symmetrical it is only necessary to compute the blocks C i,j W (θ) for i = 1, . . . , t and j = i, . . . t.

Computation of the multi-level matrix blocks ofC W
Given an octree multi-level tree domain decomposition in R 3 , as shown in [8,9], the authors described how to apply a Kernel Independent Fast Multipole Method (KIFMM) by [32] to compute all the blocksC i,i W (θ) ∈ RÑ ×Ñ for i = −1, . . . , t in O(Ñ (t + 1) 2 ) computational steps to a fixed accuracy ε F M M > 0.
For the random projection tree it is not possible to determine a-priori the sparsity of the sparse blocksC i,j W (θ). But given a value for τ by running Algorithm 7 on every cell B i k ∈ B i , at level . . .
. . . i, with the Targetdepth corresponding for level j it is possible to determine the computational cost of constructing the sparse blocksC i,j W (θ). Suppose that maximum number of cells B j k ∈ B j given by Algorithm 7 is bounded by some γ i,j ∈ N + .
Proof. Let us look at the cost of computing all the interactions between any two cells B i k ∈ B i and B j l ∈ B j . Without loss of generality assume that i ≤ j. For the cell B l k there is at most p multi-level vectors and from Lemma 2 2 t−i+1p non zero entries. Similarly for B j l . All the interactions (ψ ĩ k ) T C(θ)ψ j l now have to be computed, where ψ ĩ k ∈ B i k and ψ j l ∈ B j l . The term C(θ)ψ j l is computed by using an FMM with 2 t−j+1p sources and 2 t−i+1p targets at a cost of O((2 t−j+1p + 2 t−i+1p ) α ). Since there is at mostp multi-level vectors in B i k and B j l then the cost for computing all the interactions (ψ ĩ k ) T C(θ)ψ j l is O(p(2 t−j+1p + 2 t−i+1p ) α + 2 t−i+1p ). Now, at any level i there is at most 2 i cells, thus the result follows.

Multi-level estimator and predictor
The multi-level random projection tree can be exploited in such a way to significantly reduce the computational burden and to further increase the numerical stability of the estimation and prediction steps. This is a direct extension of the multi-level estimator and predictor formulated in [8] for higher dimensions.

Estimator
The multi-level likelihood function, l W (θ) (see equation (12) ), has the clear advantage of being decoupled from the vector β. Furthermore, the multi-level covariance matrix C W (θ) will be more numerically stable than C(θ) thus making it easier to invert and to compute the determinant. However, it is not necessary to perform the MLE estimation on the full covariance matrix C W (θ), instead construct a series of multi-level likelihood functionsl n W (θ), n = −1, . . . t, by applying the partial transform [W T t , . . . , W T n ] to the data Z wherẽ where Z n W := [W T t , . . . , W T i ] T Z,Ñ is the length of Z n W ,C n W (θ) is theÑ ×Ñ upper-left sub-matrix ofC W (θ) and C n W (θ) is theÑ ×Ñ upper-left sub-matrix of C W (θ). A consequence of this approach is that the matrices C n W (θ), n = −1, . . . , t will be increasingly more numerically stable, thus easier to solve computationally, as shown in the following theorem.

Theorem 5. Let κ(A) → R be the condition number of the matrix
Proof. This proof is a simple extension of the proof in Theorem 3.
Remark 11. If C(θ) is symmetric positive definite then for n = −1, . . . , t the matrices C n W (θ) are symmetric positive definite. The proof is immediate.
Remark 12. If the matrixC n W (θ) is close to C n W (θ), for n = 1, . . . , d, in some matrix norm sense, the condition number ofC n W (θ) will be close to C n W (θ). In section 7, for a class of covariance functions, it can be shown that for sufficiently large τ and/or w, a ∈ N + [with the index set Λ(w) m,g orΛ(w) m,g ]C n W (θ) will be close to C n W (θ). ThusC n W (θ) will be symmetric positive definite.

Predictor
An alternative formulation for obtaining the estimateẐ(x 0 ) is by solving the following problem In [25], the authors show that the solution of this problem leads to equation (4) andγ(θ) = C −1 (θ)(Z − Mˆβ(θ)). The best unbiased predictor can be evaluated aŝ and the Mean Squared Error (MSE) at the target point x 0 is given by From (15) it is observed that M Tγ (θ) = 0. This implies thatγ ∈ R n \P p (S) and can be uniquely rewritten asγ = W T γ W for some γ W ∈ R N −p . Now, rewrite C(θ)γ + Mβ = Z as Now apply the matrix W to equation (17) and obtain W{C(θ)W Tγ W + Mβ} = WZ. Since The preconditioner P W is formed from the entries of the matrix C W . For i = −1, . . . , t and each cell B i k ∈ B i form a matrix P i,k of all the interactions (ψ ĩ k ) T C(θ)ψ ĩ k , where ψ ĩ k ∈ B i k . Thus, where ψ i 1 , . . . , ψ i s ∈ B i k and s is the total number of ML basis in the cell B i k . The preconditioner is formed as a block diagonal matrix of all the matrices P i,k for every cell B i k ∈ B, i.e. P W := diag(P −1,0 , P 0,0 , P 1,0 , P 1,1 , . . . ).
Applying the preconditioner P −1 W (θ) the following system of equations Note that in some cases C W (θ) will have very small condition numbers. For this case we can set P W := I, i.e. no preconditioner.
Theorem 6. If the covariance function φ : Γ d × Γ d → R is positive definite, then the matrix P W (θ) is always symmetric positive definite.
Proof. Immediate. Although in practice the combined NESDIS and sparse Cholesky factorization is highly efficient, as shown by our numerical results, a worse case complexity bound can be obtained only for d = 2 or d = 3 dimensions (see [8]).
Two choices for the computation of (Z n W ) TCn W (θ) −1Zn W are open to us: i) a Cholesky factorization ofC n W (θ), or ii) a Preconditioned Conjugate Gradient (PCG). The PCG choice requires significantly less memory and allows more control of the error. However, the sparse Cholesky factorization ofC n W (θ) has already been used to compute the determinant. Thus we can use the same factors to compute (Z n W ) TCn W (θ) −1Zn W .

Predictor: Computation ofC Wγ W =Z W
For the predictor stage a different approach is used. Instead of inverting the sparse matrix C W (θ) a Preconditioned Conjugate Gradient (PCG) method is employed to computeγ W = C W (θ) −1 Z W .
Recall that C W = WC(θ)W T ,γ W = Wγ and Z W = WZ. Thus the matrix vector products C W (θ)γ n W in the PCG iteration are computed within three steps: where γ 0 W is the initial guess and γ n W is the n th iteration of the PCG. (1) Transformation from multi-level representation to single level. This is done in at most O(Nt) steps.
(2) Perform matrix vector product using a fast summation method. For d = 2, 3 a KIFMM is used to compute the matrix vector products. Alternatively, for d > 3 an ASKIT method is used to compute the matrix-vector products.
The matrix-vector products C W (θ)γ n W , where γ n W ∈ R N −p , are computed in O(N α + 2Nt) computational steps to a fixed accuracy ε F M M > 0. Note that α ≥ 1 is dependent on the efficiency of the fast summation method. The total computational cost is O(kN α + 2Nt), where k is the number of iterations needed to solveC W (θ)γ W (θ) =Z W to a predetermined accuracy ε P CG > 0.
It is important to point out that the introduction of a preconditioner can degrade the performance of the PCG, in particular, if the preconditioner is ill-conditioned. The accuracy of the PCG method ε P CG has to be set such that the accuracy of the unpreconditioned system C W (θ)γ W (θ) = Z W is below a user given tolerance ε > 0.

Error analysis
Consider the full solution x W = C n W (θ) −1 Z n W and sparse solutionx W =C n W (θ) −1 Z n W for n = −1, . . . , t, then the error can be bounded as This section is divided into three parts: i) decay of the multi-level matrix C W (θ) coefficients ii) analyticity of the Matérn covariance function and finally iii) a-posteriori error estimates of

Multi-level covariance matrix decay
The decay of the coefficients of the matrix C W (θ) will depend directly on the choice of the multivariate index set Λ m,g (w) and the analytic regularity of the covariance function. In this section two lemmas are derived that provide bounds on the magnitude of the coefficients of C W (θ).
For the case that ν → ∞ the Matérn covariance function converges to a Gaussian and is analytic on Γ d × Γ d . For this case error bounds that depend on the domain Γ d × Γ d are derived in Lemma 3 with respect to the subspace Pp(S) generated by the Smolyak index set Λ m,g (w + a), where w, a ∈ N + .
In general, the Matérn covariance function will not be regular at the origin. If that is the case, care has to be taken that the sparse grid domain does not cross the origin. Error bounds that depend on the extended Smolyak index setΛ m,g (w + a), where w, a ∈ N + are derived in Lemma 4.
has an analytic extension on Eσ ,...,σ × Eσ ,...,σ and is bounded byM (φ). Let Pp(S) ⊥ be the subspace in R N generated by the Smolyak index set Λ m,g (w + a) for some w, a ∈ N + . Let S m,g w+a,d be a sparse grid generated from the index set Λ m,g (w + a) with η(w + a) collocation nodes. For i, j = 0, . . . , t consider any multi-level vector ψ i m ∈ Pp(S) ⊥ , with n m non-zero entries, from the cell B i m ∈ B i and any multi-level vector ψ j q ∈ Pp(S) ⊥ , with n q non-zero entries, from the cell B j q ∈ B j , then: is the identity operator, and ǫ 2 (x, y; θ) := (I d ⊗ I d − I d ⊗ S m,g w+a,d )ǫ 1 (x, y; θ). Consider any multilevel vector ψ i m ∈ Pp(S) ⊥ in the cell B m and any vector ψ j q ∈ Pp(S) ⊥ from the cell B q , where m, q = 0, . . . , t, then The last equality follows from the fact that from Proposition 1 ( The last equality is due to the fact that ψ i m and ψ j q are orthonormal. The next step is to bound the error ǫ 2 (x, y; θ) L ∞ (Γ d ×Γ d ) . Now, extend y → z ∈ Eσ ,...,σ . For this case ǫ 1 (·, z; θ) becomes a complex number. If w > d/ log (2) then from Theorem 2 and equation (11) for all z ∈ Eσ ,...,σ whereM (φ) = max w,z∈Eσ ,...,σ ×Eσ ,...,σ |φ(w, z)|. From Theorem 7 and equation (10) it follows that where Now, if w ≤ d/ log (2) then from Theorem 2 for all z ∈ Eσ ,...,σ

From equation (11)
There are many cases where it is not possible, or hard to prove that the covariance function φ(x, y; θ) has a bounded analytic extension on a well defined region in C d × C d . For example, for the Matérn covariance function (ν = 0.5), there is a derivative discontinuity at the origin. We look at a class of covariance functions where it is possible to prove decay rates of the covariance matrix C W (θ) whenever any two cells B i m ∈ B i and B j q ∈ B j are well separated. Let f be a R d valued vector, where for n = 1, . . . , d, f n := θ n (x n − y n ) 2 . The idea is to form a sparse grid representation of the covariance function φ(f) : R d → R that avoids the analytic singularity at the origin.
Consider any two cells B i m ∈ B i and B j q ∈ B j . Assume that the minimum Euclidean distance between any two observations x * ∈ B i m and y * ∈ B j q is τ i,j > 0 i.e. x * − y * ≥ τ i,j . Conversely the maximal distance between all the observations in B i m and B j q is x * − y * ≤ 2 √ d. Let Let t be a R d valued vector and for n = 1, . . . , d let Thus t ∈ Γ d and the covariance function is recast as φ(t; θ) : Γ d → R. Now, for this case the multi-level basis constructed from the observations S and the monomial setQ a Λ m,g (w) , for some a, w ∈ N + , is not enough to show decay of the coefficients of C W (θ). Instead, assume that the extended multi-level basis (see Remark 9) is constructed from the set of monomialsQ ã Λ m,g (w) , with extended Smolyak index setΛ m,g (w + a) Lemma 4. Suppose that φ(t; θ) ∈ C 0 (Γ d ; R) has an analytic extension on Eσ ,...,σ and is bounded byM (φ). Let Pp(S) ⊥ be the subspace in R N generated by the extended Smolyak index setΛ m,g (w+ a) for some w, a ∈ N + . Let S m,g w+a,d be a sparse grid generated from the index setΛ m,g (w + a) with η(w + a) collocation nodes. For m, q = 0, . . . , t consider any extended multi-level vector ψ i m ∈ Pp(S) ⊥ , with n m non-zero entries, from the cell B i m and any extended multi-level vector ψ j q ∈ Pp(S) ⊥ , with n q non-zero entries, from the cell B j q .

Analyticity of the Matérn covariance function
In this section the analytic extension for the Matérn covariance function are positive constants, diag(θ) ∈ R d×d is a zero matrix with the vector θ on the diagonal, and x, y ∈ Γ d . Recall from Remark 6 that for different values of ν ∈ R + different shapes of the covariance function are obtained, i.e. for ν = 1/2 + n, n ∈ N + φ(r; ν, θ) = exp(− √ 2νr(θ)) Γ(n + 1) The Gaussian covariance function is first analyzed, which is the case for the Matérn covariance function when ν → ∞. However, the covariance function is recast in a slightly different form: φ(x, y; θ) = exp −(x − y) T diag(θ)(x − y) .

Lemma 5.
There exists an analytic extension of the Gaussian correlation function φ(x, y; θ) : Proof. First, the following facts from complex analysis are stated: i) The exponential function is holomorphic on C.
ii) The polynomial functions are holomorphic on C.
iii) The sum of two holomorphic functions on C are holomorphic on C.
iv) The composition of two holomorphic functions on C are holomorphic on C.
From ii) and iii) for all x, y ∈ C d it follows that (x − y) T diag(θ)(x − y) : C d × C d → C is holomorphic. From i) and iv) it is concluded that φ(x, y; θ) : The next step is to bound the function φ(x, y; θ) : C d × C d → C. However, for the isotropic sparse grid representation used in this paper it is sufficient to show that the analytic extension of φ(x, y; θ) is valid on d n=1 E n,σ × d n=1 E n,σ . Now, recall the shape of the polyellipse E n,σn = z ∈ C,σ ≥ δ n ≥ 0; Re z = e δn + e −δn 2 cos(θ), Im z = e δn − e −δn 2 sin(θ), θ ∈ [0, 2π) .
Therefore any z ∈ E n,σ can be written in the form (e 2δn +e −2δn ) 1/2 2 e iΩn , i = √ −1,σ ≥ δ n ≥ 0 and Ω n ∈ R. Thus for all x, y ∈ d n=1 E n,σ Now, the analyticity of the Matérn covariance function is analyzed. Since, in general, the Matérn covariance function is not differentiable at the origin then there is no analytic extension on d n=1 E n,σ × d n=1 E n,σ , thus Lemma 4 will not be valid. However, it is not necessary to show that an analytic extension of φ(x, y; θ) : Γ d × Γ d → R exists for the entire domain Γ d × Γ d , but only on a sub-domain that does not cross the origin.
Consider any two cells B i m ∈ B i and B j q ∈ B j , at levels i and j, and assume that the minimum Euclidean distance between any two observations x * ∈ B i m and y * ∈ B j q is τ i,j > 0 i.e. x * − y * ≥ τ i,j . Conversely the maximal distance between all observations in B i m and B j q is x * − y * ≤ 2 √ d. Let {x min , y min } be the minimal distance pair observations and similarly, let {x max , y max } be the maximal. Let t ∈ R d , where for n = 1, . . . , d, t n is equal to equation (19), and recast the Matérn covariance function as φ(t; ν, θ) : t ∈ Γ d and r(t, θ) = d n=1 θ n ((t n + 1)(M 2 n − m 2 n )/2 + m 2 n ). The next step is to find a suitable value forσ > 0 such that the Matérn function is analytic and bounded on the region d n=1 E n,σ . In the following Lemma an estimate of the region of analyticity is obtained. Lemma 6. Given any two cells B i m ∈ B i and B j q ∈ B j , such that the minimum Euclidean distance between any two observations x * ∈ B i m and y * ∈ B j q is τ i,j > 0, then there exists an analytic extension of φ(t; ν, θ) : Moreover, for all t ∈ d n=1 E n,σ ⊂ C d it follows that |r(t; θ)| ≤ α(θ,σ i,j ) := θ 2 2(3 + √ 2eσ i,j ) 1 2 and for n ∈ N + |φ(r; 1/2 + n, θ)| ≤ exp(−(4(3 + √ 2)) 1/2 ν θ 2 ) Γ(n + 1) Γ(2n + 1) Proof. From equation (18) for n = 1, . . . , d thus, d n=1 m 2 n ≥ τ i,j and for n = 1, . . . , d, m 2 n ≤ M 2 n ≤ 4. Recall that the Matérn function can be rewritten as where r(f, θ) = This condition is satisfied if for n = 1, . . . , d, Re t n > −1 − 2d , ∞). The next step is to embed a polyellipse E log ρ ⊂ (−1 − τ i,j 2d , ∞) × (−∞i, ∞i), as shown in Figure 6. This is achieved by setting ρ < 2d . The embedding is achieved for a polyellipse with ρ < The result follows.

A-posteriori error estimates
Given that the region of analyticity of the Matérn covariance function and decay rates of the covariance matrix haven been derived, estimates for the error between the full x W = C n W (θ) −1 Z n W and sparse solutionx W =C n W (θ) −1 Z n W are shown in this section. For n = t, . . . , −1 suppose that σ min (C n W (θ)) Note that the following singular value inequalities, which can be found in [18], where applied: If A, B ∈ C n,n then Remark 13. The terms σ min (C n W (θ)) and σ max (C n W (θ)) can be estimated by using a Lanzcos sparse Singular Value Decomposition (SVD) type algorithm [5,20].

Remark 14.
It is not too hard to see that for τ i,j ≥ 0 Algorithm 5 leads to E i,j = 0 if i = −1, 0 or j = −1, 0.
The contributions of the errors E i,j 2 , for i, j = max {n, 1}, . . . , t, can now be estimated, but first, the following assumption is made: Assumption 2. For i, j, = 1, . . . , t let τ i,j := τ 2 t−(i+j)/2 , for some user defined parameter τ > 0. Now, for i, j = max {1, n}, . . . , t if w + a > d/ log (2) then from Assumption 2 and Lemma 1, 2 and 3 and thus Following a similar argument for the case that w + a ≤ d/ log (2) the following theorem is proved: Suppose the same conditions as in Lemma 3 are satisfied, then from Assumption 2 where: (a) If w + a > d/ log (2) then Remark 15. By combining Theorem 7 and Lemma 5 the a-posteriori error bound is computed for the Gaussian covariance function. For this bound the user can set the size of the polyellipse parameterσ to one, thus σ = 1/2,M (φ) = d n=1 exp(θ n ) 2d(e 2σ +e −2σ )) . Now one must deal with the more general Matérn covariance function error estimates . It is straightforward to prove the following theorem.

Numerical results
The performance of the multi-level solver for estimation and prediction formed from random datasets in high dimensions is tested. The results show that the computational burden is significantly reduced. In particular, it is possible to now solve ill-conditioned problems efficiently. The implementation of the code is done with a combination of MATLAB and C++ packages: i) Matlab: The RP/kD binary tree, multi-level basis construction, formation of the sparse matrix C W , estimation and prediction components are written and executed on Matlab [23].
The computational bottlenecks are executed on the C++ software packages, dynamic shared libraries, and the highly optimized BLAS and LAPACK packages contained in MATLAB.
ii) Direct and fast summation: The matlab code estimates the computational cost between the direct and fast summation methods and chooses the most efficient approach. For the direct method the highly optimized MATLAB arithmetic functions are used. For the fast summation method the KIFMM (d = 3) or ASKIT (d > 3) c++ codes are used. The KIFMM is modified to include a Hermite interpolant approximation of the Matérn covariance function, which is implemented with the intel MKL package [1] (see [8] for details). The ASKIT code is only used for the Gaussian case since the performance for the Matérn kernel is only slightly better than the direct approach.
iii) Dynamic shared libraries: These are produced with the GNU gcc/g++ packages. These libraries implement the Hermite interpolant with the intel MKL package (about 10 times faster than Matlab Matérn interpolant) and link the MATLAB code to the KIFMM. iv) Cholesky and determinant computation: The Suite Sparse 4.2.1 package ( [11,[13][14][15][16]) is used for the determinant computation of the sparse matrixC W (θ).
The code is executed on a single CPU (4 core Intel i7-3770 CPU @ 3.40GHz.) with Linux Ubuntu 14.04 and 32 GB memory. To test the effectiveness of the Multi-Level solver the following data sets are generated: a) Random hypercube data set: The set of random observation locations C d 0 , . . . , C d 10 vary from 1,000, 2000, 4,000 to 512,000 knots generated on the hypercube [−1, 1] d for d dimensions.
The observations locations are also nested, i.e. C d 0 ⊂ · · · ⊂ C d 10 . b) Random n-sphere data set: The set of nested random observations S d 0 ⊂ · · · ⊂ S d 10 varies from 1,000, 2000, 4000 to 128,000 knots generated on the n-sphere S d−1 : c) Gaussian data set The set of observations values Z d 0 , Z d 1 , . . . Z d 5 are formed from the Gaussian random field model (1) for 1,000, 2,000, . . . 32, 000 observation locations. The data set Z d n is generated from the set of nodes C d n or S d n , with the covariance parameters (ν, ρ) and the corresponding set of monomials Q Λ m,g (w) . Due to memory limitations 32,000 was the maximum size of the generated data set. d) All the numerical test are done assuming thep = p.

Numerical stability and sparsity of the covariance multi-level matrix C W
The covariance matrix C(θ) becomes increasingly ill-conditioned as the number of observations are increased, leading to instability of the numerical solver. It is now shown how effective Theorem 3 becomes in practice. In Figure 7 the condition number of the multi-level covariance matrix C W (θ) is plotted with respect to the cardinality p of the index set Λ(w) for different w levels. The multilevel covariance matrix C W (θ) is built from the random cube C d 4 or n-sphere S d 4 observations. The covariance function is set to Matérn with ν = 1 and ρ = 1, 10, 100. As the plots confirm the covariance matrix condition number significantly improves with increasing level w for the TD, SM and HC index sets. This is in contrast with the large condition numbers of the original covariance matrix C(θ).
In Figure 8(a) and (b) the log determinant relative error | log detC W (θ)−log det C W (θ)|| | log det C W (θ)| is plotted with respect to the sparsity ofC W (θ) by choosing an increasingly larger distance criterion parameter τ . As it is observed even for highly sparse multi-level covariance matrices the relative error is about 1% and further improves as the sparsity is increased. Note that sparsity equal to one means that the matrix is fully populated.
The first observation to notice is that all the sparse matricesC i W (θ), i = t, t − 1, . . . are very well conditioned, thus numerically stable. This is in contrast to the original covariance matrices that are poorly conditioned, as are shown in Table 6. The sparsity ofC i W (θ) and the Cholesky factor G are shown in columns 8 and 10. The construction time t con is shown in column 9.
It is noted that the sparse matrices in Table 3 are built with a direct summation method. Although it appears that the multi-level covariance matrices are constructed close to linear time, the computational cost should approach quadratic time for a large enough number of observation nodes and 0 ≤ i << t. This computational cost will be further improved as better high dimensional fast summation techniques are developed. Table 3: Sparsity test on the sparse matricesC i W , i = t, t − 1, . . . . The polynomial space of the index set Λ(w) is restricted to the Hyperbolic Cross on a n-Sphere with d = 50 dimensions. The domain decomposition is formed from a RP tree. The level of the index set is w = 5, which corresponds p = 1426. The kernel function is Matérn with ν = 1.25, ρ = 10. The first column is the number of random n-Sphere nodes. The second is the maximum level of the RP tree and i is the level of the sparse matrixC i W . The forth column is the size of the matrixC i W and τ is the distance criterion parameter. The sixth column, t M L , is the total time for the construction of the multi-level basis. The seventh column is the sparsity ofC i W . The eighth column, t con is the total time for the construction of the matrixC i W . The ninth column is the sparsity of the Cholesky factor G (with nested dissection reordering) of the sparse matrixC i W . The last column is the total time to compute the Cholesky factor G. Total degree log det relative error   Table 4: Estimation of parametersν andρ with i) Total Degree polynomial index set Λ(w), ii) (a) kD or (b) RP tree, iii) n-Sphere with (a) d = 10 or (b) d = 50 dimensions, iv) The parameters of the observation data Z are formed with the following Matérn parameters: ν = 1.25 and ρ = 1, v) The number of realizations of the Gaussian random field model is set to M = 100. and vi) The distance criterion τ is set to 10 −7 . The first to forth columns are self-explanatory. The fifth column is the mean error ofν with M realization. The sixth column is the mean error ofρ. The last two columns are the standard deviation of M realizations of the parametersν andρ.

Prediction
The multi-level kriging prediction approach is now tested on the numerical solution of the system of equationsC W (θ)γ W =Z W for d = 3 and d = 50 dimensions. Pre-conditioned results for computingC W (θ)γ W =Z W for the hypercube data set with d = 3 dimensions, kD tree, and the Total Degree index set Λ(w) are shown in Table 5. The Matérn covariance coefficients θ = (ν, ρ) are set to (3/4,1/6). The relative error of the residual of PCG method for the unpreconditioned system is set to ε = 10 −3 . The KIFMM is set to high accuracy.
The construction of the preconditioner P W grows initially with quadratic time. However, around 128,000 observations the KIFMM starts to kick in and the rate is reduced to almost linear time, as expected. The condition number of the covariance matrices are fairly large, making this problem somewhat hard to solve numerically. The results show that 512,000 size problems with good accuracy are feasible with a single 4-core processor. Table 5: Pre-conditioned results for computingC W (θ)γ W =Z W for the hypercube data set with d = 3 and the Total Degree index set Λ(w). The Matérn covariance coefficients θ = (ν, ρ) are set to (3/4,1/6). The relative error of the residual of PCG method for the unpreconditioned system is set to ε = 10 −3 . The KIFMM is set to high accuracy. The second column is the condition number of the covariance matrix C. The third column is the number of iterations needed to obtain 10 −3 relative error of the unpreconditioned system with C W (θ). Let itr(C W ) be the number of CG iterations needed for convergence for ǫ = 10 −3 residual accuracy. The fourth column presents the wall clock times for the preconditioner computation. The PCG iteration wall clock times for C W are given in the fifth column. The last column represents the total wall clock time to compute γ W =C W (θ) −1Z W . In Table 6 the multilevel prediction method is tested on a Gaussian exp(− r 2 2h 2 ) and Matérn kernels for d = 50 dimensions. The first test is for the Gaussian covariance with h = √ 10 using the ASKIT for the matrix vector products. This case gives rise to badly condition covariance matrices that are difficult to solve with an iterative method as shown from the very large condition numbers of the second column. However, the multi-level prediction method is feasible to solve large problems in high dimensions. Note the wall clock time for convergence grows quadratically. This is mostly due to the quadratic performance of ASKIT as the number of observations increases above 32,000. However, it is still better than using a direct approach.
The second test is for the Matérn covariance function with ν = 5/4 and ρ = 1. The direct method is used for the matrix vector products as the performance of ASKIT is only just slightly better. The covariance matrices C(θ) still have large condition numbers. However, convergence of the multi-level covariance matrix requires only a few iterations due to the excellent condition numbers. For this case it is not necessary to compute the pre-conditioner.
i) The multilevel method is numerically stable. Hard estimation and prediction of large dimensional problems are now feasible.
ii) The method is efficiently implemented by using a combination of MATLAB, c++ software packages and dynamic libraries.
iii) Numerical results of up to 50 dimensional problems for estimation and prediction. These problems are difficult to solve directly due to the large condition numbers, but feasible with the multi-level method.
iv) A-posteriori error estimates for the solution of C W (θ)γ W = Z W are derived for extended Smolyak index sets.
v) The efficiency of this approach will be further improved as high dimensional fast summation methods are developed.
vi) A-posteriori estimates for Hyperbolic Cross index sets will be studied in a future paper. This requires deriving L ∞ (Γ d ) error estimates for the HC sparse grid.