A note on sampling recovery of multivariate functions in the uniform norm

We consider the problem of the recovery of complex-valued multivariate functions on a compact domain $D\subset\mathds{R}^{d}$ in the uniform norm. We are allowed to use function samples on a suitable set of nodes $\mathbf{X}\mathrel{\mathop{\mathchar 58\relax}}=\{\mathbf{x}^{1},...,\mathbf{x}^{n}\}\subset D$. The functions are modeled as elements from some reproducing kernel Hilbert space $H(K)$ with kernel $K\colon D\times D\to\mathds{C}$. The nodes $\mathbf{X}$ are drawn independently according to a tailored probability measure depending on the spectral properties of the embedding of $H(K)$ in $L_{2}(D,\varrho_{D})$. This “random information” is used for the whole class of functions and our main focus in this paper is on worst-case errors. We additionally apply a recently introduced sub-sampling technique, see [24] to further reduce the sampling budget. Note that in one scenario the measure $\varrho_{D}$ is at our disposal and represents a certain degree of freedom. In another scenario the measure $\varrho_{D}$ is a fixed measure, namely the one where the data is coming from.

Authors distinguish two main paradigms for the recovery of functions from “samples”. Here “samples” refer to the available information. The first one uses general “linear” information about functions, i.e., evaluations given by $n$ arbitrary linear continuous functionals as a number of inner products of the Hilbert space, e.g., Fourier or wavelet coefficients. The corresponding quantity that measures the minimal worst case error is given by the Gelfand numbers/widths, see the monographs [26, 27, 29]. In case of a Hilbert source space (like in our setting) these quantities equal the so-called approximation numbers which are usually denoted by $a_{n}$ and represent the worst-case error obtained by linear recovery algorithms. Precisely, let $H(K)$ be a reproducing kernel Hilbert space (RKHS) of multivariate complex-valued functions $f(\mathbf{x})=f(x_{1},\dots,x_{d})$, $\mathbf{x}\in D\subset\mathds{R}$ with the kernel $K\colon D\times D\rightarrow\mathds{C}$. Consider the identity operator $\mathrm{Id}\colon H(K)\rightarrow F$. Further, by $\mathcal{L}(H(K),F)$ we denote the set of linear continuous operators $A\colon H(K)\rightarrow F$. The approximation numbers $a_{n}(\mathrm{Id})$ are defined as follows

If $F$ is a Hilbert space, then $a_{n}$ represents the singular numbers $\sigma_{n}$ of the corresponding embedding.

Several practical applications are modeled with “standard information”, i.e., a finite number of function values at some points. Let $f(\boldsymbol{\rm x}^{1}),\dots,f(\boldsymbol{\rm x}^{n})$ be the given data, where $\boldsymbol{\rm X}=\{\boldsymbol{\rm x}^{1},\dots,\boldsymbol{\rm x}^{n}\}$, $\boldsymbol{\rm x}^{i}\subset D$, $i=1,\dots,n$, are the nodes. The corresponding sampling numbers are defined as follows:

Note that we always have $a_{n}(\mathrm{Id})\leq g_{n}(\mathrm{Id})$.

Our goal is to recover functions using a (weighted) least squares algorithm and measure the error in $\ell_{\infty}$. Function values are taken at independently drawn random nodes in $D$ according to a tailored probability measure $\varrho$. The main feature of this approach is that the nodes are drawn once for the whole class. This is usually termed as “random information”. Note that we do not investigate Monte Carlo algorithms in this paper.

