Sparse residual tree and forest

TL;DR

This paper introduces Sparse Residual Trees and Forests, adaptive methods for multivariate data approximation that are efficient, stable, and capable of identifying regions needing data refinement, with proven convergence and promising numerical results.

Contribution

It proposes a novel hierarchical parallel SRT algorithm combined with RBF explorations, and introduces SRF to enhance approximation accuracy and stability.

Findings

Achieves $ ext{O}(N ext{log}_2N)$ time complexity for initial work

Provides convergence results for SRTs and SRFs

Demonstrates effectiveness through numerical experiments

Abstract

Sparse residual tree (SRT) is an adaptive exploration method for multivariate scattered data approximation. It leads to sparse and stable approximations in areas where the data is sufficient or redundant, and points out the possible local regions where data refinement is needed. Sparse residual forest (SRF) is a combination of SRT predictors to further improve the approximation accuracy and stability according to the error characteristics of SRTs. The hierarchical parallel SRT algorithm is based on both tree decomposition and adaptive radial basis function (RBF) explorations, whereby for each child a sparse and proper RBF refinement is added to the approximation by minimizing the norm of the residual inherited from its parent. The convergence results are established for both SRTs and SRFs. The worst case time complexity of SRTs is $\mathcal{O}(N\log_2N)$ for the initial work and $\mathcal{O}(\log_2N)$ for each prediction, meanwhile, the worst case storage requirement is $\mathcal{O}(N\log_2N)$, where the $N$ data points can be arbitrary distributed. Numerical experiments are performed for several illustrative examples.

Equations

$$h:=h_{X,\Omega }:=\underset{x\in \Omega }{\sup }\underset{x_i\in X}{\min }\Vert x-x_i{\Vert }_2,$$

$$s_f(x,\alpha )=\sum _{j=1}^{N^{\prime \prime }}{\alpha }_{j}G(x-x_j^{\prime \prime }),x\in {\Omega }^{\prime }$$

$$\min_{\alpha\in\mathbb{R}^{N^{\prime\prime}}}\sum_{i=1}^{N^{\prime}}\Big{(}s_{f}(x^{\prime}_{i},\alpha)-f(x^{\prime}_{i})\Big{)}^{2}.$$

$$\underset{\alpha \in {\mathbb{R}}^{N^{\prime \prime }}}{\min }\Vert {\Phi }_{X^{\prime },X^{\prime \prime }}\alpha -f_{X^{\prime }}{\Vert }_2^2,$$

$$|f{|}_{W_p^k(\Omega )}={\left(\sum _{|\gamma |=k}\Vert D^{\gamma }f{\Vert }_{L_p(\Omega )}^p\right)}^{1/p}\text{and}\Vert f{\Vert }_{W_p^k(\Omega )}={\left(\sum _{|\gamma |⩽k}\Vert D^{\gamma }f{\Vert }_{L_p(\Omega )}^p\right)}^{1/p}.$$

$${\mathcal{N}}_G({\mathbb{R}}^d)=\left\{f\in C({\mathbb{R}}^d)\cap L_2({\mathbb{R}}^d):\Vert f{\Vert }_{{\mathcal{N}}_G}:={\left({\int }_{{\mathbb{R}}^d}|\hat{f}(\omega ){|}^2{e}^{\frac{\Vert \omega {\Vert }_2^2}{4{\delta }^2}}\mathrm{d}\omega \right)}^{\frac{1}{2}}<\infty \right\},$$

$${\mathcal{N}}_G(\Omega )=\left\{f{|}_{\Omega }:f\in {\mathcal{N}}_G({\mathbb{R}}^d)\text{and}(f,g{)}_{{\mathcal{N}}_G({\mathbb{R}}^d)}=0,\forall g\in {\mathcal{N}}_G({\mathbb{R}}^d)\text{s.t.}g{|}_{\Omega }=0\right\},$$

$${\mathcal{N}}_G(\Omega )\subset W_2^k(\Omega )\text{with}\Vert f{\Vert }_{W_2^k(\Omega )}⩽C_G^k{k}^{k/2}\Vert f{\Vert }_{{\mathcal{N}}_G},\forall f\in {\mathcal{N}}_G(\Omega ),$$

$$(f,g{)}_{{\mathcal{N}}_M({\mathbb{R}}^d)}={\int }_{{\mathbb{R}}^d}\hat{f}(\omega )\overline{\hat{g}(\omega )}{\hat{M}}^{-1}(\omega )\mathrm{d}\omega ,\forall f,g\in {\mathcal{N}}_M({\mathbb{R}}^d),$$

$${\mathcal{N}}_M(\Omega )\subset W_2^k(\Omega )\text{with}\Vert f{\Vert }_{W_2^k(\Omega )}⩽C_M^k{k}^k\Vert f{\Vert }_{{\mathcal{N}}_M},\forall f\in {\mathcal{N}}_M(\Omega ),$$

$$\mathrm{RAE}=\frac{{\max }_i|s(x_i)-f(x_i)|}{{\max }_i|f(x_i)|}.$$

$${\mathcal{P}}_{X^{\prime }(I)}(x)=\frac{{\mathcal{P}}_{X^{\prime },I}(x,i)}{{\mathcal{P}}_I(i)}=\frac{{\mathcal{P}}_{X^{\prime }}(x){\mathcal{P}}_I(i)}{{\mathcal{P}}_I(i)}={\mathcal{P}}_{X^{\prime }}(x),x\in {\Omega }^{\prime },i\in I;$$

$$u_1^{\prime }=\arg \underset{x^{\prime }\in X^{\prime }(I)}{\min }\Vert x^{\prime }-\overline{X^{\prime }(I)}{\Vert }_2,\text{where}\overline{X^{\prime }(I)}=\frac{1}{N_I^{\prime }}\sum _{i=1}^{N_I^{\prime }}{x}_{I(i)}^{\prime }.$$

$$u_{j+1}^{\prime }=\arg \underset{x^{\prime }\in X^{\prime }(I)}{\max }\left(\underset{1⩽l⩽j}{\min }\Vert x^{\prime }-u_l^{\prime }{\Vert }_2\right),$$

$$\omega =({\omega }_1,{\omega }_2,{\omega }_3,{\omega }_4{)}^{\mathrm{T}},$$

$${\delta }^{\prime }=\sqrt{-\frac{\ln ({\omega }_3)}{{\max }_{x^{\prime }\in X^{\prime }(I)}\Vert x^{\prime }-\overline{X^{\prime }(I)}{\Vert }_2^2}},\text{where}\overline{X^{\prime }(I)}=\frac{1}{N_I^{\prime }}\sum _{i=1}^{N_I^{\prime }}{x}_{I(i)}^{\prime };$$

