Optimal sampling for least-squares approximation

After a short literature review (§2), this article commences with a formulation and review of (weighted) least-squares approximation (§3). We then discuss multivariate polynomial approximation (§4), this being one of the main motivating examples for this work. The next two sections contain the core developments of this article. We describe the theory of least-squares approximation with random sampling and introduce the so-called Christoffel function, which plays a key role in its analysis (§5). We then show that sampling from a measure proportional to this function leads to provably near-optimal sampling (§6). The power of such a result lies in its generality: this strategy is near optimal for any linear approximation space. Next, we consider the matter of how much can be gained through this approach in comparison to Monte Carlo sampling, i.e., i.i.d. random sampling from some underlying probability measure (§7). Monte Carlo sampling is ubiquitous in applications, especially high-dimensional approximation tasks. Yet, as we discuss, sample complexity bounds for this naïve sampling strategy can be arbitrarily bad. Once more, we see the Christoffel function plays a key role in analyzing the sample complexity. Having done this, we then conclude this part of the article by discussing series of further topics (§8). In particular, we describe very recent advances of optimal (as opposed to near-optimal) sampling and its connections to sampling numbers in information-based complexity and the study of sampling discretizations in approximation theory. We also discuss connections to matrix sketching via leverage score sampling, as well as various practical considerations.

The majority of this article considers linear approximation spaces, i.e., finite-dimensional subspaces of functions. However, modern applications increasingly make use of nonlinear spaces. Moreover, in many applications the object to recover may not be a scalar-valued function, and samples may not be simple pointwise evaluations. We conclude this article by describing a recent framework for optimal sampling with general linear samples and nonlinear approximation spaces (§9). We discuss how many of the key ideas seen in linear spaces, such as Christoffel functions, naturally extend to this general setting. Finally, we end with some concluding thoughts (§10).

In this article, we focus on the foundational techniques and theory. After a brief review in the next section, we largely omit applications, although this is currently an active area of interest. This article is intended to be accessible to nonspecialists. We build most concepts up from first principles, relying on basic knowledge only. In order to make it as self-contained as possible, proofs of most of the results shown in this work are given in an appendix.

Let $\mathcal{P}\subset L^{2}_{\varrho}(D)$ be an arbitrary $n$-dimensional subspace, where $n\leq m$, in which we seek to approximate the unknown $f$ from the measurements (3.3). We term $\mathcal{P}$ the approximation space. In this work, we consider general approximation spaces. In particular, this means that interpolation – which generally requires carefully-constructed points sets – may not be possible [47]. Instead, we consider the weighted least-squares approximation

$$\hat{f}\in{\underset{p\in\mathcal{P}}{\operatorname{argmin}}}\frac{1}{m}\sum^{m}_{i=1}w(x_{i})|y_{i}-p(x_{i})|^{2}.$$

(3.4)

Note that the loss function is well defined almost surely (for fixed $f\in L^{2}_{\varrho}(D)$ and weight function $w$ as above), since pointwise evaluations of $f$ and $w$ are well-defined almost surely.

• Remark 3.2 (The scaling factor)

  The scaling factors in (3.4) are motivated by noticing that

  $$\mathbb{E}\frac{1}{m}\sum^{m}_{i=1}w(x_{i})|g(x_{i})|^{2}=\frac{1}{m}\sum^{m}_{i=1}\int_{D}w(x)|g(x)|^{2}\,\mathrm{d}\mu_{i}(x)=\int_{D}|g(x)|^{2}\,\mathrm{d}\rho(x)={\left\|g\right\|}^{2}_{L^{2}_{\varrho}(D)},$$

  (3.5)

  where the second equality is due to (3.1). Thus, in the noiseless case, (3.4) can be considered as a empirical approximation to the continuous least-squares approximation

  $$\hat{f}={\underset{p\in\mathcal{P}}{\operatorname{argmin}}}{\left\|f-p\right\|}^{2}_{L^{2}_{\varrho}(D)},$$

  (3.6)

  i.e., the best approximation to $f$ from $\mathcal{P}$ in the $L^{2}_{\varrho}(D)$-norm. In particular, if $\mu_{1}=\cdots=\mu_{m}=\mu$, then the minimizers of (3.4) converge almost surely to the minimizer of (3.6) as $m\rightarrow\infty$ [64].