A recent breakthrough was made by D. Krieg and M. Ullrich [14] in establishing a new connection between sampling numbers and singular values for $L_{2}$-approximation. This has been extended by L. Kämmerer, T. Ullrich and T. Volkmer [9] and bounds with high probability have been obtained, see also M. Moeller and T. Ullrich [23] and M.Ullrich [38] for further refinements. In all the mentioned contributions we need a logarithmic oversampling, i.e, $n\asymp m\log m$. A new sub-sampling technique has been applied by N. Nagel, M. Schäfer and T. Ullrich in [24], where it is shown that $n=\mathcal{O}(m)$ samples are sufficient. For non-Hilbert function classes we refer to D. Krieg and M. Ullrich [13] and V.N. Temlyakov [35], where a connection to Kolmogorov widths is pointed out. There is also a relation in the $\ell_{\infty}$-setting, see [32] for details.

In the present paper we obtain worst case recovery guarantees with high probability for the weighted least squares recovery operator $\widetilde{S}_{\mathbf{X}}^{m}$ (see Algorithm 1) where the error is measured in the supremum norm of the space $\ell_{\infty}(D)$ of bounded on $D$ functions. The operator $\widetilde{S}_{\mathbf{X}}^{m}$ uses a density function $\varrho_{m}$ tailored to the spectral properties of the corresponding embedding into $L_{2}(D,\varrho_{D})$, see (3.3) and (3.4). We need to impose stronger conditions, namely we consider a compact set $D\subset\mathds{R}^{d}$, a continuous kernel $K(\mathbf{x},\mathbf{y})$ such that $\sup_{\mathbf{x}\in D}\sqrt{K(\mathbf{x},\mathbf{x})}<\infty$ and a Borel measure $\varrho_{D}$ with full support (such that Mercer’s theorem is applicable). Then we get for $m\mathrel{\mathop{\mathchar 58\relax}}=\lfloor\frac{n}{c_{1}r\log n}\rfloor$, $r>1$, the following error bound with probability larger than $1-c_{2}n^{1-r}$

with three universal constants $c_{1},c_{2},c_{3}>0$, where $N_{K,\varrho_{D}}(m)$ denotes the supremum over the spectral function (Christoffel function) which is defined in Section 2. When applying this theorem to periodic spaces with mixed smoothness $H^{s}_{{\rm mix}}(\mathds{T}^{d})$ we obtain the following result (see Theorem 7.8).

If $s>(1+\log_{2}d)/2$, $r>1$, $m=\lfloor n/(10r\log n)\rfloor$ and $\beta=2s/(1+\log_{2}d)>1$ then

with probability larger than $1-3n^{1-r}$ . Clearly, we have to determine the norm in $H^{s}_{{\rm mix}}(\mathds{T}^{d})$. In the above statement we used the $\#$-norm, see (7.7). This is a preasymptotic statement.

Note, that the sampling operator uses $n\asymp m\log m$ many sample points such that we obtain a bound for $g_{\lfloor cm\log m\rfloor}$. By a recently introduced sub-sampling technique (see [24]) we obtain for the sampling numbers another result which may be better in many cases, namely

where $c_{4},c_{5}>0$ are universal constants. We evaluate this bound for several relevant examples to obtain new bounds for sampling numbers. Here the question for reasonable assumptions arises to ensure the polynomial order $q^{*}$ of the error decay, i.e., the supremum over all $q$ such that the worst case error is of the order $n^{-q}$. We improve a result by F. Y. Kuo, G. W. Wasilkowski and H. Woźniakowski [21] and show that $q^{\rm std}_{\infty}=p-1/2$ in case $N_{K,\varrho_{D}}(k)=\mathcal{O}(k)$ and $\sigma_{n}=\mathcal{O}(n^{-p})$. We actually prove more: If there is a measure $\varrho_{D}$ such that $N_{K,\varrho_{D}}(k)=\mathcal{O}(k)$ we obtain as a consequence of (1.3) and (1.5) together with [4, Lem. 3.3] the relation