$$Q_v{R}_v={\Phi }_{X^{\prime }(I),\{{\chi }_1^{\prime },\cdots ,{\chi }_v^{\prime }\}}\in {\mathbb{R}}^{N_I^{\prime },v},1⩽v⩽j,$$

$$r_j^{\prime }(X^{\prime }(I))=r^{\prime }(X^{\prime }(I))-\sum _{1⩽i⩽j}{\alpha }_j^{(i)}{e}^{-{\delta }^{\prime 2}\Vert X^{\prime }(I)-{\chi }_i^{\prime }{\Vert }_2^2}.$$

$${\chi }_{j+1}^{\prime }=u_{l^{\ast }}^{\prime }\in U_{j+d+1}^{\prime }-\{{\chi }_1^{\prime },\cdots ,{\chi }_j^{\prime }\},$$

$$l^{\ast }=\arg \underset{l\in \Gamma }{\max }\sum _{x^{\prime }\in {\Lambda }_l\cap X^{\prime }(I)}\frac{|r_j^{\prime }(x^{\prime }){|}^2}{n_l},$$

$$\kappa (R)>{\omega }_1\text{or}{ϵ}_j-{ϵ}_{j+1}<{\omega }_2\text{or}j+1=N_I^{\prime },$$

$${ϵ}_j=\sqrt{\frac{1}{N_I^{\prime }}\sum _{x^{\prime }\in X^{\prime }(I)}(r_j^{\prime }(x^{\prime }){)}^2}.$$

$$r^{\prime \prime }(X^{\prime })=r^{\prime }(X^{\prime })-\sum _{1⩽i⩽N^{\prime \prime }}{\alpha }_i^{\prime }{e}^{-{\delta }^{\prime 2}\Vert X^{\prime }-{\chi }_i^{\prime }{\Vert }_2^2}.$$

$${\bar{n}}_c=\frac{1}{M}\sum _{i=1}^M{n}_c^{(i)},$$

$$u_1^{\prime }=\arg \underset{x^{\prime }\in X^{\prime }(I)}{\max }\Vert x^{\prime }-\overline{X^{\prime }(I)}{\Vert }_2,\text{where}\overline{X^{\prime }(I)}=\frac{1}{N_I^{\prime }}\sum _{i=1}^{N_I^{\prime }}{x}_{I(i)}^{\prime }.$$

$$x_a^{\prime }=u_{l^{\ast }}^{\prime },$$

$$x_b^{\prime }=\arg \underset{x^{\prime }\in X^{\prime }(I)}{\max }\Vert x^{\prime }-x_a^{\prime }{\Vert }_2.$$

$$P_{{\vec{n}}^{\prime }}(X^{\prime })=X^{\prime }{\vec{n}}^{\prime },\text{where}{X}^{\prime }\in {\mathbb{R}}^{N^{\prime }\times d}\text{and}{\vec{n}}^{\prime }\in {\mathbb{R}}^{d\times 1};$$

$$X_1^{\prime }=\{x^{\prime }\in X^{\prime }:P_{{\vec{n}}^{\prime }}(x^{\prime })⩽c^{\prime }\}\text{and}{X}_2^{\prime }=X^{\prime }-X_1^{\prime };$$

$${\Omega }_1^{\prime }=\{x^{\prime }\in {\Omega }^{\prime }:x^{\prime }{\vec{n}}^{\prime }⩽c^{\prime }\}\text{and}{\Omega }_2^{\prime }={\Omega }^{\prime }-{\Omega }_1^{\prime }.$$

Full text

††journal:

1 Introduction

Multivariate scattered data approximation problems arise in many areas of engineering and scientific computing. In the last five decades, radial basis function (RBF) methods have gradually become an extremely powerful tool for scattered data. This is not only because they possess the dimensional independence and remarkable convergence properties (see, e.g., [1, 2, 3, 4, 5]), but also because a number of techniques, such as multipole (far-field) expansions [6, 2], multilevel methods of compactly supported kernels [7, 8, 2, 9] and partition of unity methods [10, 11, 2], have been proposed to reduce both the condition number of the resulting interpolation matrix and the complexity of calculating the interpolant. These techniques are, of course, very important in practice, however, in contrast to the stability and efficiency, maybe the later question is the most crucial one for a general representation of functions, that is, how to accurately capture and represent the intrinsic structures of a target function, especially in high dimensional space.

More specifically, when the data and the expected accuracy are given, we usually do not know at all whether the data is redundant or insufficient for the target function. So it is necessary to consider the following three questions:

Whether the current data is just right to reach the expected accuracy?

• 2.

How to establish a sparse approximation by ignoring the possible redundancy?

• 3.

How to update the approximation by replenishing the possible insufficiency?

It is often difficult to distinguish between data insufficiency and redundancy, and they could in fact exist simultaneously in different local regions.

Sparse residual tree (SRT) is developed for the purpose of representing the intrinsic structure of arbitrary dimensional scattered data. SRT is based on both tree decomposition and adaptive radial basis function (RBF) explorations. For each child a concise and proper RBF refinement, whose shape parameter is related to the current regional scale, is added to the approximation by minimizing the $2$-norm of the residual inherited from its parent; then the tree node will be further split into two according to the updated residual; and this process finally stops when the data is insufficient or the expected accuracy is reached.

The word “sparse” here has two meanings: (i) the RBF exploration applies only to a sparse but sufficient subset of the current data, which is to ensure the efficiency of the training process; (ii) the centers of the RBF refinement are also sparse relative to the sparse subset, which is to ensure the efficiency of the prediction process. Thus, on the one hand, SRT provides sparse approximations in areas where the data is sufficient or redundant, and on the other hand, SRT points out the possible local regions where data refinement is needed. In order to ensure stability, the condition number is strictly controlled for every refinement. Furthermore, SRT also yields the excellent performance in terms of efficiency. Similar to most typical tree-based algorithms [12, 13], the worst case time complexity of SRTs is $\mathcal{O}(N\log_{2}N)$ for the initial training work and $\mathcal{O}(\log_{2}N)$ for each prediction; and the worst case storage requirement is $\mathcal{O}(N\log_{2}N)$, where the $N$ data points can be arbitrary distributed. The training process can be accelerated using multi-core architectures. This hierarchical parallel algorithm allows one to easily handle ten millions of data points on a personal computer, or much more on a computer cluster.

Although there are some different attempts to combine tree structures and RBF methods in the field of machine learning (see, e.g., [14, 15, 16]), they have not paid any attention to their convergence. In fact, similar to multilevel methods [2], these combinations do not always guarantee convergence. Most of the previously used error estimates for RBF interpolation depend on the so-called power function [17, 18, 2, 4]. But recently, sampling inequalities have become a more powerful tool in this respect, and not limited to the case of interpolation [19, 20, 1, 3]. Sampling inequalities describe the fact that a differentiable function whose derivatives are bounded cannot attain large values if it is small on a sufficiently dense discrete set. Together with the stability of the least squares framework for residual trees, we prove that a SRT based on arbitrary basis functions leads to algebraic convergence orders for finitely smooth functions. Further combining the appropriate embeddings of certain native spaces, we also prove that the Gaussian or inverse multiquadric based SRT leads to exponential convergence orders for infinitely smooth functions.