The objective of this article is to describe how to choose the measures $\mu_{1},\ldots,\mu_{m}$ to achieve the most sample-efficient approximation. We shall compare such approaches against the standard approach of Monte Carlo sampling, i.e., i.i.d. random sampling from $\varrho$. This is equivalent to setting

$$\mu_{1}=\cdots=\mu_{m}=\varrho,$$

which leads, via (3.1), to $\nu\equiv 1$. Thus, (3.4) becomes an unweighted least-squares approximation.

• Remark 3.3 (Hierarchical approximation)

  Often, rather than a fixed subspace $\mathcal{P}$, one may wish to construct a sequence of approximations in a nested collection of subspaces

  $$\mathcal{P}^{(1)}\subseteq\mathcal{P}^{(2)}\subseteq\cdots,$$

  of dimension $\dim(\mathcal{P}^{(k)})=n_{k}$. In this case, given numbers $1\leq m_{1}\leq m_{2}\leq\cdots$ satisfying $m_{k}\geq n_{k}$, $\forall k$, one aims to design nested collections of sample points

  $$\{x^{(1)}_{i}\}^{m_{1}}_{i=1}\subseteq\{x^{(2)}_{i}\}^{m_{2}}_{i=1}\subseteq\cdots.$$

  We write $\hat{f}^{(1)},\hat{f}^{(2)},\ldots$ for the ensuing (weighted) least-squares approximations, where $\hat{f}^{(k)}$ is constructed from the sample points $\{x^{(k)}_{i}\}^{m_{k}}_{i=1}$. Nested implies that samples are recycled at each iteration – a highly desirable property in the setting of limited data. We term such a procedure a hierarchical (also known as a progressive [37] or sequential [19]) approximation scheme.

Given a basis $\{\phi_{i}\}^{n}_{i=1}$ for $\mathcal{P}$, the problem (3.4) is readily reformulated as an algebraic least-squares problem for the coefficients $\hat{c}=(\hat{c}_{i})^{n}_{i=1}\in\mathbb{C}^{n}$ of $\hat{f}=\sum^{n}_{i=1}\hat{c}_{i}\phi_{i}$. This takes the form

$$\hat{c}\in{\underset{c\in\mathbb{C}^{n}}{\operatorname{argmin}}}{\left\|Ac-b\right\|}^{2}_{2},$$

(3.7)

where

$$A=\left(\sqrt{\frac{w(x_{i})}{m}}\phi_{j}(x_{i})\right)^{m,n}_{i,j=1}\in\mathbb{C}^{m\times n},\qquad b=\left(\sqrt{\frac{w(x_{i})}{m}}(f(x_{i})+e_{i})\right)^{m}_{i=1}\in\mathbb{C}^{m}.$$

(3.8)

To be precise, every minimizer $\hat{f}$ satisfying (3.4) has coefficients $\hat{c}$ that satisfy (3.7) and vice versa. Classical least-squares analysis asserts that any vector $\hat{c}$ satisfying (3.7) is also a solution of the normal equations

$$A^{*}Ac=A^{*}b,$$

(3.9)

and vice versa. Rewriting the normal equations in terms of functions also leads to the following variational form of (3.4).

$$\text{Find $\hat{f}\in\mathcal{P}$ such that }\langle\hat{f},p\rangle_{\mathsf{disc},w}=\langle f,p\rangle_{\mathsf{disc},w}+\frac{1}{m}\sum^{m}_{i=1}w(x_{i})e_{i}\overline{p(x_{i})},\ \forall p\in\mathcal{P}.$$

(3.10)

This is equivalent to (3.9) in the same sense as before. Here we wrote

$$\langle g,h\rangle_{\mathsf{disc},w}=\frac{1}{m}\sum^{m}_{i=1}w(x_{i})\overline{g(x_{i})}h(x_{i}),$$

(3.11)

for the discrete semi-inner product induced by the sample points and the weight function (whenever defined). For convenience, we shall denote the corresponding seminorm as

$${\|g\|}^{2}_{\mathsf{disc},w}=\frac{1}{m}\sum^{m}_{i=1}w(x_{i})|g(x_{i})|^{2}.$$

(3.12)

In the noiseless case $e=0$, the formulation (3.10) asserts that $\hat{f}$ is precisely the orthogonal projection of $f$ onto $\mathcal{P}$ with respect to the discrete semi-inner product (3.11). Since (3.11) is an empirical approximation to the continuous inner product $\langle\cdot,\cdot\rangle_{L^{2}_{\varrho}(D)}$ (recall (3.5)), this sheds further light on why minimizers of (3.4) converge to (3.6) (the orthogonal projection in the $L^{2}_{\varrho}$-inner product).