where $c_{6}>0$ is an absolute constant and $C_{\varrho_{D},K}>0$ depends on the kernel $K$ and $\varrho_{D}$, see Corollary 4.3. This improves on a recent result by L. Kämmerer, see [8, Theorem 4.11], which was stated under stronger conditions (that is, in turn, an improvement of the previous result by L. Kämmerer and T. Volkmer [10])). It shows, that in case $N_{K,\varrho_{D}}(k)=\mathcal{O}(k)$ we only loose a $\sqrt{\log n}$ or sometimes even less (assuming a polynomial decay of the $a_{n}$). In case that $N_{K,\varrho_{D}}(k)=\mathcal{O}(k^{u})$, with $u>1$, we have at least that $q^{\rm std}_{\infty}\geq p-u/2$.

As to spaces with mixed smoothness $H^{s}_{{\rm mix}}(\mathds{T}^{d})$, $s>1/2$, we know by the work of V.N. Temlyakov [34] that

where the optimal algorithm is a (deterministic) Smolyak type interpolation operator. Our algorithm yields the bound

which is only worse by a small $d$-independent $\log$-term. However, as indicated in (1.4) we also obtain preasymptotic bounds for $g_{n}(\operatorname{I}_{s})$ from (1.3) and (1.5). The open question remains whether random iid nodes are asymptotically optimal for this framework. Such a question has been recently answered positively for isotropic Sobolev spaces in [11].

Let us point out that the above general results are also applicable to the more general periodic Sobolev type spaces $H^{w}(\mathds{T}^{d})$ with a weight sequence $(w(\mathbf{k}))_{\mathbf{k}\in\mathds{Z}^{d}}$. If $(\tau_{n})_{n}$ represents the non-increasing rearrangement of the square summable weight sequence $(1/w(\mathbf{k}))_{\mathbf{k}\in\mathds{Z}^{d}}$ then we have

with three universal constants $c_{7},c_{8},c_{9}>0$. In most of the cases we do not know whether there is a deterministic method which may compete with this bound.

Notation. By $\mathds{N}$ we denote the set of natural numbers putting ${\mathds{N}}_{0}=\mathds{N}\cup\{0\}$, by $\mathds{Z}$ - the set of integer, $\mathds{R}$ - real, and $\mathds{C}$ - complex numbers. As usual, $a\in\mathds{R}$ is decomposed into $a=\lfloor a\rfloor+\{a\}$, where $0\leq\{a\}<1$ and $\lfloor a\rfloor\in\mathds{Z}$. Let further $\mathds{N}^{d}$, ${\mathds{N}}_{0}^{d}$, $\mathds{Z}^{d}$, $\mathds{R}^{d}$, $\mathds{C}^{d}$, $d\geq 2$, be the spaces of $d$-dimensional vectors $\mathbf{x}=(x_{1},\dots,x_{d})$ with coordinates, respectively, from $\mathds{N}$, ${\mathds{N}}_{0}$, $\mathds{Z}$, $\mathds{R}$, and $\mathds{C}$. The notation $\mathbf{x}\leq\mathbf{y}$, where $\mathbf{x},\mathbf{y}\in\mathds{C}^{d}$, stands for $x_{i}\leq y_{i}$, $i=1,\dots,d$ (relations $\geq$, $<$, $>$ and $=$ are defined analogically), $\overline{\mathbf{x}}$ is complex conjugate to $\mathbf{x}$, $\mathbf{x}^{\top}$ is transpose to $\mathbf{x}$. For $\mathbf{k}\in\mathds{Z}^{d}$ and $\mathbf{x}\in\mathds{C}^{d}$ we write $\mathbf{x}^{\mathbf{k}}\mathrel{\mathop{\mathchar 58\relax}}=\prod_{i=1}^{d}x_{i}^{k_{i}}$ assuming that $0^{0}\mathrel{\mathop{\mathchar 58\relax}}=1$. Let also $|\mathbf{k}|=(|k_{1}|,\dots,|k_{d}|)$, $|\mathbf{k}|_{0}$ be the number of nonzero components (sparsity) of $\mathbf{k}$. The symbols $(\mathbf{x},\mathbf{y})$ or $\mathbf{x}\cdot\mathbf{y}$ mean the Euclidean scalar product in $\mathds{R}^{d}$ or $\mathds{C}^{d}$. We use the notation $\log$ for the natural logarithm. By $\mathds{C}^{m\times n}$ we denote the set of all $m\times n$ matrices $\boldsymbol{\rm L}$ with complex entries, where, in what follows, $\|\boldsymbol{\rm L}\|$ or $\|\boldsymbol{\rm L}\|_{2\to 2}$ is their spectral norm (the largest singular value), $\boldsymbol{\rm L}^{*}$ is adjoint matrix to $\boldsymbol{\rm L}$, i.e., the matrix obtained by taking transpose and then taking complex conjugate of each entry of $\boldsymbol{{\rm L}}$. For an operator $A\in\mathcal{L}(X,Y)$ between two Banach spaces we use the symbol $\|A\|$ for its operator norm, i.e.,