Since the SRT approximation is actually piecewise smooth, the error of each piece is significantly larger near the boundary. And the sparse residual forest (SRF), which is a combination of SRT predictors with different tree decompositions, is specifically designed to improve this situation. For all SRTs in the SRF, the splitting method of each SRT depends on the values of a random vector sampled independently and with the same distribution. This provides an opportunity to avoid those predictions with large squared deviations and to use the average value of the remaining predictions to enhance both stability and convergence. In practice, SRFs composed of a small number of SRTs perform quite well than individual SRTs; and in theory, similar to random forests [21], the error for SRFs converges with probability $1$ to a limit as the number of SRTs in the SRF becomes large. It is more efficient and accurate than the traditional partition of unity method for overcoming the boundary effect of the error.

The remainder of the paper is organized as follows. After appropriate notation and preliminaries are introduced in section 2, section 3 and section 4 give the frameworks of the SRTs and SRFs, respectively, and the stability, convergence and complexity of the SRT algorithm are discussed in section 5. A series of numerical experiments is given in section 6. In section 7, we draw some conclusions on the new method presented in this work and discuss possible extensions.

2 Notation and Preliminaries

Throughout the paper, $e$ denotes Euler’s constant, the space dimension $d\in\mathbb{N}$, the domain $\Omega\subset\mathbb{R}^{d}$ is convex, $f:\Omega\to\mathbb{R}~{}\textrm{or}~{}\mathbb{C}$ is a given target function, $X=\{x_{i}\}_{i=1}^{N}\in\Omega$ is a set of pairwise distinct interpolation points with the fill distance

[TABLE]

and $f_{X}=(f(x_{1}),\cdots,f(x_{N}))^{\mathrm{T}}$ are known function values.

Remark  2.1**.**

It is worth noting that $\Omega$ can also be extended to a finite union of convex domains, thereby $\Omega$ is bounded with Lipschitz boundary and satisfies an interior cone condition. In this case, we can first deal with these convex domains separately and then combine them into a meaningful whole by a suitable partition of unity, see section 6 for examples.

We will focus mainly on the Gaussian kernel $G(x)=G_{\delta}(x):=e^{-\delta^{2}\|x\|_{2}^{2}}$, where $\delta>0$ is often called the shape parameter. Suppose that $\Omega^{\prime}\subseteq\Omega$ is also convex, then for the subset $X^{\prime}=\{x^{\prime}_{i}\}_{i=1}^{N^{\prime}}=X\cap\Omega^{\prime}$ and selected centers $X^{\prime\prime}=\{x^{\prime\prime}_{i}\}_{i=1}^{N^{\prime\prime}}\subseteq X^{\prime}$, where $N^{\prime\prime}\leqslant N^{\prime}\leqslant N$, an Gaussian RBF approximation $s$ is required to be of the form

[TABLE]

with unknown coefficients $\alpha=(\alpha_{1},\cdots,\alpha_{N^{\prime\prime}})^{\mathrm{T}}$. Consider the following least squares (LS) problem

[TABLE]

It is worth noting here that we consider the case of $N^{\prime\prime}\ll N^{\prime}$ as a sparse approximation, and (2) can be rewritten in matrix form as

[TABLE]

where the matrix $\Phi_{X^{\prime},X^{\prime\prime}}\in\mathbb{R}^{N^{\prime}\times N^{\prime\prime}}$ is generated by the Gaussian kernel $G(x)$. Suppose $\Phi_{X^{\prime},X^{\prime\prime}}$ have a $QR$ decomposition $\Phi_{X^{\prime},X^{\prime\prime}}=QR$, where $Q\in\mathbb{R}^{N^{\prime}\times N^{\prime\prime}}$ has orthonormal columns and $R\in\mathbb{R}^{N^{\prime\prime}\times N^{\prime\prime}}$ is upper triangular, then the problem (2) has a unique solution $\alpha^{*}=\Phi_{X^{\prime},X^{\prime\prime}}^{-1}f_{X^{\prime}}=R^{-1}Q^{\mathrm{T}}f_{X^{\prime}}$, where $R$ and $Q^{\mathrm{T}}f_{X^{\prime}}\in\mathbb{R}^{N^{\prime\prime}}$ can be recursively obtained without computing $Q$ by Householder transformations [22].

We shall consider functions from certain Sobolev spaces $W_{p}^{k}(\Omega)$ with $1\leqslant p<\infty$ and native spaces of Gaussians, $\mathcal{N}_{G}(\Omega)$, respectively. The Sobolev space $W_{p}^{k}(\Omega)$ consists of all functions $f$ with distributional derivatives $D^{\gamma}f\in L_{p}(\Omega)$ for all $|\gamma|\leqslant k$, $\gamma\in\mathbb{N}_{0}^{d}$. Associated with these spaces are the (semi-)norms

[TABLE]

For the Gaussian kernel $G(x)=e^{-\delta^{2}\|x\|_{2}^{2}}$ the native space on $\mathbb{R}^{d}$ is given by

[TABLE]

further, the native space $\mathcal{N}_{G}(\Omega)$ on a bounded domain $\Omega$ is defined as

[TABLE]

where $(f,g)_{\mathcal{N}_{G}(\mathbb{R}^{d})}=\int_{\mathbb{R}^{d}}\hat{f}(\omega)\overline{\hat{g}(\omega)}e^{\frac{\|\omega\|_{2}^{2}}{4\delta^{2}}}\mathrm{d}\omega$. For any $\Omega\subseteq\mathbb{R}^{d}$ and all $k\geqslant 0$,

[TABLE]

where $C_{G}=\sqrt{\max(\delta^{-d},1)(\frac{8\delta^{2}}{e}+2)}$ depends only on the shape parameter $\delta$ and the space dimension $d$, see Theorem $7.5$ of [1] for details.

We can also consider inverse multiquadrics $M(x)=M_{\delta}(x)=(1/\delta^{2}+\|x\|_{2})^{-\beta}$ for $\beta>\frac{d}{2}$, and the inner product of native spaces $\mathcal{N}_{M}(\mathbb{R}^{d})$ can be defined as

[TABLE]

where $\widehat{M}(\omega)=\frac{2^{1-\beta}}{\Gamma(\beta)}(\delta\|\omega\|_{2})^{\beta-d/2}K_{d/2-\beta}(\|\omega\|_{2}/\delta)$ and $K_{v}$ is the modified Bessel functions; and similarly to Gaussian kernels, for any $\Omega\subseteq\mathbb{R}^{d}$ and all $k\geqslant 0$,