• Remark 3.4 (Numerical considerations)

  Fast numerical computations are not the primary concern of this article. However, we note that (3.7) can be solved using standard linear algebra techniques. Since the matrix $A$ is generally dense and unstructured, each matrix-vector multiplication involves $\mathcal{O}(mn)$ floating-point operations (flops). Hence, when using an iterative method such as conjugate gradients, the number of flops that suffice to compute $\hat{c}$ to an error of $\eta>0$ (in the norm ${\|A\cdot\|}_{2}$ )is roughly $\mathrm{cond}(A)\cdot m\cdot n\cdot\log(1/\eta)$, where $\mathrm{cond}(A)$ is the condition number of $A$. In §5 we see that the sufficient conditions that ensure accuracy and stability of the approximation $\hat{f}$ also guarantee that $A$ is well conditioned.

We now introduce some key terminology that will be used from now on. First, we say that the approximation $\hat{f}$ is $L^{2}_{\varrho}$-quasi-optimal or $(L^{2}_{\varrho},L^{\infty}_{\varrho})$-quasi-optimal if, in the absence of noise,

$${\|f-\hat{f}\|}_{L^{2}_{\varrho}(D)}\lesssim\inf_{p\in\mathcal{P}}{\|f-p\|}_{L^{2}_{\varrho}(D)},\quad\text{or}\quad{\|f-\hat{f}\|}_{L^{2}_{\varrho}(D)}\lesssim\inf_{p\in\mathcal{P}}{\|f-p\|}_{L^{\infty}_{\varrho}(D)},$$

respectively (note that the term instance optimality also sometimes used [47]). Obviously the former is stronger – therefore, achieving it will be our main objective. We say that the recovery is stable if, in the presence of noise, the recovery error depends on ${\|e\|}_{2}$, i.e.,

$${\|f-\hat{f}\|}_{L^{2}_{\varrho}(D)}\lesssim e_{\mathcal{P}}(f)+{\|e\|}_{2},$$

where $e_{\mathcal{P}}(f)$ is some best approximation error term. Finally, we say that a sampling strategy (i.e., a collection of measures $\mu_{1},\ldots,\mu_{m}$) has near-optimal sample complexity or optimal sample complexity if, respectively, $m\geq cn\log(n)$ or $m\geq cn$ samples suffice for obtaining a quasi-optimal and stable approximation, for some constant $c>0$. This is sometimes referred to as rate optimality [47].

Let $D\subseteq\mathbb{R}^{d}$ be a domain and $\varrho$ a measure. Let $S\subset\mathbb{N}^{d}_{0}$ be a finite set of multi-indices with $|S|=n$. We consider the polynomial space

$$\mathcal{P}=\mathbb{P}_{S}:=\mathrm{span}\left\{x\mapsto x^{\nu}:\nu\in S\right\}\subset L^{2}_{\varrho}(D).$$

(4.1)

Here, we use the notation $x=(x_{i})^{d}_{i=1}$ for the $d$-dimensional variable, $\nu=(\nu_{1},\ldots,\nu_{d})$ for a multi-index and $x^{\nu}=x^{\nu_{1}}_{1}\cdots x^{\nu_{d}}_{d}$. There are several standard choices for the index set $S$. In low dimensions, one may consider the (isotropic) tensor-product or total degree

$$S=S^{\mathsf{TP}}_{p}=\left\{\nu\in\mathbb{N}^{d}_{0}:\max_{k=1,\ldots,d}\nu_{k}\leq p\right\},\qquad S=S^{\mathsf{TD}}_{p}=\left\{\nu\in\mathbb{N}^{d}_{0}:\sum^{d}_{k=1}\nu_{k}\leq p\right\}$$

(4.2)

index sets of order $p\in\mathbb{N}_{0}$. Unfortunately, the cardinality of these index sets scales poorly with respect to $d$ (for fixed $p$). A better choice is moderate dimensions is the hyperbolic cross index set

$$S=S^{\mathsf{HC}}_{p}=\left\{\nu\in\mathbb{N}^{d}_{0}:\sum^{d}_{k=1}(\nu_{k}+1)\leq p+1\right\}.$$

(4.3)