If $X,Y$ are Hilbert spaces, by $A^{*}$ we denote adjoint to $A$ operator $A^{*}\colon Y\rightarrow X$, such that for $f\in X$, $g\in Y$ the following equality for the corresponding scalar products holds: $(Af,g)_{Y}=(f,A^{*}g)_{X}$. The notation $A^{-1}$ stands for inverse operator. For two sequences $(a_{n})_{n=1}^{\infty},(b_{n})_{n=1}^{\infty}\subset\mathds{R}$ we write $a_{n}\lesssim b_{n}$ if there exists a constant $c>0$, such that $a_{n}\leq cb_{n}$ for all $n$. We will write $a_{n}\asymp b_{n}$ if $a_{n}\lesssim b_{n}$ and $b_{n}\lesssim a_{n}$. If the constant $c$ depends on the dimension $d$ and smoothness $s$, we indicate it by $\lesssim_{s,d}$ and $\asymp_{s,d}$.

We will work in the framework of reproducing kernel Hilbert spaces. The relevant theoretical background can be found in [1, Chapt. 1] and [3, Chapt. 4].

Let $D\subset\mathds{R}^{d}$ be an arbitrary subset and $\varrho_{D}$ be a measure on $D$. By $L_{2}(D,\varrho_{D})$ denote the space of complex-valued square-integrable functions with respect to $\varrho_{D}$, where

and, respectively,

With $\ell_{\infty}(D)$ we denote the space of bounded on $D$ functions, for which

Note, that one do not need the measure $\varrho_{D}$ for this embedding. In fact, here we mean “boundedness” in the strong sense (in contrast to essential boundedness w.r.t. the measure $\varrho_{D}$).

We further consider a reproducing kernel Hilbert space $H(K)$ with a Hermitian positive definite kernel $K(\mathbf{x},\mathbf{y})$ on $D\times D$. The crucial property is the identity

for all $\mathbf{x}\in D$. It ensures that point evaluations are continuous functionals on $H(K)$. We will use the notation from [3, Chapt. 4]. In the framework of this paper, the boundedness of the kernel $K$ is assumed, i.e.,

The condition (2.1) implies that $H(K)$ is continuously embedded into $\ell_{\infty}(D)$, i.e.,

The embedding operator

is Hilbert-Schmidt under the finite trace condition of the kernel

see [7], [33, Lemma 2.3], which is less strong than (2.1). In view of considering in this paper a recovery of functions in the uniform norm, we assume that the kernel $K$ is bounded from now. We additionally assume that $H(K)$ is at least infinite dimensional.