[TABLE]

where $C_{M}>0$ depends only on $\beta$ and $d$, see Theorem $7.6$ of [1] for details.

3 Sparse residual tree

Sparse residual tree is based on both tree decomposition and adaptive RBF explorations. Suppose $\epsilon_{\mathrm{E}}>0$ is the expected relative absolute error (RAE) for an approximation $s$ of the target function $f$ on the interpolation dataset $X$, where

[TABLE]

For each child, for example, $X^{\prime}\subset\Omega^{\prime}$, which is $X\subset\Omega$ itself in the beginning, we need to (i) explore a sparse and proper RBF approximation $s_{r^{\prime}}$ to minimize the $2$-norm of the current residual $r^{\prime}(X^{\prime})$; and then, (ii) split the dataset $X^{\prime}$ into two proper subsets $X^{\prime}_{1}$ and $X^{\prime}_{2}$ as well as the domain $\Omega^{\prime}$ into two proper subdomains $\Omega^{\prime}_{1}$ and $\Omega^{\prime}_{2}$, as shown in the following diagram. We call it an exploration-splitting process. {diagram}

As mentioned above, the RBF exploration applies only to a sparse but sufficient subset of $X^{\prime}$ and the centers of $s_{r^{\prime}}$ are also sparse relative to the sparse subset. Hence, let us start with a sparsification of the dataset $X^{\prime}$ when its number is large.

3.1 Sparsification of datasets

Except for updating the residual, we hope to improve efficiency by replacing $X^{\prime}$ with its subset which has the same distribution of $X^{\prime}$ if the number of $X^{\prime}$ is large. Since $s_{r^{\prime}}$ is only used to refine the relative global component of the current residual $r^{\prime}$, it is not necessary to use all the data. Let $I$ be an index vector containing $N^{\prime}_{I}(\leqslant N^{\prime})$ unique integers selected randomly from $1$ to $N^{\prime}$ inclusive, then $X^{\prime}(I)$ is exactly what we need. Actually, from the independence of $X^{\prime}$ and $I$, it follows that

[TABLE]

i.e., $X^{\prime}(I)$ has the same probability distribution of $X^{\prime}$. And the choice of the number $N^{\prime}_{I}$ will be discussed later.

3.2 Quasi-uniform subsequence

Now we consider a method for generating a quasi-uniform subsequence of $X^{\prime}(I)$, which is the basis for adaptive RBF explorations. To find a quasi-uniform subsequence $U^{\prime}$ from $X^{\prime}(I)=\{x^{\prime}_{I(1)},\cdots,x^{\prime}_{I(N^{\prime}_{I})}\}$, we start with the approximate mean point, that is,

[TABLE]

And for known $U^{\prime}_{j}=\{u^{\prime}_{1},\cdots,u^{\prime}_{j}\}$, the subsequent point $u^{\prime}_{j+1}$ is determined as

[TABLE]

i.e., $u^{\prime}_{j+1}\in X^{\prime}(I)$ is the point that maximizes the minimum of the set of distances from it to a point in $U^{\prime}_{j}$. By storing an $N^{\prime}_{I}$-dimensional distance vector and an $N^{\prime}_{I}$-dimensional index vector, it only takes $\mathcal{O}(jN^{\prime}_{I})$ operations to generate $j$ quasi-uniform points $U^{\prime}_{j}$ and determine the relationship between every point of $X^{\prime}(I)$ and the Voronoi diagram of $U^{\prime}_{j}$, see Fig. 1 for examples.

3.3 Adaptive RBF exploration

The purpose of this adaptive exploration is to determine the centers of the RBF refinement $s_{r^{\prime}}$ which is only used to refine the relative global component of $r^{\prime}$. We first introduce the working parameters of SRT

[TABLE]

where $\omega_{1}=\kappa>0$ is the upper bound of condition numbers, $\omega_{2}>0$ is the termination error of explorations, $\omega_{3}\in(0,1)$ is the factor of shape parameters, and $\omega_{4}>0$ is the termination factor of tree nodes. For a fixed factor $\omega_{3}$, the current shape parameter can be determined as

[TABLE]

and the meaning of the remaining parameters will be clarified more clearly later.

Suppose $\{\chi^{\prime}_{1},\cdots,\chi^{\prime}_{j^{\prime}}\}$ are the centers inherited from its father, a reasonable idea is to choose the $(j^{\prime}+1)$th center $\chi^{\prime}_{j^{\prime}+1}$ from the quasi-uniform subsequence $U^{\prime}_{j^{\prime}+d+1}$ which is generated by (8) with the initial $U^{\prime}_{j^{\prime}}=\{\chi^{\prime}_{1},\cdots,\chi^{\prime}_{j^{\prime}}\}$; for the root node, we choose $\chi^{\prime}_{1}=u^{\prime}_{1}\in U^{\prime}_{1}$ given by (7). Without loss of generality, for known $\{\chi^{\prime}_{1},\cdots,\chi^{\prime}_{j}\}\subset U^{\prime}_{j+d}$ with $j\geqslant j^{\prime}$, we determine the $(j+1)$th center $\chi^{\prime}_{j+1}$ from $U^{\prime}_{j+d+1}-\{\chi^{\prime}_{1},\cdots,\chi^{\prime}_{j}\}$ by the following procedure:

From the recursive QR decomposition (as mentioned in section 2), $R_{j}$ and $Q_{j}^{\mathrm{T}}r^{\prime}(X^{\prime}(I))$ can be recursively obtained by $R_{j-1}$ and $Q_{j-1}^{\mathrm{T}}r^{\prime}(X^{\prime}(I))$ without computing $Q_{j-1}$, where

[TABLE]

and $\Phi_{X^{\prime}(I),\{\chi^{\prime}_{1},\cdots,\chi^{\prime}_{v}\}}$ is generated by the Gaussian kernel $G_{\delta^{\prime}}$. 2. 2.

The temporary residual can be obtained by

[TABLE]

where the coefficients $\alpha_{j}=\left(\alpha^{(1)}_{j},\cdots,\alpha^{(j)}_{j}\right)^{\mathrm{T}}=R_{j}^{-1}Q_{j}^{\mathrm{T}}r^{\prime}(X^{\prime}(I))$. 3. 3.

Suppose $\{\Lambda_{l}\}_{l=1}^{j+d+1}$ be the Voronoi diagram of the set $U^{\prime}_{j+d+1}$ and $\{\Lambda_{l}\}_{l\in\Gamma}$ are Voronoi regions with respect to those elements from the complementary set $U^{\prime}_{j+d+1}-\{\chi^{\prime}_{1},\cdots,\chi^{\prime}_{j}\}$, then