However, as the dimension increases, this too may become too large to use. Since high-dimensional functions often have a very anisotropic dependence with respect to the coordinate variables, in higher dimensions one may look to consider anisotropic versions of these index sets. Given an anisotropy parameter $a=(a_{k})^{d}_{k=1}$ with $a>0$ (understood componentwise) and $p\geq 0$ (not necessarily an integer), the corresponding anisotropic index sets are defined as

$\displaystyle S^{\mathsf{TP}}_{p,a}=\left\{\nu\in\mathbb{N}^{d}_{0}:\max_{k=1,\ldots,d}a_{k}\nu_{k}\leq p\right\},\qquad S^{\mathsf{TD}}_{p,a}=\left\{\nu\in\mathbb{N}^{d}_{0}:\sum^{d}_{k=1}a_{k}\nu_{k}\leq p\right\}$

and

$$S^{\mathsf{HC}}_{p,a}=\left\{\nu\in\mathbb{N}^{d}_{0}:\sum^{d}_{k=1}(\nu_{k}+1)^{a_{k}}\leq p+1\right\}.$$

Notice that the isotropic index sets are recovered by setting $a=1$ (the vector of ones).

The choice of index set is not the focus of this paper. For more discussion, see, e.g., [64] and [5, Chpts. 2 & 5]. We remark, however, that all index sets defined above are examples of lower (also known as monotone or downward closed) sets. A set $S\subseteq\mathbb{N}^{d}_{0}$ is lower if whenever $\nu\in S$ and $\mu\leq\nu$ (understood componentwise, once more), one also has $\mu\in S$.

As we shall see in the next section, orthonormal bases play a key role in least-squares approximation from random samples and, in particular, optimal sampling. It is therefore convenient to be able to readily compute orthonormal polynomial bases for the space $\mathbb{P}_{S}$.

Fortunately, when $S$ is lower and $D$ and $\varrho$ are of tensor-product type, such polynomials are easily generated via tensor-products of univariate orthogonal polynomials. For concreteness, let

$$D=(a_{1},b_{1})\times\cdots\times(a_{d},b_{d}),\qquad\varrho=\rho_{1}\times\cdots\times\rho_{d},$$

where, for each $k=1,\ldots,d$, $-\infty\leq a_{k}<b_{k}\leq\infty$ and $\rho_{k}$ is a probability measure on $(a_{k},b_{d})$. Then, under mild conditions on $\rho_{k}$ (see, e.g., [112, §2.2]), there exists a unique sequence of orthonormal polynomials

$$\{\psi^{(k)}_{i}\}^{\infty}_{i=0}\subset L^{2}_{\rho_{k}}(a_{k},b_{k}),$$

where, for each $i$, $\psi^{(k)}_{i}$ is a polynomial of degree $i$. Using this, one immediately obtains an orthonormal basis of $L^{2}_{\varrho}(D)$ via tensor products. Specifically,

$$\{\Psi_{\nu}\}_{\nu\in\mathbb{N}^{d}_{0}}\subset L^{2}_{\varrho}(D),\quad\text{where }\Psi_{\nu}=\psi^{(1)}_{\nu_{1}}\otimes\cdots\otimes\psi^{(d)}_{\nu_{d}},\ \forall\nu=(\nu_{k})^{d}_{k=1}\in\mathbb{N}^{d}_{0}.$$

What about the subspace $\mathbb{P}_{S}$ introduced in (4.1)? Fortunately, whenever $S$ is a lower set, the functions $\Psi_{\nu}$ with indices $\nu\in S$ also form an orthonormal basis for this space. In other words,

$$\text{$S$ lower}\ \Longrightarrow\ \mathrm{span}\{\Psi_{\nu}:\nu\in S\}=\mathbb{P}_{S}.$$

See, e.g., [14, Proof of Thm. 6.5]. This property, combined with the tensor-product structure of the basis functions, makes optimal sampling and least-squares approximation in the subspaces $\mathbb{P}_{S}$ computationally feasible and, in many cases, straightforward, for tensor-product domains and measures. See §8 for some further discussion.

To conclude this section, we list several standard families of univariate measures and their corresponding orthogonal polynomials. Consider a compact interval, which without loss of generality we take to be $(-1,1)$. Given parameters $\alpha,\beta>-1$, the Jacobi (probability) measure is given by

$$\,\mathrm{d}\rho(x)=c_{\alpha,\beta}(1-x)^{\alpha}(1+x)^{\beta}\,\mathrm{d}y,\ y\in(-1,1),\qquad\text{where }c_{\alpha,\beta}=\left(\int^{1}_{-1}(1-x)^{\alpha}(1+x)^{\beta}\,\mathrm{d}x\right)^{-1}.$$