Consider the embeddings $\mathrm{Id}$ from (2.3) and $W_{\varrho_{D}}=\mathrm{Id}^{*}\circ\mathrm{Id}\colon H(K)\rightarrow H(K)$. Due to the compactness of $\mathrm{Id}$, the operator $W_{\varrho_{D}}$ provides an (at most) countable system of strictly positive eigenvalues. Let $(\lambda_{n})_{n=1}^{\infty}$ be a non-increasing rearrangement of these eigenvalues, $(\sigma_{n})_{n=1}^{\infty}$ be the sequence of singular numbers of $\mathrm{Id}$, i.e., $\sigma_{k}\mathrel{\mathop{\mathchar 58\relax}}=\sqrt{\lambda_{k}}$, $k=1,2,\dots$. Further, with $(e^{*}_{n}(\mathbf{x}))_{n=1}^{\infty}\subset H(K)$ denote the system of right eigenfunctions (that correspond to the ordered system $(\lambda_{n})_{n=1}^{\infty}$), and with $(\eta_{n}(\mathbf{x}))_{n=1}^{\infty}\subset L_{2}(D,\varrho_{D})$ the system consisting of $\eta_{k}(\mathbf{x})\mathrel{\mathop{\mathchar 58\relax}}=\sigma_{k}^{-1}e^{*}_{k}(\mathbf{x})$, $k=1,2,\dots$, which both represent orthonormal systems (ONS) in the respective spaces.

We would like to emphasize that the embedding (2.3) is not necessarily injective. As discussed in [3, p. 127], $\mathrm{Id}$ is in general not injective as it maps a function to an equivalence class. As a consequence, the system of eigenvectors $(e^{*}_{n}(\mathbf{x}))_{n=1}^{\infty}$ may not be a basis in $H(K)$ (note that $H(K)$ may not even be separable). However, there are conditions which ensure that ONS $(e^{*}_{n}(\mathbf{x}))_{n=1}^{\infty}$ is an orthonormal basis (ONB) in $H(K)$, see [3, Chapt. 4.5], which is related to Mercer’s theorem, see [3, Theorem 4.49]. Indeed, if we additionally assume that the kernel $K$ is bounded and continuous on $D\times D$ (for a compact domain $D\subset\mathds{R}^{d}$), then $H(K)$ is separable and consists of continuous functions, see [1, Theorems 16, 17]. If we finally assume that the measure $\varrho_{D}$ is a Borel measure with full support then $(e^{*}_{n}(\mathbf{x}))_{n=1}^{\infty}$ is a complete orthonormal system, i.e., ONB in $H(K)$.

In this case we have the pointwise identity (see, e.g., [1, Cor. 4])

In what follows, we assume that the kernel $K$ is continuous on a compact domain $D$ (the so-called Mercer kernel). Such kernels satisfy all the indicated above conditions (see [3, Theorem 4.49]), and we even have (2.5) with absolute and uniform convergence on $D\times D$. Taking into account the spectral theorem [3, Chapt. 4.5], we have that the system $(e^{*}_{n}(\mathbf{x}))_{n=1}^{\infty}$ of eigenfunctions of the non-negative compact self-adjoint operator $W_{\varrho_{D}}$ is an orthonormal basis in $H(K)$. Note, that $(\eta_{n}(\mathbf{x}))_{n=1}^{\infty}$ is also ONB in $L_{2}(D,\varrho_{D})$.

Hence, for every $f\in H(K)$ with a Mercer kernel $K$, it holds

where $(f,e^{*}_{k})_{H(K)}$, $k=1,2,\dots$, are the Fourier coefficients of $f$ with respect to the system $(e^{*}_{n}(\mathbf{x}))_{n=1}^{\infty}$. Let us further denote

the projection onto the space $\operatorname{span}\{e^{*}_{1}(\cdot),...,e^{*}_{m}(\cdot)\}$. Due to the Parseval’s identity, the norm of the space $H(K)$ takes the following form

Note, that for $m\in\mathds{N}$ the function

is often called “Christoffel function” in literature (see [6] and references therein). For $V_{m}=\operatorname{span}\{\eta_{1}(\cdot),...,\eta_{m-1}(\cdot)\}$, we define