[TABLE]

where

[TABLE]

and $n_{l}$ is the point number of $\Lambda_{l}\cap X^{\prime}(I)$. 4. 4.

And the termination criteria is

[TABLE]

where $\kappa(R)=\frac{\max_{l}|R_{ll}|}{\min_{l}|R_{ll}|}$ is an estimation of the condition number $\|R^{-1}\|\|R\|$ and

[TABLE]

To ensure that the centers is not too sparse, its number should usually be greater than $d+2$ (imagine a case that the domain $\Omega^{\prime}$ is a $d$-dimensional simplex). Obviously, each newly selected center is in the Voronoi region with the largest mean squared error of the temporary residual. This allows the exploration to effectively capture the global component of the residual, see Fig. 2 for examples.

The sparse RBF refinement $s_{r^{\prime}}$ is obtained when the exploration is terminated, then we update the residual $r^{\prime\prime}$ on the full set $X^{\prime}$. Let the final number of the centers is $N^{\prime\prime}$ and relevant coefficients $\alpha^{\prime}=\left(\alpha^{\prime}_{1},\cdots,\alpha^{\prime}_{N^{\prime\prime}}\right)^{\mathrm{T}}$, then

[TABLE]

In addition, assume that the number of all currently existing nodes is $M$ and $\{n_{c}^{(i)}\}_{i=1}^{M}$ is the set of the center number of each node, now define the average

[TABLE]

and we can use a certain multiple of the average $\bar{n}_{c}$, say $100$ times, as the value of $N^{\prime}_{I}$ for the sparsification of the next node. For the initial node we usually take a fixed value related to the dimension $d$.

3.4 Equal binary splitting and termination

First we consider the selection of two splitting points, then use a hyperplane, whose normal is defined by these two points, to split all the points $X^{\prime}$ into two parts as well as the domain $\Omega^{\prime}$ into two subdomains. Clealy, since the half space and $\Omega^{\prime}$ are both convex, each subdomain is also convex. In order to block the spread of error, we expect to separate the points with large errors from those with small errors. First, we generate $d+1$ quasi-uniform points $U^{\prime}_{d+1}$ of $X^{\prime}(I)$ by the method of subsection 3.2 with a different starting point:

[TABLE]

Assume that the domain $\Omega^{\prime}$ is a $d$-dimensional simplex and $X^{\prime}$ is dense enough, then $U^{\prime}_{d+1}$ can almost be viewed as its vertices. Let $\{\Lambda_{l}\}_{l=1}^{d+1}$ be the Voronoi diagram of $U^{\prime}_{d+1}$, then the first splitting point is determined as

[TABLE]

where $l^{*}=\arg\max_{l}\sum_{x^{\prime}\in\Lambda_{l}\cap X^{\prime}(I)}\frac{|r^{\prime\prime}(x^{\prime})|^{2}}{n_{l}}$ and $n_{l}$ is the point number of $\Lambda_{l}\cap X^{\prime}(I)$. Then the second splitting point is determined as

[TABLE]

Then, according to the projections of $X^{\prime}$ in the direction $x^{\prime}_{b}-x^{\prime}_{a}$ and its median, $X^{\prime}$ can be splitted into $X^{\prime}_{1}$ and $X^{\prime}_{2}$ with the sizes $\lceil\frac{N^{\prime}}{2}\rceil$ and $N^{\prime}-\lceil\frac{N^{\prime}}{2}\rceil$, respectively; where $\lceil t\rceil$ denotes the least integer greater than or equal to $t$. Specifically, let $\vec{n}^{\prime}=(x^{\prime}_{b}-x^{\prime}_{a})^{\mathrm{T}}$, then the projections

[TABLE]

let $c^{\prime}=\textrm{median}(P_{\vec{n}^{\prime}}(X^{\prime}))$, then $X^{\prime}_{1}$ and $X^{\prime}_{2}$ can be given as

[TABLE]

and similarly, $\Omega^{\prime}_{1}$ and $\Omega^{\prime}_{2}$ can be given as

[TABLE]

Since the local high-frequency error tends to propagate over the entire domain, blocking its propagation is very important for a sparse approximation, and this is the motivation for designing the above splitting, see Fig. 3.

This exploration-splitting process finally stops if the expected RAE $\epsilon_{\mathrm{E}}$ is reached or the data is insufficient at the current tree node. Another important use of the average $\bar{n}_{c}$ defined in (13) is to determine whether the data is sufficient. Obviously, a sparse approximation must be based on relatively sufficient data, so if the size of $X^{\prime}_{1}$ or $X^{\prime}_{2}$ is less than $\omega_{4}$ times the average $\bar{n}_{c}$ and the RAE of residual still does not reach the expected $\epsilon_{\mathrm{E}}$, then we consider that the relevant node is lack of data, terminate further operations and record the node. A proper $\omega_{4}$ can guarantee that the prediction does not over-fit the data.

3.5 SRT prediction and its error characteristics

Suppose $s^{\prime}$ is the current approximation on the domain $\Omega^{\prime}$ and $s^{\prime}_{r^{\prime\prime}(X^{\prime}_{i})}$ is the refinement on $\Omega^{\prime}_{i}~{}(i=1,2)$. Then the next approximation $s^{\prime}_{i}$ on $\Omega^{\prime}_{i}~{}(i=1,2)$ can be given as

[TABLE]

It is clear that the SRT prediction is actually piecewise smooth on the original domain $\Omega$, hence the error of each piece will be significantly larger near the boundary.

The following example illustrates the error characteristics of SRTs. Although the SRT prediction, as shown on the left-hand side of Fig. 4, can adaptively build a piecewise and sparse approximation according to local features of the target function, the approximation error, as shown on the right-hand side of Fig. 4, may be significantly larger near the boundary of each piece. Hence, we will introduce the sparse residual forest for overcoming this boundary effect of the error in the next section.

The partition of unity is also one of the methods to address this issue. By introducing appropriate overlapping domains and rapidly decaying weight functions, the boundary effect of the error can be alleviated to some extent. However, since the overlapping domains usually cannot be too small and the depth of the tree is often not small, its time and space costs are significantly higher than $\mathcal{O}(N\log_{2}N)$. Instead, sparse residual forests still have the same cost as SRTs. And it provides even better performance than the partition of unity based method in terms of accuracy.

4 Sparse residual forest

Sparse residual forest (SRF) is a combination of SRT predictors with different tree decompositions. It provides an opportunity to avoid those predictions near the boundary and then use the average value of the remaining predictions to enhance both stability and convergence. First, we introduce a random splitting for SRTs. It can help generate random tree decompositions.

4.1 Random binary splitting