This measure generates the Jacobi polynomials for general $\alpha,\beta$ and the ultraspherical polynomials when $\alpha=\beta$. Of particular interest are the following cases.

• $\bullet$  

  $\alpha=\beta=-1/2$, which corresponds to the arcsine measure $\,\mathrm{d}\rho(x)=(\pi\sqrt{1-x^{2}})^{-1}\,\mathrm{d}y$. This yields the Chebyshev polynomials of the first kind.

• $\bullet$  

  $\alpha=\beta=0$, which corresponds to the uniform measure $\,\mathrm{d}\rho(x)=\frac{1}{2}\,\mathrm{d}x$. This yields the Legendre polynomials.

• $\bullet$  

  $\alpha=\beta=1/2$, which corresponds to the measure $\,\mathrm{d}\rho(x)=(2/\pi)\sqrt{1-x^{2}}\,\mathrm{d}x$. This yields the Chebyshev polynomials of the second kind.

We will consider these polynomials later in this paper. We will also briefly discuss certain unbounded domains. Here two common examples are as follows.

• $\bullet$  

  $\,\mathrm{d}\rho(x)=(2\pi)^{-1/2}\mathrm{e}^{-x^{2}/2}\,\mathrm{d}x$ over $\mathbb{R}$, which yields the Hermite polynomials.

• $\bullet$  

  $\,\mathrm{d}\rho(x)=\mathrm{e}^{-x}\,\mathrm{d}x$ over $[0,\infty)$, which yields the Laguerre polynomials.

Accuracy and stability of the weighted least-squares approximation are controlled by the following discrete stability constants:

$\displaystyle\alpha_{w}=\inf\left\{{\|p\|}_{\mathsf{disc},w}:p\in\mathcal{P},\ {\left\|p\right\|}_{L^{2}_{\varrho}(D)}=1\right\},\quad\beta_{w}=\sup\left\{{\|p\|}_{\mathsf{disc},w}:p\in\mathcal{P},\ {\left\|p\right\|}_{L^{2}_{\varrho}(D)}=1\right\}.$

In other words, $\alpha_{w}$ and $\beta_{w}$ are the optimal constants in the norm equivalence

$$\alpha_{w}{\left\|p\right\|}_{L^{2}_{\varrho}(D)}\leq{\|p\|}_{\mathsf{disc},w}\leq\beta_{w}{\left\|p\right\|}_{L^{2}_{\varrho}(D)},\quad\forall p\in\mathcal{P}.$$

(5.1)

In approximation theory, this is known as a sampling discretization [68] or a (weighted) Marcinkiewicz–Zygmund inequality [116]. Squaring and writing out the discrete semi-norm, (5.1) is equivalent to

$$\alpha^{2}_{w}{\left\|p\right\|}^{2}_{L^{2}_{\varrho}(D)}\leq\frac{1}{m}\sum^{m}_{i=1}w(x_{i})|p(x_{i})|^{2}\leq\beta^{2}_{w}{\left\|p\right\|}^{2}_{L^{2}_{\varrho}(D)},\quad\forall p\in\mathcal{P}.$$

(5.2)

Hence, the existence of finite, positive constants $0<\alpha_{w}\leq\beta_{w}<\infty$ implies that ${\left\|\cdot\right\|}_{\mathsf{disc},w}$ is a norm over the $n$-dimensional space $\mathcal{P}$, with $\alpha_{w},\beta_{w}$ being the constants of the norm equivalence.

Note that if $\{\phi_{i}\}^{n}_{i=1}$ is an orthonormal basis for $\mathcal{P}$, then it is straightforward to show that

$$\alpha_{w}=\sigma_{\min}(A)=\sqrt{\lambda_{\min}(A^{*}A)},\quad\beta_{w}=\sigma_{\max}(A)=\sqrt{\lambda_{\max}(A^{*}A)},$$

(5.3)

where $A$ is the least-squares matrix (3.8).

This result is a standard exercise. In particular, the condition number statement follows immediately from (5.3). We include a short proof of the other parts in the appendix for completeness. We also observe that this result holds for arbitrary weight functions $w$ and sample points $x_{1},\ldots,x_{m}$ satisfying the stipulated assumptions. At this stage, we do not require the sample points to be random. This will be used in the next subsection to derive concrete sample complexity estimates.

