We study the kernel-based multilevel method for approximating or learning a multivariate function from scattered data, motivated in part by recent applications in sparse grid methods. For nested families of point sets, we derive a nodal representation of the multilevel interpolant that depends explicitly on the values of the target function. This representation allows us to characterize the range of the associated interpolation operator and to construct a cardinal basis for this space. We prove that the resulting basis functions exhibit exponential decay, analogous to the localization properties known for certain kernel-based Lagrange functions. We further analyze the computational cost of the resulting formulation. For non-nested families of point sets, we derive a generalized nodal representation.
Radial Basis Functions, multilevel method, Lagrange functions
46E35, 65D05 65D12, 65D15, 65J10
Kernel-based methods have become a well-established tool in scattered data approximation [1]–[3], meshless numerical methods for partial differential equations [4]–[6] and machine learning [7]–[10]. For a broader overview see also [11]. They are particularly valued for their flexibility, stability and deterministic error analysis [12]–[14].
A major difficulty arises when the number of data becomes large: the kernel interpolation matrix may become severely ill-conditioned [14]. A common remedy is to introduce a scaling parameter and decrease the correlation length if the number of points grows. If chosen appropriately, this leads to well-conditioned interpolation matrices. However, then, convergence of the approximation method is lost. This phenomenon is known as the trade-off principle: either the approximation method converges but the system matrix becomes ill-conditioned, or the system matrix is well-conditioned but convergence is no longer guaranteed [14], [15].
For compactly supported kernels, ideas have been developed to overcome this trade-off principle by representing multiple scales within the data [15], [16]. The algorithm follows a residual correction strategy. First, an approximation is computed on a coarse set of data sites using a large kernel support, capturing large-scale features of the target function. The residual between the target function and this approximation typically contains finer-scale information. This residual is then approximated on a finer set of points using a smaller kernel support. The procedure is iterated until the desired accuracy is reached. In the interpolation setting, the sum of the levelwise interpolants is again an interpolant of the target function on the finest set of points. Although the procedure was established early and was shown to converge numerically, theoretical results were established much later in [17], [18]. This method was also studied in [19]–[22] for globally supported kernels and was used for solving partial differential equations in, for example, [23]–[26].
Recently, this multilevel approach has been combined with Smolyak’s construction to obtain kernel methods for high-dimensional function reconstruction [27]–[30]. In this setting the multilevel operator is tensorized, which requires an explicit representation of its dependence on the data. However, the standard formulation expresses the operator implicitly through recursively defined residuals. Although first results in this direction were obtained in the aforementioned works, the resulting formulations are computationally infeasible in practice. In this paper we derive a nodal representation of the multilevel approximation operator that makes its dependence on the data explicit. This representation yields new structural insights into the global approximation space. In particular, for nested families of point sets we construct both a Lagrange basis and a Newton-type basis of the multiscale space. These bases allow us to characterize the range of the multilevel operator and to analyze localization properties and computational cost.
The analysis is more involved for non-nested sets of data sites. Nevertheless, the ideas developed for the nested case can still be used to construct a basis of the global approximation space and to characterize the range of the multilevel operator. In general, this basis is no longer cardinal. Such non-nested configurations arise naturally in adaptive variants of multilevel methods.
These results provide a new structural understanding of multilevel kernel methods and lead to practical representations of the approximation operator that are suitable for high-dimensional constructions.
The remainder of the paper is organized as follows. 2 introduces notation and reviews the multilevel approximation method. We also extend the convergence result of [18] to a broader range of Sobolev norms. 3 studies nested families of data sites and derives Lagrange and Newton-type bases for the global approximation space. Exponential decay of the Lagrange basis functions is established in 3.2. 4 treats the non-nested case. 5 provides a conclusion and outlook.
The setup of the kernel-based multilevel method is usually the following one, cf. [18]. For a bounded domain \(\Omega \subseteq \mathbb{R}^d\) we assume that we are given a sequence of finite point sets \(X_1, X_2, \ldots \subseteq \Omega\), where \(X_\ell = \left\{ \boldsymbol{x}_1^{(\ell)}, \dots, \boldsymbol{x}^{(\ell)}_{N_\ell} \right\}.\)
Associated to the set \(X_\ell\), \(\ell\in\mathbb{N}\), we define the fill distance \(h_\ell\) and separation radius \(q_\ell\) as usual by \[h_\ell := h_{X_\ell,\Omega} = \sup_{\boldsymbol{x}\in \Omega} \min_{1\le i\le N_\ell} \| \boldsymbol{x}- \boldsymbol{x}^{(\ell)}_{i} \|_2, \qquad q_\ell := q_{X_\ell} = \frac{1}{2} \min_{i \neq j}\| \boldsymbol{x}^{(\ell)}_i - \boldsymbol{x}^{(\ell)}_j \|_2.\] For some results it is important that the data sites are quasi-uniform, that is there exists a level-independent constant \(c_{qu}>0\) such that \[q_\ell \leq h_\ell \leq c_{qu} q_\ell, \qquad \ell \in \mathbb{N}.\] The sets \(X_\ell\) are not required to be nested, however, we assume that the fill distances decay uniformly, i.e., that there is a uniform refinement parameter \(0 < \mu < 1\) and a constant \(c \in (0,1]\), such that \[c \mu h_\ell \leq h_{\ell+1} \leq \mu h_\ell, \qquad \ell \in \mathbb{N}.\]
The second important ingredient for defining our method is a compactly supported, positive definite radial basis function (RBF) \(\Phi: \mathbb{R}^d \to \mathbb{R}\) with support in the closed unit ball. Using scaling parameter \(\delta_\ell := \nu h_\ell\) with an \(\ell\)-independent \(\nu>1\), we obtain a basis function for every level \(\ell \in \mathbb{N}\), given as \[\label{eq:RescaledKernel} \Phi_\ell := \Phi_{\delta_\ell} := \delta^{-d} \Phi(\cdot / \delta_\ell), \quad \ell \in \mathbb{N}.\tag{1}\]
Noting that \(\Phi_\ell\) has support in the closed ball about zero with radius \(\delta_\ell\), we have the local approximation spaces \[W_\ell := \operatorname{span} \{ \Phi_\ell (\cdot - \boldsymbol{x}^{(\ell)}_{i}) \; : \; 1\le i\le N_\ell \},\] and, for \(L\in N\), global approximation spaces \[\label{eq:DefApproximationSpace} V_L=W_1+\cdots +W_L.\tag{2}\] Under the above assumption, it can be shown that the sum in 2 is direct, see [17]. Moreover, we can give an alternative basis for these spaces. To be more precise for every level \(\ell\) and every \(1 \le i \le N_\ell\) there are functions \(\chi_i^{(\ell)}\in W_\ell\) satisfying \(\chi_i^{(\ell)}(\boldsymbol{x}_j^{(\ell)})=\delta_{ij}\). They are given by \(\left(\chi^{(\ell)}_{1}, \ldots, \chi^{(\ell)}_{N_\ell} \right) = \boldsymbol{r}_\ell^{\text{T}} M_\ell^{-1},\) where \(\boldsymbol{r}_\ell = (\Phi_\ell(\cdot - \boldsymbol{x}_{i}^{(\ell)}))_i^{\text{T}}: \Omega\to \mathbb{R}^{N_\ell}\) and \(M_\ell = (\Phi_\ell(\boldsymbol{x}_{i}^{(\ell)} - \boldsymbol{x}_{j}^{(\ell)})_{ij})\in\mathbb{R}^{N_\ell\times N_\ell}\) is the positive definite kernel matrix.
Definition 1. A basis of the global approximation space \(V_L\) is defined by \[\mathcal{B}_{\operatorname{full}}:=\{\chi_i^{(\ell)} : 1\le \ell \le L, 1\le i\le N_\ell\}.\]
Obviously, this is indeed a basis of \(V_L\), showing that the dimension of the global approximation space \(V_L\) is \(N_1+\cdots + N_L\). In the case of nested data sets \(X_\ell\subseteq X_{\ell+1}\) the basis \(\mathcal{B}_{\operatorname{full}}\) is a Newton-type basis, as we have for any \(1\le \ell\le L\) that \(\chi_i^{(\ell)}(\boldsymbol{x}_j^{(m)})=0\) for all points from \(X_m\), \(1\le m\le \ell\) but not for the later added points from \(X_{\ell+1}\) to \(X_L\), i.e. it satisfies the Newton property of a basis batch-wise.
It is well-known that the local approximation spaces \(W_\ell\) are not rich enough and that the global spaces \(V_L\) are redundant. Hence, it is important to properly choose a good subspace of \(V_L\), which is still rich enough but less redundant. We will achieve this by taking the image of a multilevel method, which are describing now.
The kernel-based multilevel method depicts one way of successively choosing local approximants \(s_\ell\in W_\ell\) to form a global approximant \(f_\ell=s_1+\cdots+s_\ell\in V_\ell\). To be more precise, the method starts with an initial approximation \(f_0=0\) and an initial error \(e_0=f\) and then proceeds for \(\ell\in\mathbb{N}\) by computing a local approximation \(s_\ell\in W_\ell\) to the error \(e_{\ell-1}\). Then, the global approximation and the error are updated via \[\begin{align} f_\ell & = f_{\ell-1}+s_\ell,\\ e_\ell &= e_{\ell-1}-s_\ell. \end{align}\] Numerically, the procedure stops after, say, \(L\) levels, but for theoretical results we may assume that \(L\) tends to infinity.
In this paper, we are particularly interested in interpolation as the means of determining the local approximation \(s_\ell\in W_\ell\). Thus, we set \(s_\ell=I_{X_\ell,\Phi_\ell} (e_{\ell-1})\) for any \(\ell \in \mathbb{N}\), where the interpolation operator \(I_{X_\ell, \Phi_\ell}: C(\Omega) \to W_\ell\) is defined by \[I_{X_\ell,\Phi_\ell}(g) = \sum_{i=1}^{N_\ell} g(\boldsymbol{x}_i^{(\ell)}) \chi_i^{(\ell)}, \qquad g \in C(\Omega).\]
Definition 2. With the notation and assumptions made above, the (interpolatory) multilevel operator* \(A_L: C(\Omega) \to V_L\) is defined by \[A_L(f) := \sum_{\ell=1}^{L} I_{X_\ell,\Phi_\ell}(e_{\ell-1}).\] Its image will be denoted by \(\widetilde{V}_L:=A_L(C(\Omega))\).*
The following simple observation will be crucial in our analysis.
Lemma 1. The operator \(A_L\) is linear and, for any \(f\in C(\Omega)\), \(A_L f\) interpolates \(f\) on \(X_L\), i.e. \(A_L(f)(\boldsymbol{x}) = f(\boldsymbol{x})\) for \(\boldsymbol{x}\in X_L\).
Proof. Linearity of the operator follows from an easy induction on the number of levels and the observation that the concatenation of linear operators is linear. The proof of the interpolation property can be found in, e.g., [3]. ◻
We end this section by extending the error result of [18] to more general Sobolev norms. While this is obviously a straight-forward generalization, it seems that it has not been stated before and is important for itself. It also allows us to specify the choice of the RBF \(\Phi\) in more details.
A Hilbert space \(H\) of functions \(f: \Omega \to \mathbb{R}\) is called a reproducing kernel Hilbert space (RKHS) if there is a kernel \(K: \Omega \times \Omega \to \mathbb{R}\) such that \(K(\cdot, \boldsymbol{x}) \in H\) and \(f(\boldsymbol{x}) = \langle f, K(\cdot, \boldsymbol{x}) \rangle_H\) for all \(\boldsymbol{x}\in \Omega\). \(K\) is then called the reproducing kernel of \(H\). If \(\Omega\subseteq\mathbb{R}^d\) has a Lipschitz boundary or if \(\Omega=\mathbb{R}^d\) then the Sobolev space \(H^\sigma(\Omega)\) with \(\sigma>d/2\) is a RKHS. If the Fourier transform of the RBF \(\Phi: \mathbb{R}^d \to \mathbb{R}\) satisfies \[\label{eq:AlgebraicDecay} c_1(1+ \| \boldsymbol{\omega}\|_2^2)^{-\sigma} \leq \widehat{\Phi} (\boldsymbol{\omega}) \leq c_2(1+ \| \boldsymbol{\omega}\|_2^2)^{-\sigma}, \qquad \boldsymbol{\omega}\in \mathbb{R}^d,\tag{3}\] with constants \(c_1, c_2 > 0\), then \(K: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}\) defined by \(K(\boldsymbol{x},\boldsymbol{y}) := \Phi(\boldsymbol{x}- \boldsymbol{y})\) is the reproducing kernel of \(H^{\sigma}(\mathbb{R}^d)\) if equipped with the norm \[\label{eq:NativeSpaceNorm} \| f \|_{\Phi} := \int_{\mathbb{R}^d} \frac{ | \widehat{f}(\boldsymbol{\omega})|^2}{\widehat{\Phi}(\boldsymbol{\omega})} \;d \boldsymbol{\omega}.\tag{4}\] It is easy to see that this norm is equivalent to the usual norm on \(H^{\sigma}(\mathbb{R}^d)\). Furthermore, a simple Fourier argument shows that also our scaled RBFs \(\Phi_\ell\) are reproducing kernels of \(H^\sigma(\mathbb{R}^d)\). However, then, the equivalence constants depend on \(\delta_\ell\), i.e. we have \(\widetilde{c}_1 \|f\|_{\Phi_{\delta_\ell}} \le \|f\|_{H^\sigma(\mathbb{R}^d)}\le \widetilde{c}_2 \delta_\ell^{-\sigma} \|f\|_{\Phi_{\delta_\ell}}\) for all \(f\in H^\sigma(\mathbb{R}^d)\). For more details, we refer to [18].
The next theorem is a sampling inequality, taken from [31], which is a slight generalization of the sampling inequality originally derived in [12].
Theorem 1. Let \(\Omega \subseteq \mathbb{R}^d\) be a bounded Lipschitz domain. Let \(q \in [1,\infty]\) and \(\sigma > d/2\). Then there exist constants \(h_0 > 0\) and \(C > 0\) such that for all \(X = \{\boldsymbol{x}_1, \dots, \boldsymbol{x}_N \} \subseteq \Omega\) with \(h = h_{X, \Omega} \leq h_0\) and all \(u \in H^{\sigma} (\Omega)\) with \(u|_X = 0\), the bound \[\| u \|_{W_q^{\tau}(\Omega)} \leq C h^{\sigma - \tau - d(\frac{1}{2} - \frac{1}{q})_+} \| u \|_{H^{\sigma}(\Omega)}\] holds for all \(0 \leq \tau \leq \widetilde{\sigma}\), where \(\widetilde{\sigma} = \widetilde{\sigma}_0 := \sigma - d (1/2 - 1/q)_+\) in the case \(\sigma \in \mathbb{N}\) and either \(q > 2\) and \(\widetilde{\sigma}_0 \in \mathbb{N}\) or \(q = 2\). Otherwise, \(\widetilde{\sigma} = \lceil \widetilde{\sigma}_0 \rceil - 1\). Moreover, \(\tau \in \mathbb{N}_0\) in the case \(q = \infty\).
Now we are in the position to state the easy extension of the convergence result for the multilevel method. As in the original result in [18], which only gives \(L_2\)-error bounds, the sequence of data sites \(X_\ell\) is not required to be quasi-uniform.
Theorem 2. Let \(\Omega \subseteq \mathbb{R}^d\) be a bounded Lipschitz domain. Let \(X_1, X_2, \ldots \subseteq\Omega\) be a sequence of point sets in \(\Omega\) with fill distances \(h_1, h_2, \dots\) satisfying \(c \mu h_\ell \leq h_{\ell+1} \leq \mu h_\ell\) for \(\ell\in\mathbb{N}\) with fixed \(\mu \in (0,1)\), \(c \in (0,1]\) and \(h_1\) sufficiently small. Let \(\sigma > d/2\) and \(\Phi\) be the reproducing kernel of \(H^{\sigma}(\mathbb{R}^d)\), i.e., its Fourier transform satisfies 3 , and let \(\Phi_\ell\) be defined by 1 with scale factor \(\delta_\ell = \nu h_\ell\). Assume \(1 /h_1 \geq \nu \geq \gamma / \mu\) with a fixed \(\gamma > 0\). Let \(q \in [1, \infty]\). Then there exist constants \(C,C_1 > 0\) such that for all \(f \in H^{\sigma}(\Omega)\), with \(\alpha = C_1 \mu^{\sigma - \tau - d (\frac{1}{2} - \frac{1}{q})_+}\), \[\label{eq:ErrorEstimate} \| f - A_L(f) \|_{W^{\tau}_q(\Omega)} \leq C \alpha^{L} \| f \|_{H^{\sigma}(\Omega)},\qquad{(1)}\] holds for all \(0 \leq \tau \leq \widetilde{\sigma}\), where \(\widetilde{\sigma} = \widetilde{\sigma}_0 := \sigma - d (1/2 - 1/q)_+\) in the case \(\sigma \in \mathbb{N}\) and either \(q > 2\) and \(\widetilde{\sigma}_0 \in \mathbb{N}\) or \(q = 2\). Otherwise, \(\widetilde{\sigma} = \lceil \widetilde{\sigma}_0 \rceil - 1\). Moreover, \(\tau \in \mathbb{N}_0\) in the case \(q = \infty\).
Proof. We employ a specific recursion derived in [18] Proof of Theorem 1. If \(E: H^{\sigma}(\Omega) \to H^{\sigma}(\mathbb{R}^d)\) denotes Stein’s extension operator for Sobolev functions (see [32]) then this recursion is given by \[\label{eq:RecursionEstimate} \| E e_\ell \|_{\Phi_{\ell+1}} \leq C_1 \mu^{\sigma} \| E e_{\ell-1} \|_{\Phi_\ell},\tag{5}\] where we use the norm of 4 induced by \(\Phi_{\ell+1}\) and \(\Phi_\ell\), respectively.
The error on level \(L\) satisfies \(e_L |_{X_L} = 0\). This allows us to use the sampling inequality from 1 to obtain, in complete analogy to [18], the bound \[\begin{align} \| e_L \|_{W^{\tau}_q(\Omega)} &\leq C h_{L}^{\sigma - \tau - d (\frac{1}{2} - \frac{1}{q})_+} \| e_L \|_{H^{\sigma} (\Omega)} \\ &\leq C h_{L}^{\sigma - \tau - d (\frac{1}{2} - \frac{1}{q})_+} \| Ee_L \|_{H^{\sigma} (\mathbb{R}^d)} \\ &\leq C h_{L}^{\sigma - \tau - d (\frac{1}{2} - \frac{1}{q})_+} \delta^{-\sigma}_{L+1} \| Ee_L \|_{\Phi_{L+1}} \\ & = C \nu^{-L} h_L^{- \tau - d (\frac{1}{2} - \frac{1}{q})_+} \| Ee_L \|_{\Phi_{L+1}} \\ &\leq C C_1^L (\mu^L)^{- \tau - d (\frac{1}{2} - \frac{1}{q})_+} (\mu^L)^{\sigma} \| Ef \|_{\Phi_{1}} \\ &= C \left(C_1 \mu^{\sigma - \tau - d (\frac{1}{2} - \frac{1}{q})_+} \right)^L \| Ef \|_{\Phi_1}, \end{align}\] where we used the recursion 5 \(L\)-times. The stated estimate then follows from \(\| Ef \|_{\Phi_1} \leq C \| Ef \|_{H^{\sigma}(\mathbb{R}^d)} \leq C \| f \|_{H^{\sigma}(\Omega)},\) using the norm equivalence and the properties of the extension operator \(E\). ◻
In particular, these error estimates imply that the multilevel operator \(A_L\) is bounded.
Corollary 1. With the notation and under the assumptions of 2, the multilevel operator \(A_L\) is bounded with norm satisfying \[\|A_L\|_{H^{\sigma}(\Omega) \to W^{\tau}_q(\Omega)} \leq C \left(\alpha^L + 1 \right).\]
Proof. We see that \(\| A_L(f) \|_{ W^{\tau}_q(\Omega)} \leq \| A_L(f) - f \|_{ W^{\tau}_q(\Omega)} + \| f \|_{ W^{\tau}_q(\Omega)}.\) Using the error result ?? and the Sobolev embedding theorem in the form \(\| f \|_{W^{\tau}_q(\Omega)} \leq C \| f \|_{H^{\sigma} (\Omega)}\) completes the proof. ◻
Throughout this section, we will assume that the data sets \(X_\ell\) are nested, i.e. that we have \(X_{\ell}\subseteq X_{\ell+1}\). We will also assume that the points are ordered in such a way that, for \(1\le \ell\le L\), the first \(N_\ell\) points of \(X_L\) belong to \(X_\ell\), i.e. \[X_{\ell} = \{\boldsymbol{x}_1,\ldots,\boldsymbol{x}_{N_\ell}\},\] allowing us to drop upper indices to indicate the current level. This also means that \(X_\ell\setminus X_{\ell-1} = \{\boldsymbol{x}_{N_{\ell-1}+1},\ldots, \boldsymbol{x}_{N_\ell}\}\) are those points which are newly introduced on level \(\ell\). Moreover, if we set \(N_0=0\) then for every \(1\le i\le N_L\) there is an \(1\le \ell \le L\) such that \(N_{\ell-1}+1 \le i \le N_{\ell}\), meaning that the point \(\boldsymbol{x}_i\) belongs to \(X_{\ell}\setminus X_{\ell-1}\), i.e. it is first introduced into the data set at level \(\ell\).
In this situation, we now want to derive new ways of representing the multilevel operator \(A_L\), i.e. we are looking for good bases of \(\widetilde{V}_L=A_L(C(\Omega))\). First ideas in this direction have been introduced in [30] and [27], which we will refine here.
Though we have introduced the multilevel operator \(A_L\) as an operator \(A_L:C(\Omega)\to V_L\), it actually uses only \(f|X_L\) to compute \(A_L(f)\) for any \(f\in C(\Omega)\) if the data sets \(X_\ell\) are nested. For every vector \(\boldsymbol{y}\in\mathbb{R}^{N_L}\) we can obviously find a function \(f\in C(\Omega)\) with \(f|X_L=\boldsymbol{y}\), which means that we can, with a slight abuse of notation, also consider the multilevel operator as a mapping \(A_L:\mathbb{R}^{N_L}\to V_L\). This leads to the following first definition of a basis for \(\widetilde{V}_L=A_L(C(\Omega))= A_L(\mathbb{R}^{N_L})\).
Definition 3. Denote the \(i\)-th unit vector in \(\mathbb{R}^{N_L}\) by \(\boldsymbol{e}_i^{(L)}\). Then, the multilevel cardinal basis \(\mathcal{B}=\{b_1,\ldots, b_{N_L}\}\) for the approximation space \(\widetilde{V}_L\) is defined by \[b_i = A_L\left(\boldsymbol{e}_i^{(L)}\right), \qquad 1\le i\le N_L.\]
Note that we also have \(b_i=A_L(\chi_i^{(L)})\), using the previously introduced Lagrange functions, as we have \(\chi_i^{(L)}|X_L=\boldsymbol{e}_i^{(L)}\). Writing \(f|X_L\) as \(f|X_L = \sum_{i=1}^{N_L} f(\boldsymbol{x}_i) \boldsymbol{e}_i^{(L)}\) and the linearity of \(A_L\) immediately yield \[A_L f = \sum_{j=1}^{N_L} f(\boldsymbol{x}_i) A_L\left(\boldsymbol{e}_i^{(L)}\right) = \sum_{j=1}^{N_L} f(\boldsymbol{x}_i) b_i,\] showing that \(\mathcal{B}\) spans \(\widetilde{V}_L\). It is also indeed a cardinal or Lagrange basis.
Theorem 3. The multilevel cardinal basis is indeed a basis of \(\widetilde{V}_L = A_L(C(\Omega))\), showing particularly \(\dim \widetilde{V}_L=N_L\). Moreover, the basis functions are cardinal in the sense that \(b_i(\boldsymbol{x}_j) = \delta_{ij}\) for \(1\le i,j\le N_L\).
Proof. The cardinal condition follows from the fact that \(A_Lf\) interpolates \(f\) on \(X_L\). Here, this means \(b_i(\boldsymbol{x}_j) = A_L(\chi_i^{(L)})(\boldsymbol{x}_j) = \chi_i^{(L)}(\boldsymbol{x}_j) = \delta_{ij}.\) This condition also guarantees that the functions \(b_i\) are linearly independent. This means \(\dim \operatorname{span}(\mathcal{B}) = N_L\) and the above considerations have shown \(\widetilde{V}_L\subseteq \operatorname{span}(\mathcal{B})\). Finally, any \(f\in \operatorname{span}(\mathcal{B})\) can be written as \(f=\sum \alpha_i b_i\). This, however, means \(f=A_L(\boldsymbol{\alpha})\) and thus \(f\in \widetilde{V}_L\) showing \(\widetilde{V}_L = \operatorname{span}(\mathcal{B})\). ◻
The above result shows that \(\widetilde{V}_L\) is indeed a proper subspace of \(V_L\), as the latter has dimension \(N_1+\cdots + N_L\), while the former has only dimension \(N_L\).
To better understand this multilevel cardinal basis, we will now derive several other representations. To this end, we will use an alternative way of expressing the multilevel operator, see [27], [30]. In what follows, we will always assume that subsets \(\mathfrak{u}= \{u_1,\ldots, u_{|u|}\} \subseteq \{1,\ldots, L\}\) are ordered as \(u_1 < u_2 < \cdots < u_{|\mathfrak{u}|}\).
Theorem 4. The multilevel operator \(A_L\) can be represented as \[A_L = \sum_{\substack{\mathfrak{u}\subseteq \{1, \ldots, L \} \\ 1 \leq |\mathfrak{u}| \leq L}} (-1)^{|\mathfrak{u}|+1} \mathcal{I}_{\mathfrak{u}} = \sum_{j = 1}^{L} (-1)^{j+1} \sum_{\substack{\mathfrak{u} \subseteq \{1, \dots, L \}\\ |\mathfrak{u}| = j}} \mathcal{I}_{\mathfrak{u}},\] where \(\mathcal{I}_{\mathfrak{u}}\) is the identity if \(\mathfrak{u}=\emptyset\) or otherwise \(\mathcal{I}_{\mathfrak{u}}:C(\Omega)\to W_{|\mathfrak{u}|}\) is the combined operator \[\label{eq:CombinedInterpolationOperator} \mathcal{I}_{\mathfrak{u}} = I_{X_{\mathfrak{u}_{|\mathfrak{u}|}}, \Phi_{\mathfrak{u}_{|\mathfrak{u}|}}} I_{X_{\mathfrak{u}_{|\mathfrak{u}|-1}}, \Phi_{\mathfrak{u}_{|\mathfrak{u}|-1}}} \cdots I_{X_{u_1}, \Phi_{u_1}}.\qquad{(2)}\]
The ordering in (?? ) is important, as the interpolation operators do not commute.
To apply this representation, we note that for any index \(N_{m-1} +1 \le i \le N_m\) with \(1\le m\le L\) and any \(1\le u_1\le m-1\) we have \(\chi_i^{(L)}|X_{u_1}=0\), showing particularly \(I_{X_{u_1},\Phi_{u_1}} (\chi_i^{(L)}) = 0.\) This immediately yields the following result, which is essentially a reformulation of [27] and provides a first alternative representation of \(\mathcal{B}=\{b_i : 1 \leq i \leq N_L\}\).
Corollary 2. Assume that the index \(i\) satisfies \(N_{m-1}+1\le i\le N_m\) for some \(1\le m \le L\). Then, the corresponding basis function has the representation \[\label{eq:DefBi} b_i = \sum_{\emptyset \neq \mathfrak{u}\subseteq \{m, \dots, L\}} (-1)^{|\mathfrak{u}| +1} \mathcal{I}_{\mathfrak{u}} \left( \chi^{(L)}_{i} \right), \quad 1 \leq i \leq N_L,\qquad{(3)}\] with \(\mathcal{I}_{\mathfrak{u}}\) as in ?? .
Proof. This immediately follows from \[b_i = A_L\left(\boldsymbol{e}_i^{(L)}\right) = \sum_{\substack{\mathfrak{u}\subseteq \{1, \ldots, L \} \\ 1 \leq |\mathfrak{u}| \leq L}} (-1)^{|\mathfrak{u}|+1} \mathcal{I}_{\mathfrak{u}}\left(\boldsymbol{e}_i^{(L)}\right)\] and the fact that for all \(\mathfrak{u}\) with \(u_1\le m-1\) we have \(\mathcal{I}_{\mathfrak{u}}(\chi_i^{(L)}) = 0\). ◻
While this is already an improvement, it still sums over too many sets \(\mathfrak{u}\). As before, we can alternatively use \(\mathcal{I}_{\mathfrak{u}}(\boldsymbol{e}_i^{(L)})\) instead of \(\mathcal{I}_{\mathfrak{u}}(\chi_i^{(L)})\) in ?? .
For a more expressive representation of the function \(b_i\), we need the following lemma concerning interpolation of Lagrange functions. The idea of its proof has already been laid out above.
Lemma 2. Let \(X_\ell \subseteq X_L \subseteq \Omega\) be two sets of data sites with associated cardinal functions \(\{\chi^{(\ell)}_{i} : 1 \leq i \leq N_\ell\} \subseteq W_\ell\) and \(\{\chi^{(L)}_{i} : 1 \leq i \leq N_L\} \subseteq W_L\), respectively. Then, for any \(1\le i\le N_L\) we have \[I_{X_\ell, \Phi_\ell} \left(\chi^{(L)}_{i} \right) = \begin{cases} \chi^{(\ell)}_{i},& \quad 1 \leq i \leq N_\ell, \\ 0,& \quad N_{\ell}+1 \leq i \leq N_L. \end{cases}\]
Proof. The Lagrange property of \(\chi^{(L)}_{i}\) implies that \(\chi^{(L)}_{i}|{X_\ell} = \boldsymbol{e}_i \in \mathbb{R}^{N_\ell}\) if \(i \leq N_\ell\), i.e. if \(\boldsymbol{x}_{i} \in X_\ell\) and \(\chi^{(L)}_{i}|{X_\ell} = \boldsymbol{0} \in \mathbb{R}^{N_\ell}\) if \(i > N_\ell\), i.e. if \(\boldsymbol{x}_{i} \notin X_\ell\). The claim then follows from the uniqueness of the interpolant. ◻
In essence, 2 states that interpolating a fine-level Lagrange function on a coarse level yields the corresponding coarse-level Lagrange function if its anchor point lies in the coarse set of data sites. Otherwise, the result is zero.
We need another idea, which is very similar to the one that we used a few times before. If \(f|X_\ell = 0\) for some \(1\le \ell\le L\) then the first \(\ell\) steps of the multilevel method produce the zero function as the approximation and the actual approximation starts with step \(\ell+1\). We will use this and hence introduce the notation \(A_{\{\ell+1,\dots,L\}}:C(\Omega) \to V_L\) for the multilevel operator starting at level \(\ell+1\). It is easy to see that the representation from Theorem 4 becomes in this situation \[\label{restrictedmultilevel} A_{\{\ell+1, \dots,L\}}(f) := \sum_{\emptyset \neq \mathfrak{u}\subseteq \{\ell +1, \dots, L\}} (-1)^ {|\mathfrak{u}|+1} \mathcal{I}_{\mathfrak{u}} (f)\tag{6}\]
With this we can prove another representation of the multilevel cardinal bases \(\mathcal{B}= \{b_i : 1 \leq i \leq N_L\}\).
Theorem 5. With the assumptions and notation from above, let \(1 \leq m \leq L\) and \(N_{m-1}+1 \leq i \leq N_m\). Then the function \(b_i\) has the representation \[\label{eq:NewRepresentationBi} b_i = \chi^{(L)}_{i} + \sum_{\ell = m} ^{L-1} \left[\chi^{(\ell)}_{i} - A_{\{\ell+1, \dots,L\}} \left(\chi^{(\ell)}_{i} \right) \right].\qquad{(4)}\]
Proof. We start with the representation of \(b_i\) from ?? in the form \[b_i = \sum_{\emptyset \neq \mathfrak{u}\subseteq \{m, \dots, L\}} (-1)^{|\mathfrak{u}| +1} \mathcal{I}_{\mathfrak{u}} \left( \chi^{(L)}_{i} \right).\] Next, we order the sets \(\emptyset \neq \mathfrak{u}\subseteq \{m, \dots, L \}\) according to their smallest element, which we denote by \(\ell\) and which obviously satisfies \(\ell \ge m\). Thus, we set \(\mathfrak{u}= \{\ell \} \cup \widetilde{\mathfrak{u}}\), where \(\widetilde{\mathfrak{u}} \subseteq \{\ell +1, \dots, L \}\). In particular, we have for \(\ell = L\) that \(\mathfrak{u}= \{L\}\), \(\widetilde{\mathfrak{u}} = \emptyset\) and \[\sum_{\mathfrak{u}= \{L\}} (-1)^{| \mathfrak{u}| +1} \mathcal{I}_{\mathfrak{u}} (\chi^{(L)}_{i}) = I_{X_L,\Phi_L} (\chi^{(L)}_{i}) = \chi^{(L)}_{i},\] since \(\chi^{(L)}_{i} \in W_L\) and the uniqueness of the interpolant. This leads to \[\begin{align} b_i &= \chi^{(L)}_{i} + \sum_{\ell = m}^{L-1} \sum_{\widetilde{\mathfrak{u}} \subseteq \{\ell +1,\dots, L \}} (-1)^{| \widetilde{\mathfrak{u}}|} \mathcal{I}_{\widetilde{\mathfrak{u}}} \left(I_{X_\ell,\Phi_\ell}(\chi^{(L)}_{i})\right) \\ &= \chi^{(L)}_{i} + \sum_{\ell = m}^{L-1} \sum_{\widetilde{\mathfrak{u}} \subseteq \{\ell+1,\dots, L \}} (-1)^{| \widetilde{\mathfrak{u}}|} \mathcal{I}_{\widetilde{\mathfrak{u}}} \left(I_{X_\ell,\Phi_\ell}(\chi^{(\ell)}_{i}) \right), \end{align}\] where we used 2 to go from \(I_{X_\ell,\Phi_\ell}(\chi^{(L)}_{i})\) to \(I_{X_\ell,\Phi_\ell}( \chi^{(\ell)}_{i})\) in the last step, which we were allowed to do as we have \(i\le N_m\le N_\ell\). To arrive at the claim, we see that, for fixed \(m \leq \ell \leq L-1\), \[\begin{align} \sum_{\widetilde{\mathfrak{u}} \subseteq \{\ell +1,\dots, L \}} (-1)^{| \widetilde{\mathfrak{u}}|} \mathcal{I}_{\widetilde{\mathfrak{u}}}\left( I_{X_\ell,\Phi_\ell} (\chi^{(\ell)}_{i})\right) &= I_{X_\ell,\Phi_\ell} (\chi^{(\ell)}_{i}) - A_{\{\ell+1, \dots, L\}} \left(I_{X_\ell,\Phi_\ell} (\chi^{(\ell)}_{i}) \right) \\ &= \chi^{(\ell)}_{i} - A_{\{\ell+1, \dots, L\}} (\chi^{(\ell)}_{i}), \end{align}\] using also (6 ). ◻
As it turns out, we can even prove a relation of the image spaces \(\widetilde{V}_\ell = A_\ell(C(\Omega))\) of multilevel operators on different levels.
Corollary 3. For all \(1 \leq \ell \leq L\) we have \(A_L(\widetilde{V}_\ell) = \widetilde{V}_\ell\). To be more precise, for any \(f\in C(\Omega)\), we have \(A_L(A_\ell(f)) = A_\ell(f)\). This means in particular \(A_L^2 = A_L\) and the nestedness of the the approximation spaces, i.e. \[\widetilde{V}_1 \subseteq \widetilde{V}_2 \subseteq \cdots \subseteq \widetilde{V}_L.\]
Proof. Let \(f\in C(\Omega)\) be given. We prove \(A_L(A_\ell(f)) = A_\ell(f)\) by induction on \(m:=L-\ell\in\mathbb{N}_0\) for \(0\le m\le L\). For \(m=0\) we have \(\ell=L\). The multilevel method interpolates on the finest level, i.e. \((A_L f)|X_L=f|X_L\) and as it only depends on the values at \(X_L\) we thus have \(A_L(A_L(f))=A_L(f)\).
Next, assume that the statement is true for any \(m\le L-1\). To show it for \(m+1\le L\), we define \(\ell\) via \(m+1=L-\ell\). The recursive definition of the multilevel method and the induction hypothesis for \(m=L-1-\ell\) yield \[\begin{align} A_L(A_\ell(f)) &= A_{L-1}(A_\ell(f)) + \mathcal{I}_{X_L,\Phi_L}(A_\ell(f)-A_{L-1}(A_\ell(f)))\\ &= A_\ell(f) + \mathcal{I}_{X_L,\Phi_L} (A_\ell(f)-A_\ell(f))\\ & = A_\ell(f), \end{align}\] completing the induction. The statement \(A_L^2=A_L\) immediately follows from this. Moreover, we have \(\widetilde{V}_\ell = A_\ell(C(\Omega)) \subseteq C(\Omega)\) and thus \(\widetilde{V}_\ell = A_{\ell+1}(\widetilde{V}_\ell) \subseteq A_{\ell+1}(C(\Omega)) = \widetilde{V}_{\ell+1}\), proving the nestedness of the approximation spaces. ◻
This nestedness and the finite dimensionality of the approximation spaces \(\widetilde{V}_\ell\) allows us to decompose them in the form \[\label{decomp2} \widetilde{V}_{\ell} = \widetilde{V}_{\ell-1} + \widetilde{W}_{\ell}\tag{7}\] and we will now prove that this sum is again a direct sum and introduce another basis for \(\widetilde{V}_L\), which also gives another representation of the multilevel operator.
Definition 4. A Newton-type basis for the approximation space \(\widetilde{V}_L\) is defined by \[\widetilde{\mathcal{B}} = \left\{ \chi_i^{(\ell)} : 1\le \ell \le L, N_{\ell-1}+1\le i\le N_\ell\right\}.\]
The next theorem will show that this is indeed a basis for \(\widetilde{V}_L\). It then follows that in the above decomposition (7 ), the detail space \(\widetilde{W}_\ell\) has the basis \(\{\chi_i^{(\ell)} : N_{\ell-1}+1 \le i\le N_{\ell}\}\). This new bases demonstrates more accurately the inductive structure of the multilevel method. Again, we see that the additional basis elements \(\{\chi^{(\ell)}_{i} : N_{\ell-1} + 1 \leq i \leq N_\ell \}\) vanish on \(X_{\ell-1}\), showing a Newton-type behavior.
Theorem 6. With \(N_0:=0\) and \(A_0(f):= 0\), the multilevel approximation at level \(L\) to \(f \in C(\Omega)\) can be written as \[\label{rep2} A_L(f) = \sum_{\ell=1}^{L} \sum_{i=N_{\ell-1}+1}^{N_\ell} \left[f(\boldsymbol{x}_i)-A_{\ell-1} f(\boldsymbol{x}_i)\right] \chi_i^{(\ell)}.\qquad{(5)}\] Consequently, the approximation space \(\widetilde{V}_L=A_L(C(\Omega))\) is recursively given by \[\widetilde{V}_{L} = \begin{cases} \operatorname{span}\left\{\chi_i^{(1)} : 1\le i\le N_1\right \} &ifL=1, \\ \widetilde{V}_{L-1} \oplus \operatorname{span} \left\{\chi_i^{(L)} : N_{L-1}+1 \le i\le N_L\right\} &ifL\ge 2. \end{cases}\]
Proof. We prove the representation (?? ) again by induction on \(L\). For \(L =1\) we have \(A_1(f) = I_{X_1,\Phi_1}f\) and hence the statement.
For \(L \geq 2\) we again use the definition of the multilevel operator to conclude \[\begin{align} A_L(f) &= A_{L-1}(f) + I_{X_L, \Phi_L}(f - A_{L-1}(f)) \nonumber\\ &= A_{L-1}(f) + \sum_{i=1}^{N_L} \left( f(\boldsymbol{x}_i) - A_{L-1}(f)(\boldsymbol{x}_i) \right) \chi^{(L)}_{i}\nonumber\\ &= A_{L-1}(f) + \sum_{i=N_{L-1}+1}^{N_L} \left(f(\boldsymbol{x}_i)- A_{L-1}(f)(\boldsymbol{x}_i) \right) \chi^{(L)}_{i},\label{eq:RecursiveRepresentationAj} \end{align}\tag{8}\] where we used that \(A_{L-1}f|X_{L-1}= f|X_{L-1}\) in the last step. The induction hypothesis now leads to the representation (?? ). Moreover, this shows that \(\widetilde{\mathcal{B}}\) spans \(\widetilde{V}_L\) and a comparison of the dimensions shows that \(\widetilde{\mathcal{B}}\) is even a basis of \(\widetilde{V}_L\). This also means that the stated sum of subspaces is direct. ◻
Additionally, 8 allows us to derive a recursive representation of \(b_i\), where the recursion is done over the levels. To emphasize the level-dependence, we will write \(b_i^{(\ell)}\) for the \(i\)-th basis function of the multilevel cardinal basis \(\mathcal{B}_\ell := \{b_i^{(\ell)} : 1\le i\le N_\ell\}\) of \(\widetilde{V}_{\ell}\).
Theorem 7. The basis functions \(b_i^{(L)}\in\mathcal{B}_L\) of the multilevel cardinal basis satisfy \[b_i^{(L)} = \begin{cases} \displaystyle b_i^{(L-1)} - \sum_{j = N_{L-1}+1}^{N_L} b_i^{(L-1)}(\boldsymbol{x}_j) \chi^{(L)}_j & for1 \leq i \leq N_{L-1},\\ \chi_i^{(L)}& forN_{L-1}+1 \leq i \leq N_L. \end{cases}\]
Proof. Let \(f \in C(\Omega)\) be fixed. We use the representation 8 of \(A_L(f)\) to derive \[\begin{align} A_L(f) &= A_{L-1}(f) + \sum_{j=N_{L-1}+1}^{N_L} \left(f(\boldsymbol{x}_j)- A_{L-1}(f)(\boldsymbol{x}_j) \right) \chi^{(L)}_{j} \\ &= \sum_{i=1}^{N_{L-1}} f(\boldsymbol{x}_i) b_i^{(L-1)} + \sum_{j=N_{L-1}+1}^{N_L} \left( f(\boldsymbol{x}_j) - \sum_{i=1}^{N_{L-1}} f(\boldsymbol{x}_i) b_i^{(L-1)}(\boldsymbol{x}_j) \right) \chi_j^{(L)} \\ &= \sum_{i=1}^{N_{L-1}} f(\boldsymbol{x}_i) \left( b_i^{(L-1)} - \sum_{j=N_{L-1}+1}^{N_L} b_i^{(L-1)}(\boldsymbol{x}_j) \chi_j^{(L)} \right) + \sum_{i=N_{L-1}+1}^{N_L} f(\boldsymbol{x}_i) \chi_i^{(L)}. \end{align}\] A comparison to \(A_L(f) = \sum_{i=1}^{N_L} f(\boldsymbol{x}_i) b_i^{(L)}\) finishes the proof. ◻
It is well-established that the cardinal functions \(\chi_i^{(\ell)}\) of \(W_\ell\) decay exponentially with growing \(\| \boldsymbol{x}- \boldsymbol{x}_i \|_2/q_\ell\). A precise formulation for nested data sets is in the next theorem, though the nestedness is not really necessary. For its proof we refer to [33].
Theorem 8. Let \(\Phi\) be a compactly supported reproducing kernel of \(H^{\sigma}(\mathbb{R}^d)\), i.e., its Fourier transform satisfies 3 . Let \(X_1\subseteq X_2\subseteq \cdots \subseteq X_L \subseteq \mathbb{R}^d\) be a family of quasi-uniform sets of data sites with fill distances \(h_1, \ldots, h_L\). For \(1 \leq \ell \leq L\), let \(\Phi_{\ell} = \Phi_{\delta_\ell}\) be defined by 1 with \(\delta_\ell = \nu h_\ell\), with \(\nu \geq 1\). Then, there are constants \(C_1 > 0\) and \(\eta > 0\) that are independent of \(\ell\), such that the bound \[| \chi_i^{(\ell)}(\boldsymbol{x}) | \leq C_1 e^{- \eta \frac{\| \boldsymbol{x}- \boldsymbol{x}_i \|_2}{q_\ell}}, \qquad \boldsymbol{x}\in\mathbb{R}^d\] holds for all \(1 \leq \ell \leq L\) and \(1 \leq i \leq N_\ell\).
In particular, the basis functions of the Newton-type basis \(\mathcal{B}_{\operatorname{full}}\) of \(V_L\) and the basis functions from the Newton-type basis \(\widetilde{\mathcal{B}}\) of \(\widetilde{V}_{L}\) enjoy this type of exponential decay.
Following ?? , the multilevel cardinal basis \(\mathcal{B}=\{b^{(L)}_i : 1\le i\le N_L\}\) of \(\widetilde{V}_L\) is closely related to the Newton-type basis \(\widetilde{\mathcal{B}}\). In particular, for \(N_{L-1}+1 \leq i \leq N_L\), the basis function \(b_i^{(L)}\) is even equal to \(\chi^{(L)}_i\). Hence, we expect to obtain a similar exponential decay for the basis functions from \(\mathcal{B}\).
Theorem 9. With the notation and assumptions of 8, fix \(L \in \mathbb{N}\). Then, there is a constant \(C=C_L>0\) such that for all \(1\le m\le L\) and \(N_{m-1}+1\le i\le N_m\) the estimate \[| b^{(L)}_i(\boldsymbol{x})| \le C_L e^{-\frac{\eta}{q_m} \| \boldsymbol{x}- \boldsymbol{x}_i \|_2}, \qquad \boldsymbol{x}\in\mathbb{R}^d,\] holds, where \(\eta>0\) is the constant from 8.
Proof. The proof is once again by induction on the level \(L\). For \(L=1\) we note that we must have \(m=1\) and \(1\le i\le N_1\). Thus, we have \(b_i^{(1)} = \chi_i^{(1)}\) and thus 8 shows the statement.
For the inductions step, we assume that the result is correct for \(L-1\) and thus \(1\le m\le L-1\). For level \(L\) we might have \(m=L\), which then again leads to \(b_{i}^{(L)} = \chi_i^{(L)}\) for \(N_{L-1}+1 \le i\le N_L\) and Theorem 8 gives the desired statement. If \(m\le L-1\), we can use the recursion \[\label{bi1} b_i^{(L)}(\boldsymbol{x}) = b_i^{L-1}(\boldsymbol{x}) - \sum_{j=N_{L-1}+1}^{N_L} b_i^{(L-1)}(\boldsymbol{x}_i) \chi_j^{(L)}(\boldsymbol{x})\tag{9}\] together with the induction assumption for bounding \(b_i^{(L-1)}\) and 8 for bounding \(\chi_{j}^{(L)}\). This yields first of all \[\begin{align} \left| b_i^{(L-1)}(\boldsymbol{x}_i) \chi_j^{(L)}(\boldsymbol{x})\right| & \le C_1C_{L-1} e^{-\frac{\eta}{q_m}\|\boldsymbol{x}_i-\boldsymbol{x}_j\|_2} e^{-\frac{\eta}{q_L}\|\boldsymbol{x}-\boldsymbol{x}_j\|_2}\\ &= C_1C_{L-1} e^{-\frac{\eta}{q_m}\left(\|\boldsymbol{x}_i-\boldsymbol{x}_j\|_2+\|\boldsymbol{x}-\boldsymbol{x}_j\|_2\right)} e^{-\|\boldsymbol{x}-\boldsymbol{x}_j\|_2\left(\frac{\eta}{q_L}-\frac{\eta}{q_m} \right)}\\ &\le C_1C_{L-1}e^{-\frac{\eta}{q_m}\|\boldsymbol{x}-\boldsymbol{x}_i\|_2} e^{-\|\boldsymbol{x}-\boldsymbol{x}_j\|_2\left(\frac{\eta}{q_L}-\frac{\eta}{q_m} \right)}. \end{align}\] Next, we turn to the sum over these terms. The above estimates immediately lead to \[\sum_{j=N_{L-1}+1}^{N_L} \left| b_i^{(L-1)}(\boldsymbol{x}_i) \chi_j^{(L)}(\boldsymbol{x})\right| \le C_1C_{L-1} e^{-\frac{\eta}{q_m}\|\boldsymbol{x}-\boldsymbol{x}_i\|_2} \sum_{j=N_{L-1}+1}^{N_L} e^{-\|\boldsymbol{x}-\boldsymbol{x}_j\|_2\left(\frac{\eta}{q_L}-\frac{\eta}{q_m} \right)}\] and we now need to show that the latter sum can be bounded by a constant. To this end, we introduce \(E_n = \{\boldsymbol{y}\in\mathbb{R}^d : nq_L \le \|\boldsymbol{x}-\boldsymbol{y}\|_2<(n+1) q_L\}\) for \(n\in \mathbb{N}_0\) and note that any \(\boldsymbol{x}_i\in X_L\) must be contained in exactly one of these \(E_n\). Moreover, a simple volume argument shows \(|X_L\cap E_0|\le 2^d\) and, as in the proof of [3] Theorem 12.3, \(|X_L\cap E_n|\le 3^d n^{d-1}\) for \(n\ge 1\). Thus, we find \[\begin{align} \sum_{j=N_{L-1}+1}^{N_L} e^{-\|\boldsymbol{x}-\boldsymbol{x}_j\|_2\left(\frac{\eta}{q_L}-\frac{\eta}{q_m}\right)} & \le \sum_{n=0}^\infty \sum_{\boldsymbol{x}_j\in E_n} e^{-nq_L \left(\frac{\eta}{q_L}-\frac{\eta}{q_m}\right)}\\ & \le 2^d + 3^d \sum_{n=1}^\infty n^{d-1} e^{-n\eta\left(1 - \frac{q_L}{q_m}\right)}\\ &\le 2^d +3^d \sum_{n=1}^\infty n^{d-1} e^{-n\eta\left(1 - \frac{q_L}{q_{L-1}}\right)} =: \widetilde{C}_L. \end{align}\] Note that our assumptions on \(q_\ell\) and \(h_\ell\) give the bound \(q_L/q_{L-1} \le \mu c_{qu}\). Thus, if \(\mu c_{qu}<1\) then we can choose the constant \(\widetilde{C}_L\) even independently of \(L\).
In any case, plugging this all into 9 , and using the induction hypothesis on the first term again yields \[\begin{align} |b_i^{(L)}(\boldsymbol{x})| &\le \left| b_i^{L-1}(\boldsymbol{x})\right| + \sum_{j=N_{L-1}+1}^{N_L} \left| b_i^{(L-1)}(\boldsymbol{x}_i) \chi_j^{(L)}(\boldsymbol{x})\right| \\ &\le C_{L-1} e^{-\frac{\eta}{q_m}\|\boldsymbol{x}-\boldsymbol{x}_i\|_2} + C_1C_{L-1} \widetilde{C}_L e^{-\frac{\eta}{q_m}\|\boldsymbol{x}-\boldsymbol{x}_i\|_2} = C_L e^{-\frac{\eta}{q_m}\|\boldsymbol{x}-\boldsymbol{x}_i\|_2} \end{align}\] with \(C_L:=C_{L-1}(1+C_1 \widetilde{C}_L)\). ◻
While the basis functions of both \(\mathcal{B}\) and \(\widetilde{\mathcal{B}}\) are now proven to have an exponential decay there is an important difference between them. For the basis functions \(\chi_i^{(\ell)}\) of \(\widetilde{\mathcal{B}}\) the constant \(C_1>0\) is independent of the level, while the constant \(C_L>0\) for the basis functions \(b_i^{(L)}\) from \(\mathcal{B}\) depends on the highest level \(L\). As this does not show up in numerical calculations, it remains an open problem whether this constant can also be chosen independently of \(L\).
An exponential decay of the basis functions immediately leads to uniform boundedness of their \(\ell_1\)-sum, as shown in the next corollary.
Corollary 4. The \(\ell_1\)-norm of the basis functions in \(\mathcal{B}\) and in \(\widetilde{\mathcal{B}}\) are uniformly bounded. To be more precise there exist a level independent constant \(\widetilde{C}_1>0\) and a level dependent constant \(\widetilde{C}_L>0\) such that \[\sum_{j=1}^N |b_j^{(L)}(\boldsymbol{x})| \le \widetilde{C}_L, \qquad \sum_{\ell=1}^L \sum_{i=N_{\ell-1}+1}^{N_\ell}|\chi_i^{(\ell)}(p\boldsymbol{x})|\le \widetilde{C}_1\] for all \(\boldsymbol{x}\in \mathbb{R}^d\).
Proof. As the proof is the same for both bases, we let \(\phi_i\) be either \(\chi_i^{(\ell)}\) or \(b_i^{(L)}\) for some \(N_{\ell-1}+1\le i\le N_\ell\) and \(1\le \ell\le L\). Then, the exponential decay gives \[|\phi_i(\boldsymbol{x})| \le C e^{-\frac{\eta}{q_\ell}\|\boldsymbol{x}-\boldsymbol{x}_i\|_2} \le C e^{-\frac{\eta}{q_1}\|\boldsymbol{x}-\boldsymbol{x}_i\|_2} , \qquad \boldsymbol{x}\in\mathbb{R}^d.\] This time, letting \(E_n = \{\boldsymbol{y}\in\mathbb{R}^d : nq_1 \le \|\boldsymbol{x}-\boldsymbol{y}\|_2< (n+1)q_1\}\) yields, as in the proof of Theorem 9, \[\sum_{i=1}^{N_L} |\phi_i(\boldsymbol{x})| \le C \sum_{n=0}^\infty \sum_{\boldsymbol{x}_i\in E_n} e^{-\eta n} \le C\left(2^d + \sum_{n=0}^\infty n^{d-1}e^{-\eta n}\right),\] where the sum is clearly a level-independent constant. Hence, the overall constant is only level-dependent for the multilevel cardinal basis but not the Newton-type basis. ◻
Note that the result on \(\mathcal{B}\) means that the associated Lebesgue function, which is simply the \(\ell_1\)-norm of the basis functions in the case of a cardinal basis and which coincides with \(\|A_L\|_{L_\infty\to L_\infty}\) is bounded by \(\widetilde{C}_L\). Unfortunately, again the question arises whether this can be done with a constant that is independent of the level. Only in that case, it immediately follows that \(A_L(f)\) also converges to \(f\) point-wise for even any function \(f\in C(\overline{\Omega})\).
Similarly, it is possible to show that the basis functions of both bases are Lipschitz continuous. Following [34], this then leads to the first step in the direction of showing that \(\mathcal{B}\) is a stable basis. To be more precise, for an \(f = \sum_{j=1}^{N_L} a_jb_j^{(L)}\) let \(P_\ell f = \sum_{j={N_{\ell-1}+1}}^{N_\ell}a_jb_j^{(L)}\) and \(\boldsymbol{a}_\ell=(a_{N_{\ell-1}+1},\ldots, a_{N_\ell})^\text{T}\) then the general theory of [34] shows that \[c_1q_\ell^{d/p} \|\boldsymbol{a}_\ell\|_{\ell_p} \le \|P_\ell f\|_{L_p(\Omega)}\le c_2q_\ell^{d/p} \|\boldsymbol{a}_\ell\|_{\ell_p}.\] Unfortunately, here the constants \(c_1,c_2>0\) depend so far also on the level \(L\). Thus, we leave out the details.
From a numerical point of view, computing the cardinal bases \(\chi_{i}^{(\ell)}\) should be avoided if possible. However, using the multilevel cardinal basis has advantages when it comes to tensorizing such bases for higher dimensional problems. Hence, we shortly describe how this can be done and analyze the computational cost. To this end, we start for a fixed index \(i\) with a representation of the from \[b_i = A_L(\boldsymbol{e}^{(L)}_i) = A_L(\chi_i^{(L)}) = \sum_{\ell=1}^{L} \sum_{j = 1}^{N_{\ell}} \alpha_j^{(\ell)} \Phi_\ell( \cdot - \boldsymbol{x}_j).\] The coefficient vectors \(\boldsymbol{\alpha}^{(m)}\), \(1 \leq \ell \leq L\), are the unique solutions of the block linear system \[\begin{pmatrix} M_1 \\ B_{21} & M_{2} & & & \\ B_{31} & B_{32} & M_{3} & & \\ \vdots & \vdots & \cdots & \ddots & \\ B_{L1} & B_{L2} & \cdots & B_{L(L-1)} & M_L \end{pmatrix} \begin{pmatrix} \boldsymbol{\alpha}^{(1)} \\ \boldsymbol{\alpha}^{(2)} \\ \boldsymbol{\alpha}^{(3)} \\ \vdots \\ \boldsymbol{\alpha}^{(L)} \end{pmatrix} = \begin{pmatrix} \chi^{(L)}_i|_{X_1} \\ \chi^{(L)}_i|_{X_2} \\ \chi^{(L)}_i|_{X_3} \\ \vdots \\ \chi^{(L)}_i|_{X_L} \end{pmatrix},\]where \(M_\ell = ( \Phi_\ell(\boldsymbol{x}_{i} - \boldsymbol{x}_{j}))_{1 \leq i,j \leq N_\ell}\) is the interpolation matrix on level \(\ell\), and \(B_{\ell m} = (\Phi_m(\boldsymbol{x}_{i} - \boldsymbol{x}_{j}))_{\substack{1 \leq i \leq N_{\ell} \\ 1 \leq j \leq N_m}}\) is the transportation matrix from level \(m\) to level \(\ell > m\).
For a fixed level \(L\), the nestedness of the data sites enforces a clear recursive structure. If the anchor \(\boldsymbol{x}_i\) of \(b_i\) appears on the finest level, i.e. \(N_{L-1} +1 \leq i \leq N_L\), then we only have to solve the subsystem \(M_L \boldsymbol{\alpha}^{(L)} = \boldsymbol{e}_i.\) If \(\boldsymbol{x}_i\) appears on the second to last level, i.e., \(N_{L-2} + 1 \leq i \leq N_{L-1}\) we have to solve the subsystems on levels \(L-1\) and \(L\). In general, we only have to solve the whole block linear system only for those \(\boldsymbol{x}_i\) that appear in the first level. Thus, the complexity at level \(L\) is governed by this recursive coupling, as summarized in the following corollary.
Corollary 5. Assume that the sets \(\{X_\ell : 1\le \ell \le L\}\) are quasi-uniform. Then, the computation of all necessary coefficients of \(\{b_i^{(L)} : 1 \leq i \leq N_L\}\) can be done in \(\mathcal{O}(N_L^2 \log(N_L))\) time, if the systems can be solved in linear time. Keeping them in memory needs \[\mathcal{O} \left(\sum_{\ell=1}^{L} (N_\ell - N_{\ell-1}) \sum_{m=\ell}^L N_m \right)\] space. A single point evaluation of \(A_L(f)\) in nodal representation takes \(\mathcal{O} (N_L \log(N_L))\) time.
Proof. Assuming that the matrices are build, we have to solve \(N_L\) many multilevel problems. Each takes \(\mathcal{O}\left(N_L \log(N_L)\right)\) time. The second claim follows directly from the ideas above.
To see the cost of a single point-evaluation, we write \(A_L(f)(\boldsymbol{x})\) as \[A_L(f)(\boldsymbol{x}) = \sum_{i=1}^{N_L} f(\boldsymbol{x}_{i}))b_i^{(L)}(\boldsymbol{x}) = \sum_{i=1}^{N_L} f(\boldsymbol{x}_i) \sum_{\ell=1}^{L} \sum_{j = 1}^{N_{\ell}} \alpha_j^{(\ell)} \Phi_\ell( \boldsymbol{x}- \boldsymbol{x}_j).\] The innermost summation is actually only over those indices \(j\) with \(\|\boldsymbol{x}-\boldsymbol{x}_j\|_2\le \delta_\ell\). As the sets are quasi-uniform this can be done in constant time. Together with \(L \sim \log(N_L)\), we arrive at the claimed cost. ◻
Both the computation of the coefficients and the point evaluation is by a factor of \(N_L\) more expensive than the representation of the multilevel operator in 2. Hence, the multilevel cardinal basis \(\mathcal{B}\) should only be used when a representation with separated sampling and evaluation points is required.
We now shift our focus on non-nested point sets \(\{X_\ell :1 \leq \ell \leq L\}\), i.e., we do not assume \(X_{\ell}\subseteq X_{\ell+1}\) but we also do not assume that \(X_\ell\cap X_{m} = \emptyset\) for \(\ell\ne m\).
One way of handling this situation is to introduce a new sequence \(\{Y_\ell : 1\le \ell\le L\}\) of nested point sets \[Y_\ell:= X_1 \cup \cdots \cup X_\ell =: \{\boldsymbol{y}_1,\ldots,\boldsymbol{y}_{n_\ell}\},\] having \(n_\ell:= |Y_\ell|\) points.
Then, we can proceed as in the nested case and compute a multilevel approximation which we denote by \(A_L(Y;\cdot)\). Obviously, all previous results hold for this new sequence of data sets, i.e. we have a multilevel cardinal basis \(b_i:=A_L(Y;\boldsymbol{e}_i^{(L)})\), which is a cardinal basis with respect to \(Y_L\), and, similarly, a Newton-type basis, as well. However, these are not bases for the approximation spaces \(V_\ell\) built with the sequence \(\{X_\ell\}\) but rather for the spaces \(V_\ell(Y)\) built with the sequence \(\{Y_\ell\}\).
If we want to use the multilevel algorithm with the original sequence of data sets, we can still derive a basis of the space \(\widetilde{V}_L\) and a representation of \(A_L=A_L(X;\cdot)\). We still have that \(A_L\) is linear and uses only the data \(f|X_\ell\) for \(1\le \ell\le L\). The latter means that it uses the data \(f|Y_L\). Thus, writing once again \(f|Y_L = \sum_{i=1}^{n_L} f(\boldsymbol{y}_i)\boldsymbol{e}_i\) with the \(i\)-th unit vector \(\boldsymbol{e}_i\in\mathbb{R}^{n_L}\) we see \[A_L(X;f) = \sum_{j=1}^{n_L} f(\boldsymbol{y}_i) A_L(X;\boldsymbol{e}_i).\] Thus, we again see that \(\widetilde{V}_L = A_L(X;C(\Omega))\) is spanned by \(\mathcal{B}(Y_L) = \{b_i=A_L(X;\boldsymbol{e}_i) : 1\le i\le n_L\}\). We can not expect that the so-defined functions \(\{b_i\}\) are cardinal. Nevertheless, we can still prove some properties of \(\widetilde{V}_L\).
Theorem 10. In the case of non-nested sets \(\{X_\ell\}\) let \(Y_L=X_1\cup\cdots\cup X_L\) and let \(\mathcal{B}(Y_L) = \{b_i=A_L(X;\boldsymbol{e}_i) : 1\le i\le n_L\}\). Then, the multilevel approximation \(A_L(X;f)\) to \(f \in C(\Omega)\) is given by \[A_L(X;f) = \sum_{i=1}^N f(\boldsymbol{y}_i) b_i\] Moreover, \(\mathcal{B}(Y_L)\) forms a basis of \(\widetilde{V}_L\), showing \(\dim \widetilde{V}_L=n_L\).
Proof. We already know that \(A_L(X;f)\) can be written in the stated form and that \(\mathcal{B}(Y_L)\) spans \(\widetilde{V}_L\). Next, we will show that \(\mathcal{B}(Y_L)\) is linearly independent. To this end, assume that \(0 = \sum_{i=1}^{n_L} \alpha_i b_i(\boldsymbol{x})\) for all \(\boldsymbol{x}\in \Omega\).
Next, to simplify the notation, we choose an \(f\in C(\Omega)\) with \(f|Y_L=\boldsymbol{\alpha}\). Then, returning to the original form of the multilevel method with local approximations \(s_\ell\) and errors \(e_\ell\) we find \[0 = A_L(X;f)(\boldsymbol{x}) = \sum_{\ell=1}^L s_\ell(\boldsymbol{x}), \qquad \boldsymbol{x}\in\Omega.\] As the global approximation space is a direct sum of the local ones, \(V_L = W_1 \oplus \cdots \oplus W_L\), we immediately can conclude that \(s_\ell = 0\) for \(1\le \ell \le L\), which implies \(e_\ell = f - (s_1 + \cdots + s_\ell ) = f\) for all \(1 \leq \ell \leq L\). This, in turn, leads to \[0 = s_\ell = I_{X_\ell,\Phi_\ell}(e_{\ell-1}) = I_{X_\ell,\Phi_\ell}(f),\] showing \(f|_{X_\ell} = 0\) for all \(1 \leq j \leq L\). Consequently, this means \(\boldsymbol{\alpha}= \boldsymbol{0} \in \mathbb{R}^{n_L}\) and hence the linear independence. of \(\mathcal{B}(Y_L)\).
To see that \(\widetilde{V}_L\) is not a proper subspace of \(\operatorname{span}(\mathcal{B}(Y_L))\), we can proceed exactly as in the proof of Theorem 3. ◻
Non-nested point sets appear naturally in adaptive multilevel approximations, where points at the next level are used if the current residual is larger than a given threshold, see [35]. This has obviously a connection to a greedy selection of points, which is also a possible extension of the multilevel method introduced in 2.
In this paper we studied the structure of multilevel kernel approximation spaces. Our main contribution is a nodal representation of the multilevel operator that makes its dependence on the data explicit. This representation reveals the structure of the associated global approximation space and allows the construction of basis functions adapted to the multilevel setting.
For nested families of data sites we derived both a Lagrange basis and a Newton-type basis of the multilevel space. These bases provide a direct description of the range of the multilevel operator and enable a detailed analysis of structural properties such as localization and computational complexity. In particular, we established exponential decay of the Lagrange basis functions, which indicates that the multilevel representation inherits favorable localization properties from the underlying kernel.
We also investigated the more general case of non-nested sets. Although a cardinal basis is no longer available in this setting, the ideas developed for the nested case still allow the construction of a basis for the global approximation space and yield a characterization of the range of the multilevel operator.
The nodal representation derived in this work provides a practical description of the multilevel operator that is suitable for high-dimensional constructions, for instance when combined with sparse tensor techniques such as Smolyak-type algorithms.
Several directions for future research remain open. From a theoretical perspective it would be desirable to obtain sharper localization and stability estimates for the basis functions. Another interesting direction is the analysis of adaptive strategies which dynamically generate non-nested data sites. From a computational viewpoint, the explicit representation of the multilevel operator may enable the development of efficient algorithms for high-dimensional approximation based on sparse tensor constructions and kernel methods. Finally, it would be interesting to investigate extensions of the present framework to more general approximation settings, including anisotropic kernels and operator learning problems.