To get a random splitting, we only need to replace the median with a random percentile in (14). Let $p_{r}$ be a randomly selected integer from $37$ to $62$ inclusive, then $c^{\prime}$ can be redefined as

[TABLE]

where $\textrm{percentile}(Z,p_{r})$ denotes the percentile of the values in a data vector $Z$ for the percentage $p_{r}$. Note that $0.618$ is the golden ratio and this method depends on the values of a random vector sampled independently and with the same distribution.

4.2 SRF prediction

Suppose $n_{t}$ is the number of SRTs in the SRF, we usually apply the equal splitting to generate the first SRT and the random splitting to create the remaining $n_{t}-1$ SRTs. SRF helps us to avoid those predictions with large squared deviations and to use the average value of the remaining predictions to enhance both stability and convergence.

For any $x\in\Omega$, let $s_{\textrm{SRT}}^{(i)}(x)$ be the $i$th SRT prediction ($1\leqslant i\leqslant n_{t}$), then the squared deviation

[TABLE]

further, let the indicator set

[TABLE]

then the SRF prediction

[TABLE]

The indicator set $I_{F}$ here is used to avoid those predictions near the boundaries. In practice, as shown in Fig. 5, SRFs composed of a small number of SRTs perform quite well than individual SRTs; and in theory, similar to random forests [21], the error for SRFs converges with probability $1$ to a limit as $n_{t}$ becomes large, see Fig. 6 for examples and subsection 5.3 for details.

Although SRF predictions usually have smaller errors when the SRT number $n_{t}$ is larger, we usually do not recommend choosing a large $n_{t}$, which means $n_{t}$ times the storage and computational cost.

5 Theory

5.1 Stability properties

Suppose $\Omega_{L-1}$ is a leaf node, that is at the lowest level in a SRT, and $L$ levels of approximation, then there exists a domain sequences $\Omega_{0}\supset\Omega_{1}\supset\cdots\supset\Omega_{L-1}$ and a relevant dataset sequences $X_{0}\supset X_{1}\supset\cdots\supset X_{L-1}$ with relevant sizes $N_{0}>N_{1}>\cdots>N_{L-1}$ and shape parameters $\delta_{0}<\delta_{1}<\cdots<\delta_{L-1}$, where $\Omega_{0}=\Omega$ is convex, $X_{0}=X$ and $N_{0}=N$; and then, the SRT prediction of the target function $f$ is

[TABLE]

and the final residual

[TABLE]

where $s_{i}(x)=\sum_{j=1}^{N^{\prime}_{i}}\alpha_{i}^{(j)}G_{\delta_{i}}(x-\chi_{i}^{(j)})\in\mathcal{N}_{G_{\delta_{i}}}(\Omega_{i})$ is the LS approximation of the residual $r_{i}(X_{i})$ with respect to the centers $\chi_{i}=\{\chi_{i}^{(j)}\}_{j=1}^{N^{\prime}_{i}}\in X_{i}$, and $r_{i+1}=r_{i}-s_{i}$ with $r_{0}=f$. Then, for any $1\leqslant i\leqslant L-1$, it follows that

[TABLE]

where $Q_{i}R_{i}$ is the QR decomposition of the current matrix $\Phi_{i}=\Phi_{X_{i},\chi_{i}}$ generated by the kernel $G_{\delta_{i}}$. If $\tau_{i}$ is the smallest singular value of $R_{i}$, then

[TABLE]

According to the orthogonality of $s_{i}(X_{i})$ and $r_{i+1}(X_{i})$, we can obtain the following recurrence relations

[TABLE]

and

[TABLE]

thus, it follows that

[TABLE]

Together with (20), we proved the following theorem.

Theorem  5.1**.**

Suppose $s_{\mathrm{SFT}}$ is a SRT prediction of a function $f$ on a leaf node $\Omega_{L-1}\subset\Omega$ with respect to the data $(X,f_{X})$, as defined in (18). Let $\alpha_{i}$ be the coefficients of the $i$th level LS approximation $s_{i}$, then

[TABLE]

where $\tau=\min_{1\leqslant i\leqslant L-1}\tau_{i}$ and the constants $\tau_{i}$ comes from (20).

Note that this theorem obviously holds for our SRTs with sparsification processes introduced in subsection 3.1. And now we can prove the following theorem.

Theorem  5.2**.**

Under the supposition of Theorem 5.1. For all $1\leqslant p<\infty$, $k\in\mathbb{N}_{0}$, $\delta\geqslant\delta_{L-1}$, and any leaf node $\Omega_{L-1}$ of the prediction $s_{\mathrm{SRT}}$, it holds that

[TABLE]

where the constant $\tau$ comes from Theorem 5.1, the constant $C_{W}$ depends only on $\delta_{0},\delta_{L-1},d,p$ and $k$, and the constant $C_{\mathcal{N}}$ depends only on $\delta_{L-1}$ and $d$.

Proof.

To prove the first inequality, observe that

[TABLE]

where $M_{\delta_{i}}^{p,k}=\left(\sum_{|r|=k}\|D^{r}G_{\delta_{i}}\|^{p}_{L_{p}(\mathbb{R}^{d})}\right)^{1/p}$, and for any $0\leqslant i\leqslant L-1$, $M_{\delta_{i}}^{p,k}<M_{\delta_{L-1}}^{p,k}$ when $k>1$; or $M_{\delta_{i}}^{p,k}<M_{\delta_{0}}^{p,k}$ when $k<1$; or $M_{\delta_{i}}^{p,k}=M^{p,k}$ is independent of $\delta_{i}$ when $k=1$. Together with Theorem 5.1, we have

[TABLE]

where $C_{W}=C_{1}M_{\delta_{L-1}}^{p,k}$ when $k>1$, or $C_{W}=C_{1}M_{\delta_{0}}^{p,k}$ when $k<1$, or $C_{W}=C_{1}M^{p,k}$ when $k=1$.

To prove the second inequality, observe that for any $s_{i}\in\mathcal{N}_{G_{\delta_{i}}}(\Omega_{i})$, there is a natural extension $\mathcal{E}s_{i}\in\mathcal{N}_{G_{\delta_{i}}}(\mathbb{R}^{d})$ with $\|\mathcal{E}s_{i}\|_{\mathcal{N}_{G_{\delta_{i}}}(\mathbb{R}^{d})}=\|s_{i}\|_{\mathcal{N}_{G_{\delta_{i}}}(\Omega_{i})}$. From the definition of native spaces of Gaussians, we see that $\mathcal{E}s_{i}\in\mathcal{N}_{G_{\delta}}(\mathbb{R}^{d})$ with

[TABLE]

where $\delta\geqslant\delta_{L-1}>\cdots>\delta_{0}$; and further, the restriction $\mathcal{E}s_{i}|\Omega_{L-1}=s_{i}|\Omega_{L-1}$ of $\mathcal{E}s_{i}$ to $\Omega_{L-1}\subseteq\Omega_{i}$ is contained in $\mathcal{N}_{G_{\delta}}(\Omega_{L-1})$ with