• Remark 5.2 (The noise bound)

  On the face of it, the noise term ${\left\|e\right\|}_{2,w}$ is undesirable since coefficients $e_{i}$ corresponding to large values of $w(x_{i})$ are more heavily weighted than others. We will take this into account later when we construct near-optimal sampling measures. Specifically, in §6.1 we construct sampling measures that lead to log-linear sample complexity and for which $w(x)\leq 2$. Hence, the noise term ${\left\|e\right\|}_{2,w}\leq\sqrt{2}{\left\|e\right\|}_{2}$ in this case.

We now return to the main setting of this paper, where the samples points are drawn randomly and independently with $x_{i}\sim\mu_{i}$, $i=1,\ldots,m$, for measures $\mu_{i}$ satisfying Assumption 3.1. Our aim is to analyze the sample complexity of weighted least-squares approximation. In view of Lemma 5.1, this involves first bounding the discrete stability constants $\alpha_{w}$ and $\beta_{w}$.

A key tool in this analysis is the Christoffel function of $\mathcal{P}$. Christoffel functions are well-known objects in approximation theory [98, 125], where they are typically considered in the context of spaces spanned by algebraic polynomials. It transpires that Christoffel functions – or, more precisely, their reciprocals – are also fundamentally associated with random sampling for least-squares approximation.

In other words, $\mathcal{K}(x)$ measures how large in magnitude an element of $\mathcal{P}$ can be at $x\in D$ in relation to its $L^{2}_{\varrho}$-norm. This function also admits an explicit expression. Let $\{\phi_{i}\}^{n}_{i=1}$ be an arbitrary orthonormal basis of $\mathcal{P}$. Then it is a short exercise to show that

$$\mathcal{K}(x)=\sum^{n}_{i=1}|\phi_{i}(x)|^{2}.$$

(5.6)

Often taken as the definition of $\mathcal{K}$, this formulation will be useful in our subsequent analysis. It also emphasizes the fact that $\mathcal{K}$ is precisely the diagonal of the reproducing kernel of $\mathcal{P}$ in $L^{2}_{\varrho}(D)$ [98, §3].

For reasons that will become clear soon, we are particularly interested in the maximal behaviour of the function $w(x)\mathcal{K}(x)$, where $w=1/\nu$ is the weight function defined by (3.1). We therefore let

$$\kappa_{w}=\kappa_{w}(\mathcal{P}):=\operatorname*{ess\,sup}_{x\sim\varrho}w(x)\mathcal{K}(\mathcal{P})(x).$$

(5.7)

To continue the connection with approximation theory, it is worth noting that $\kappa_{w}$ is the optimal constant in the (weighted) Nikolskii-type inequality (see, e.g., [92] and references therein),

$${\|\sqrt{w(\cdot)}p(\cdot)\|}_{L^{\infty}_{\varrho}(D)}\leq\sqrt{\kappa_{w}}{\|p\|}_{L^{2}_{\varrho}(D)},\quad\forall p\in\mathcal{P}.$$

(5.8)

Thus, $\kappa_{w}$ measures how large the scaled element $\sqrt{w(\cdot)}p(\cdot)$ can be uniformly in relation to the $L^{2}_{\varrho}$-norm of $p$. It is important to observe that

$$\kappa_{w}(\mathcal{P})\geq n,$$

(5.9)

for any weight function $w$ and $n$-dimensional subspace $\mathcal{P}$. This bound follows by integrating both sides with respect to the measure $\varrho$ and noticing that $\int_{D}\mathcal{K}(x)\,\mathrm{d}\varrho(x)=n$, the latter being an immediate consequence of (5.6). This gives

$$n=\int_{D}\mathcal{K}(x)\,\mathrm{d}\varrho(x)=\int_{D}w(x)\mathcal{K}(x)\frac{1}{w(x)}\,\mathrm{d}\varrho(x)\leq\kappa_{w}\int_{D}\frac{1}{w(x)}\,\mathrm{d}\varrho(x)=\kappa_{w},$$

where in the last equality we used (3.2) and the fact that $w=1/\nu$.

The following result establishes a key relationship between the Christoffel function and the sample complexity of weighted least-squares approximation.

This result is well known. In view of (5.3), its proof relies on bounding the maximum and minimum eigenvalues of $A^{*}A$. This is achieved by using what have now become quite standard matrix concentration inequalities, such as the matrix Chernoff bound [118, Thm. 1.1] (see also Theorem A.1). See the appendix for a proof.