We may use the quantity $N_{K,\varrho_{D}}(m)$ to bound the operator norm $\|\mathrm{Id}-P_{m-1}\|_{K,\infty}$ (We use the abbreviation $\|A\|_{K,\infty}\mathrel{\mathop{\mathchar 58\relax}}=\|A\|_{H(K)\to\ell_{\infty}(D)}$). It holds

This implies

Using, that $e_{k}^{*}(\mathbf{x})=\sigma_{k}\eta_{k}^{*}$, where the sequence $(\sigma_{n})_{n=1}^{\infty}$ is non-increasing, i.e., for all $2^{l}\leq k<2^{l+1}$ it holds $\sigma_{k}^{2}\leq\frac{1}{2^{l-1}}\sum\limits_{2^{l-1}\leq j<2^{l}}\sigma_{j}^{2}$, we obtain by a straight-forward calculation

Hence, in view of the relations $2^{l+1}\leq 4j$ and $1/2^{l-1}<2/j$, that are true for all $2^{l-1}\leq j<2^{l}$, we get

Note, that the system $(\eta_{k}(\mathbf{x}))_{k=1}^{\infty}$ is not necessarily a uniformly $\ell_{\infty}$-bounded ONS in $L_{2}(D,\varrho_{D})$, nevertheless the equality (2.5) is true. For systems $(\eta_{k}(\mathbf{x}))_{k=1}^{\infty}$, where for all $k\in\mathds{N}$

we have

Let further $\boldsymbol{\rm X}=\{\boldsymbol{\rm x}^{1},\dots,\boldsymbol{\rm x}^{n}\}$, $\boldsymbol{\rm x}^{i}=({\rm x}_{1}^{i},\dots,{\rm x}_{d}^{i})\subset D$, $i=1,\dots,n$, be the drawn independently nodes according to a Borel measure $\varrho_{D}$ on $D$ that has full support. In [9] a recovery operator $S^{m}_{\boldsymbol{\rm X}}$ was constructed which computes the best least squares fit $S^{m}_{\boldsymbol{\rm X}}f$ to the given data $\mathbf{f}(\boldsymbol{\rm X})\mathrel{\mathop{\mathchar 58\relax}}=(f(\boldsymbol{\rm x}^{1}),\dots,f(\boldsymbol{\rm x}^{n}))$ from the finite-dimensional space spanned by $\eta_{k}(\mathbf{x})$, $k=1,\dots,m-1$.

Namely, let $\boldsymbol{\rm f}\mathrel{\mathop{\mathchar 58\relax}}=(f(\boldsymbol{\rm x}^{1}),\dots,f(\boldsymbol{\rm x}^{n}))^{\top}$, $\boldsymbol{\rm c}\mathrel{\mathop{\mathchar 58\relax}}=(c_{1},\dots,c_{m-1})^{\top}$, $(\eta_{k}(\mathbf{x}))_{k=1}^{\infty}=(\sigma_{k}^{-1}e^{*}_{k}(\mathbf{x})))_{k=1}^{\infty}$ and

We should solve the over-determined linear system ${\boldsymbol{\rm L}}_{m}\cdot\boldsymbol{\rm c}=\boldsymbol{\rm f}$ via least squares (e.g. directly or via the LSQR algorithm [31]), i.e., compute

The output $\boldsymbol{\rm c}\in\mathds{C}^{m-1}$ gives the coefficients of the approximant

Further, let us consider a weighted least squares algorithm (see below) using the following weight/density function, which was introduced by D. Krieg and M. Ullrich in [14]. We set

This density first appeared in [14]. We will also use the below algorithm with the simpler density function

in case $\varrho_{D}$ is a probability measure (see Remark 3.4). This density function has been implicitly used by V.N. Temlyakov in [35] and by Cohen, Migliorati [5]. In many relevant cases both densities yield similar bounds. However, as we point out in Remark 7.9 it could make a big difference even in the univariate situation.