[TABLE]

hence, we have $\|s_{i}|\Omega_{L-1}\|_{\mathcal{N}_{G_{\delta}}(\Omega_{L-1})}\leqslant\|\mathcal{E}s_{i}\|_{\mathcal{N}_{G_{\delta_{i}}}(\mathbb{R}^{d})}$, and then

[TABLE]

Together with Theorem 5.1 and

[TABLE]

we finally have $\|s_{\mathrm{SFT}}\|_{\mathcal{N}_{G_{\delta}}(\Omega_{L-1})}<C_{\mathcal{N}}\cdot\tau^{-1}\cdot\|f_{X}\|_{2}$, where $C_{\mathcal{N}}=C_{2}(2\delta_{L-1})^{d/2}\pi^{d/4}$. ∎

Remark  5.1**.**

See Theorems $10.46$ and $10.47$ in [2] for details about the restriction and extension of functions from certain native spaces.

Remark  5.2**.**

The second inequality depends on the embeddings (21) of native spaces of Gaussians. As mentioned in section 2, the Fourier transform of the inverse multiquadrics is $\widehat{M}_{\delta}(\omega)=\frac{2^{1-\beta}}{\Gamma(\beta)}(\delta\|\omega\|_{2})^{\beta-d/2}K_{d/2-\beta}(\|\omega\|_{2}/\delta)$, then for any $\beta>\frac{d}{2}$ and $\delta\geqslant\delta_{i}$, $\widehat{M}_{\delta}^{-1}(\omega)\leqslant\widehat{M}_{\delta_{i}}^{-1}(\omega)$, and then, for an inverse multiquadric based $s_{i}$,

[TABLE]

hence, the second inequality also holds for native spaces of inverse multiquadrics.

5.2 Error estimates for SRTs

Theorem  5.3**.**

Under the supposition of Theorem 5.1. If $f\in W_{p}^{k}(\Omega)$ and $r_{L}$ is the residual $f-s_{\mathrm{SRT}}$ on an arbitrary leaf node $\Omega_{L-1}$, then for any $1\!\leqslant\!q\!\leqslant\!\infty$, $\gamma\in\mathbb{N}_{0}^{d}$, and $1\!\leqslant\!p\!<\!\infty$ with $k>|\gamma|+d/p$ if $p>1$, or with $k\geqslant|\gamma|+d$ if $p=1$, it holds that

[TABLE]

where $(t)_{+}=\max(t,0)$, the fill distance $h$ is assumed to be sufficiently small, the constant $C$ do not depend on $f,r_{L}$ or $h$, and the constant $C_{W}$ comes from Theorem 5.2.

Proof.

According to the sampling inequality for functions from certain Sobolev spaces on a bounded domain (see Theorem $2.6$ in [3]), we have

[TABLE]

and further,

[TABLE]

Applying the first inequality of Theorem 5.2 finishes the proof. ∎

This result also explains how the matrix $\Phi_{i}$ at each level affects the convergence. It is worth noting that this proof does not depend on the radial basis functions, so the next observation is an immediate consequence.

Corollary  5.1**.**

The result of Theorem 5.3 holds for arbitrary basis functions based SRTs provided those basis functions belongs to $W_{p}^{k}(\Omega)$.

It shows that a SRT, whose basis functions are differentiable and have bounded derivatives on $\Omega$ (regardless of polynomials, trigonometric polynomials, radial basis functions), leads to algebraic convergence orders for finitely smooth target functions. For infinitely smooth target functions, the following theorem shows that the Gaussian based SRT leads to exponential convergence orders.

Theorem  5.4**.**

Under the supposition of Theorem 5.1. If $f\in\mathcal{N}_{G_{\delta}}(\Omega)$ and $r_{L}$ is the residual $f-s_{\mathrm{SRT}}$ on an arbitrary leaf node $\Omega_{L-1}$, then for any $1\!\leqslant\!q\!\leqslant\!\infty$, $\gamma\in\mathbb{N}_{0}^{d}$, and $\delta>\delta_{L-1}$, there are constants $C$ and $h_{0}$ such that for all $h\leqslant h_{0}$, it holds that

[TABLE]

where the constant $C$ depends only on the geometry of $\Omega_{L-1}$, $h_{0}$ may depend on $d,p,q,\gamma$ and the geometry of $\Omega_{L-1}$ but not on $h$ or $f$, $C^{\prime}$ do not depend on $h$ or $r_{L}$, and the constant $C_{\mathcal{N}}$ comes from Theorem 5.2.

Proof.

According to the sampling inequality for functions from certain native spaces of Gaussians on a bounded domain (see Theorems $3.5$ and $7.5$ in [1]), we have

[TABLE]

and further,

[TABLE]

where $f|\Omega_{L-1}$ is the restriction of $f$ to $\Omega_{L-1}$ with $\|f|\Omega_{L-1}\|_{\mathcal{N}_{G_{\delta}}(\Omega_{L-1})}\leqslant\|f\|_{\mathcal{N}_{G_{\delta}}(\Omega)}$ (see Theorem $10.47$ in [2]); and applying the second inequality of Theorem 5.2 finishes the proof. ∎

Similarly, according to Remark 5.2 and the sampling inequality for functions from certain native spaces of Gaussians on a bounded domain (see Theorems $3.5$ and $7.6$ in [1]), we can also prove the convergence for the inverse multiquadric based SRTs.

Theorem  5.5**.**

Under the supposition of Theorem 5.1. If $f\in\mathcal{N}_{M_{\delta}}(\Omega)$, $s_{\mathrm{SRT}}$ is based on inverse multiquadrics, and $r_{L}$ is the residual $f-s_{\mathrm{SRT}}$ on an arbitrary leaf node $\Omega_{L-1}$, then for any $1\!\leqslant\!q\!\leqslant\!\infty$, $\gamma\in\mathbb{N}_{0}^{d}$, and $\delta>\delta_{L-1}$, there are constants $C$ and $h_{0}$ such that for all $h\leqslant h_{0}$, it holds that

[TABLE]

where the constants $C$ and $h_{0}>0$ depends only on $d,p,q,\gamma$ and the geometry of $\Omega_{L-1}$, $C^{\prime}$ do not depend on $h$ or $r_{L}$, and the constant $C_{\mathcal{N}}$ comes from Theorem 5.2.

5.3 Error estimates for SRFs

For any $x\in\Omega_{L-1}\subset\Omega$, each SRT prediction $s_{\mathrm{SRT}}^{(i)}(x)$ ($1\leqslant i\leqslant n_{t}$) in a SRF converges to the target function $f(x)$ and satisfies relevant error estimates, thus, together with the Strong Law of Large Numbers and the Lindeberg-Levy central limit theorem, it follows that:

Theorem  5.6**.**

For any $x\in\Omega_{L-1}$, there exists an expectation $m_{\mathrm{SRF}}(x)$ such that

[TABLE]

converges almost surely to $m_{\mathrm{SRF}}(x)$. Further, for any $1\leqslant q\leqslant\infty$ and $\gamma\in\mathbb{N}_{0}^{d}$, if $\|D^{\gamma}(s_{\mathrm{SRT}}^{(i)}(x)-f(x))\|_{L_{q}(\Omega_{L-1})}\leqslant\epsilon$, then there exists $\sigma\leqslant 2\epsilon$ such that the random variables $\|D^{\gamma}(s_{\mathrm{SRF}}(x;n_{t})-m_{\mathrm{SRF}}(x))\|_{L_{q}(\Omega_{L-1})}$ converge in distribution to a normal $N(0,\sigma/\sqrt{n_{t}})$, i.e., for any $\lambda_{a}>0$, the inequality

[TABLE]

holds with probability $1-a$, where $a=\frac{1}{\sqrt{2\pi}}\int_{-\lambda_{a}}^{\lambda_{a}}e^{-\frac{t^{2}}{2}}\mathrm{d}t$.

Obviously, the above result also holds for the SRF prediction defined in (17) that is more stable and is specially designed for overcoming the boundary effect of the error, as shown in Fig. 3. Combining the results of the previous subsection, one can obtain the error estimates for SRF predictions in the corresponding spaces.

5.4 Complexity analysis

Since the maximum depth of a binary tree is $\log_{2}N$ and the full data is only used for updating the residual, it is easy to see that:

Theorem  5.7**.**

Algorithm in section 3 needs $\mathcal{O}(N\log_{2}N)$ time and $\mathcal{O}(N\log_{2}N)$ space in the worst case to train a SRT for $N$ arbitrary distributed points; and needs $\mathcal{O}(\log_{2}N)$ time in the worst case to make a prediction for a new point $x$. And the costs of algorithm in section 4 are $n_{t}$ times that of the SRT for a SRF with $n_{t}$ SRTs.

This result shows that the SRT or SRF also yields the excellent performance in terms of efficiency in addition to accuracy and adaptability. It is worth pointing out that the algorithm in section 3 is designed to achieve hierarchical parallel processing so that the training process can be accelerated using multi-core architectures.

6 Numerical examples

In this section we compare the performance of both SRT and SRF with the Gaussian process regression (GPR). For an approximation $s$ of the target function $f$ on a test dataset $Z=\{z_{i}\}_{i=1}^{N_{t}}$ of size $N_{t}$, we use the relative mean absolute error (RMAE) as a measure of accuracy, i.e.,

[TABLE]

We use two test functions: one is Franke’s function, which is defined as:

[TABLE]

where $x\in[0,1]^{d}$ for $d\geqslant 2$; and the other is local oscillating and defined as:

[TABLE]

where $x\in[-7,7]^{d}$ for $d\geqslant 2$.

All our numerical tests are based on scattered data which are either randomly generated or the Halton sequence [23]. In addition, the procedure for the above two methods at each sample size is repeated $5000$ times for investigating the stability of the results. We use Matlab’s function fitrgp to generate a GPR model trained using the same sample data of proposed methods. Fit the GPR model using the subset of regressors method for parameter estimation and fully independent conditional method for prediction. Standardize the predictors. Besides, since the computational complexity of GPR is $\mathcal{O}(N^{3})$ for training work and $\mathcal{O}(N^{2})$ for each prediction, where $N$ is the sample size, it is very difficult to use GPR for large data set, so the sample size is varied from $10^{1}$ to $10^{4}$ for all numerical tests by using GPR.

6.1 Accuracy, sparsity, storage and computational time

The size of Halton points are varied from $10^{1}$ to $10^{4}$ for Franke’s function. The results are shown in Fig. 7. From the upper left of Fig. 7, as expected, the RMAE of both SRT and SRF are much lower than GPR as data point $N$ is large. Moreover, from the upper right of Fig. 7, we can find out that the average number of centers for SRT prediction at $1$ point is varying from $10$ to $182$. Besides, since the size of sample points is varied from $10^{1}$ to $10^{4}$, both the storage requirement and the computational time of the proposed methods are much lower than those of GPR.

6.2 Insufficient data report

We choose the second test function $g(x)$, $x\in[-7,7]^{2}$, to illustrate the insufficient data situation. From Fig. 8 we can find out that since $g(x)$ is complicated near the central of domain and the RAE of residual still does not reach the expected error when the sample points $N=3000$; that is, the relevant node is lack of data at this area. Further, by adding the size of sample points $N$ to $6000$, as expected, from the lower left of Fig.8 we find out that the RAE clearly decreased (with the maximum RAE decreases from $0.1926$ to $0.0784$). Besides, from the lower right of Fig. 8, the median value of centers for SRT prediction at one point for both $N=3000$ and $N=6000$ are close to $110$.

6.3 $3$-dimensional problem

The $3$-dimensional Franke’s function can be shown in Fig. 9. The size of Halton points are varied from $10^{1}$ to $10^{6}$. We find out that the RMAEs of proposed methods are not as good as that of GPR when size of sample points is less than $10^{4}$. Further, since the size of sample points is varied from $10^{4}$ to $10^{6}$, leading to the value of error varying from $5.7224\times 10^{-4}$ to $2.3126\times 10^{-7}$ by using SRT, and from $1.3037\times 10^{-4}$ to $4.7757\times 10^{-8}$ by using SRF. Besides, from low figures of Fig. 9, it is noted that both the storage requirement and the computing time of proposed methods are less than those of GPR. The average number of centers for SRT prediction at one point is varying from $10$ to $853$.

7 Conclusions

In this work, we proposed two new methods for multivariate scattered data approximation, named Sparse residual tree (SRT) and Sparse residual tree (SRF), respectively. We proved that the time complexity of SRTs is less than $\mathcal{O}(N\log_{2}N)$ for the initial work and $\mathcal{O}(\log_{2}N)$ for each prediction, and the storage requirement is less than $\mathcal{O}(N\log_{2}N)$, where $N$ is the data points. From the numerical experiments, we can find out that the proposed methods are good at dealing with cases where the data is sufficient or even redundant. For the higher dimensional problem, the proposed methods do not work as well as we expected. The possible reason is that the sample size is usually difficult to be sufficient or even redundant for higher dimensional problems, and the proposed methods tend to point out the possible local regions where data refinement is needed, rather than obtain approximations. It provides that the proposed methods can be used to solve the large data sets problems. In the following works, we will try to improve the proposed methods for solving higher dimensional problems.