• Remark 5.5 (One-sided estimates)

  The conclusions of Lemma 5.1 only rely on bounding the lower discrete stability constant $\alpha_{w}$ from below. This can be done with a slightly smaller sampling condition than (5.11). It follows readily from the proof of Theorem 5.4 that $\alpha_{w}>\sqrt{1-\delta}$ with probability at least $1-\epsilon$, whenever

  $$m\geq c^{\prime}_{\delta}\cdot\kappa_{w}(\mathcal{P})\cdot\log(n/\epsilon),\quad\text{where }c^{\prime}_{\delta}=((1-\delta)\log(1-\delta)+\delta)^{-1}.$$

  However, bounding $\beta_{w}$ from above yields a bound on the condition number of $A$ (see Lemma 5.1), which, as discussed in Remark 3.2, is important for numerical purposes.

We next combine Lemma 5.1 and Theorem 5.4 to obtain error bounds for weighted least-squares approximation. We split these bounds into two types: error bounds in probability (this subsection) and error bounds in expectation (the next subsection). In these two subsections, we will strive for generality by tracking the dependence in these bounds on the parameter $0<\delta<1$ appearing in Theorem 5.4. However, it is generally informative to think of this as a fixed scalar, e.g., $\delta=1/2$.

This result follows immediately from Lemma 5.1 and Theorem 5.4 via the estimate

$${\left\|f-p\right\|}_{\mathsf{disc},w}\leq c_{w}{\|f-p\|}_{L^{\infty}_{\varrho}(D)}.$$

Now suppose that $\delta=1/2$ (for concreteness) and assume further that $w(x)\lesssim 1$, a.e. $x\sim\varrho$. This will be the case in Section 6 when we construct near-optimal sampling measures. Then (5.12) reads

$${\|f-\hat{f}\|}_{L^{2}_{\varrho}(D)}\lesssim\inf_{p\in\mathcal{P}}{\|f-p\|}_{L^{\infty}_{\varrho}(D)}+{\left\|e\right\|}_{2}.$$

Hence the approximation is stable and $(L^{2}_{\varrho},L^{\infty}_{\varrho})$-quasi-optimal. In some problems, the difference between the $L^{2}_{\varrho}$- and $L^{\infty}_{\varrho}$-norms may be of little consequence. For example, in the case of polynomial approximation of holomorphic functions in low dimensions, the best approximation error decays exponentially fast with respect to $n$ in either norm (see, e.g., [5, §3.5-3.6]). On the other hand, for high-dimensional holomorphic functions or functions of finite regularity in any dimension, the best approximation errors decay algebraically fast, with, typically, the $L^{\infty}_{\varrho}$-norm error decaying at least $\mathcal{O}(\sqrt{n})$ slower than the $L^{2}_{\varrho}$-norm error (see, e.g., [5, §3.8-3.9]). Thus, the crude bound (5.12) may underestimate the convergence rate of the least-squares approximation. Motivated by these considerations, we next discuss how to establish $L^{2}_{\varrho}$-quasi-optimality results.

• Remark 5.7 (Uniform versus nonuniform)

  Corollary 5.6 is a uniform result, in the sense that a single random draw of the sample points suffices for all functions. We next discuss a nonuniform results, in which the error bound holds with high probability for each fixed function. Uniform bounds are desirable in many applications, as it means that the same sample points (which may correspond to, e.g., sensor locations) can be re-used for approximating multiple functions. Theoretically, it also means that one can achieve worst-case error bounds. Indeed, let $\mathcal{F}\subset L^{2}_{\varrho}(D)$ be a set of functions that are defined everywhere and for which

  $$E_{\mathcal{P}}(\mathcal{F})=\sup_{f\in\mathcal{F}}\inf_{p\in\mathcal{P}}{\|f-p\|}_{L^{\infty}_{\varrho}(D)}<\infty.$$

  Typically, $\mathcal{F}$ may be a unit ball of some Banach space – for example, the space of Sobolev functions $H^{k}_{\varrho}(D)$. Then, ignoring noise for simplicity and assuming as before that $\delta=1/2$ and $w(x)\lesssim 1$, a.e. $x\sim\varrho$, Corollary 5.6 implies the following uniform bound with high probability:

  $$\sup_{f\in\mathcal{F}}{\|f-\hat{f}\|}_{L^{2}_{\varrho}(D)}\lesssim E_{\mathcal{P}}(\mathcal{F}).$$

  As we discuss in §8, this has implications in the study of sampling numbers in information-based complexity and the efficacy of pointwise samples (standard information).

• Remark 5.8 (The term $c_{w}$)

  As an alternative to assuming that $w(x)\lesssim 1$, a.e. $x\sim\varrho$, one may also bound the term $c_{w}$ by assuming that $\mathcal{P}$ contains a function $h$ with ${\|h\|}_{L^{2}_{\varrho}(D)}=1$ and $h(x)\gtrsim 1$, a.e. $x\sim\varrho$. This holds, for example, whenever the constant function $1\in\mathcal{P}$. In this case,

  $$c_{w}\lesssim{\left\|h\right\|}_{\mathsf{disc},w}\leq\beta_{w}{\left\|h\right\|}_{L^{2}_{\varrho}(D)}\leq\sqrt{1+\delta}.$$

  However, as noted in Remark 5.1, we can always construct $w$ so that the former assumption holds.

We now present a nonuniform ‘in probability’ bound that provides $L^{2}_{\varrho}$-quasi-optimality, at the expense of a poor scaling with respect to the failure probability $\epsilon$. The proof is based on Markov inequality, which, roughly speaking, is used to bound the discrete error term arising in (5.4).

Suppose again that $\delta=1/2$ and $w(x)\lesssim 1$, $\forall x$. Then this bound implies that

$${\|f-\hat{f}\|}_{L^{2}_{\varrho}(D)}\lesssim(1+1/\sqrt{\epsilon\log(4n/\epsilon)})\inf_{p\in\mathcal{P}}{\|f-p\|}_{L^{2}_{\varrho}(D)}+{\left\|e\right\|}_{2}.$$

While stable and $L^{2}_{\varrho}$-quasi-optimal, the scaling with respect to $\epsilon$ is undesirable. To obtain an $\epsilon$-independent bound, this suggests we either need to $n$ to be exponentially large in $1/\epsilon$, or impose an additional constraint on $m$ that $m\gtrsim\kappa_{w}(\mathcal{P})/\epsilon$. Neither is a desirable outcome.

One possible way to circumvent this issue involves using Bernstein’s inequality instead of Markov’s inequality. This exploits the fact that the discrete term in (5.4) is a sum of independent random variables with bounded variance. It therefore concentrates exponentially fast (in $m$) around its mean $\mathbb{E}{\|f-p\|}^{2}_{\mathsf{disc},w}={\|f-p\|}^{2}_{L^{2}_{\varrho}(D)}$ (recall (3.5)). This leads to the following result, which is also nonuniform.

This result asserts a mixed type of quasi-optimality, involving the $L^{2}_{\varrho}$-norm and a (weighted) $L^{\infty}$-norm divided by the factor $\sqrt{k}$. Notice that the factor $\sqrt{w}$ can be removed whenever $w(x)\lesssim 1$, a.e. $x\sim\varrho$, as will be the case later. Therefore, consider, as in Remark 5.4, a case where the $L^{2}_{\varrho}$-norm best approximation error decays algebraically fast in $n=\dim(\mathcal{P})$. As we noted therein, the $L^{\infty}_{\varrho}$-norm best approximation error often decays $\mathcal{O}(\sqrt{n})$ slower than the former. Hence, one may choose $k=n$ in (5.15) to show that $\hat{f}$ achieves the same algebraic convergence rate in $L^{2}_{\varrho}$-norm as the $L^{2}_{\varrho}$-norm best approximation in $\mathcal{P}$. This approach has also been used in the related context of function approximation via compressed sensing in [110] and [5, §7.6].

To obtain error bounds in expectation, we need to modify the least-squares estimator to avoid the ‘bad’ regime where the discrete stability constants can be poorly behaved. In this section, we proceed as in [47, §2.2], which is based on [34, 35].

Let $\{\phi_{i}\}^{n}_{i=1}$ be an orthonormal basis of $\mathcal{P}$. First, we notice that (5.10) holds whenever

$${\|G-I\|}_{2}\leq\delta,$$

where $G$ is the discrete Gram matrix

$$G=\left(\langle\phi_{j},\phi_{k}\rangle_{\mathsf{disc},w}\right)^{n}_{j,k=1}\in\mathbb{C}^{n\times n}.$$

It is then possible to show the following bound.