Optimal sampling for least-squares approximation
Ben Adcock




Fair Minimum Representation Clustering via Integer Programming A preliminary version of this paper appeared in CPAIOR 2024.
    [
    
===========================================================================================================================



§ ABSTRACT 
Least-squares approximation is one of the most important methods for recovering an unknown function from data. While in many applications the data is fixed, in many others there is substantial freedom to choose where to sample. In this paper, we review recent progress on optimal sampling for (weighted) least-squares approximation in arbitrary linear spaces. We introduce the Christoffel function as a key quantity in the analysis of (weighted) least-squares approximation from random samples, then show how it can be used to construct sampling strategies that possess near-optimal sample complexity: namely, the number of samples scales log-linearly in n, the dimension of the approximation space. We discuss a series of variations, extensions and further topics, and throughout highlight connections to approximation theory, machine learning, information-based complexity and numerical linear algebra. Finally, motivated by various contemporary applications, we consider a generalization of the classical setting where the samples need not be pointwise samples of a scalar-valued function, and the approximation space need not be linear. We show that even in this significantly more general setting suitable generalizations of the Christoffel function still determine the sample complexity. This provides a unified procedure for designing improved sampling strategies for general recovery problems. This article is largely self-contained, and intended to be accessible to nonspecialists.






§ INTRODUCTION


Least-squares approximation is the process of fitting an unknown function from samples by computing a best ℓ^2-norm fit in a given subspace, which is often termed the approximation space. Least squares is a classical topic, yet it is one of the widely-used tools in applied mathematics, computer science, engineering and numerous other disciplines. For the data scientist, it is almost always ones first `go-to' method when trying to fit a function to data.

In many data-fitting problems, the samples are fixed. However, many other problems offer substantial flexibility to choose where to sample. When data is also expensive to acquire – which, despite claims about `big data' is often the case in applications in science and engineering – we are naturally led to the following questions. How many samples do we need – or, in other words, what is the sample complexity – and how should we best choose them? This is by no means a new question. It arises in many different guises in different fields, including optimal design of experiments in statistics, active learning in machine learning, optimal sensor placement in sampling theory and signal processing, and optimal (standard) information in information-based complexity.

The purpose of this article is to survey recent advances made in the last 5-10 years in optimal sampling, as we shall term it from now on, which has been motivated by certain function approximation problems in high dimensions. Such methods are in essence importance sampling techniques, where samples are drawn randomly from a probability measure chosen specifically for the given approximation space. Throughout, our aim is to ensure quasi-optimal recovery (in an appropriate sense) with near-optimal sample complexity. 



 §.§ Overview




After a short literature review (<ref>), this article commences with a formulation and review of (weighted) least-squares approximation (<ref>). We then discuss multivariate polynomial approximation (<ref>), 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 (<ref>). We then show that sampling from a measure proportional to this function leads to provably near-optimal sampling (<ref>). 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 (<ref>). 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 (<ref>). 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 (<ref>). 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 (<ref>).




 §.§ Scope and target audience


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. 











§ LITERATURE REVIEW


We commence with a short discussion of relevant literature. Additional literature on variations, extensions and further topics can be found in <ref>.

Least squares is a classical topic, with origins tracing back to the work of Gauss and Legendre <cit.>. Starting in the early 2010s, and motivated by problems in parametric and stochastic Differential Equations (DEs), there was a resurgence of research on this topic, focusing on high- and and infinite-dimensional function approximation, and typically involving polynomial spaces. Key works in this direction include <cit.>. This resurgence was based on least squares with random sampling, inspired by Monte Carlo quadrature and its ability to integrate functions without succumbing to the curse of dimensionality. However, it is worth noting that the goal of least-squares approximation is to achieve quasi-optimal rates of convergence with respect to the approximation space. Typically, the resulting rate will exceed the error rate for Monte Carlo quadrature.

As noted above, Monte Carlo sampling generically leads to suboptimal sample complexity bounds for least-squares approximation. This observation led to a concerted effort to develop practical sampling strategies with better performance (see <cit.> for an overview), culminating in the near-optimal strategies which are the basis of this work. These were developed in <cit.>, but also appeared slightly earlier in <cit.> in the case of (total degree) polynomial spaces. 

At a similar time, related techniques under the name leverage score sampling – which are based on the classical topic of statistical leverage – have become increasingly popular in machine learning and data science. In particular, leverage score sampling is an effective tool for matrix sketching <cit.>. As we comment in <ref>, it is can also be viewed as a special case of the techniques described in this article, corresponding to functions defined over a discrete domain. 

Finally, we remark on some applications. As observed, this work is closely related to optimal design of experiments and optimal sensor placement in sampling theory and signal processing – both large areas with countless applications that we shall not attempt to review. However, this specific line of research emerged out of computing polynomial approximations to high-dimensional functions arising in parametric and stochastic DEs <cit.>, and this remains a key area of application. See <cit.> and references therein. 
For other surveys focused multivariate polynomial approximation and parametric and stochastic DEs, see <cit.> and <cit.>.


Recently, these techniques have also been applied to the closely related problem of numerical integration (cubature) <cit.>.  There are also emerging applications in Trefftz methods for solving Helmoltz equations <cit.> and methods for option pricing in finance <cit.>. On the theoretical side, this line of work has also spurred recent advances in approximation theory (so-called sampling discretizations) and information-based complexity (so-called sampling numbers). We discuss these topics further in <ref>. Related ideas have also been used in sampling theory <cit.>. We also note that Christoffel functions are useful tools for empirical inference in data analysis <cit.>.

Finally, through the close connection to leverage score sampling, there are manifold applications in machine learning and data science. These include randomized numerical linear algebra <cit.>, kernel methods <cit.> and active learning <cit.>.
Moreover, the generalization we describe in <ref> opens the door to applications in many seemingly unrelated areas, such as inverse problems in imaging <cit.>.













§ PRELIMINARIES


Let (D,,ϱ) be a measure space and L^2_ϱ(D) be the Lebesgue space of square-integrable functions f : D → with respect to ϱ. Typically, in this work, D ⊆^d.  For convenience, we assume that ϱ is a finite measure (ϱ(D) < ∞) and, therefore, without loss of generality, that ϱ is a probability measure (ϱ(D) = 1). It is possible to consider infinite measures, but for ease of exposition we shall not do this.

Given m ∈, we consider sampling measures μ_1,…,μ_m. These are assumed to be such that (D,,μ_i) is a probability space for every i. We also make the following assumption.

[Absolute continuity and positivity]

The additive mixture

μ= 1/m ∑^m_i=1 μ_i

is absolutely continuous with respect to ϱ and its Radon–Nikodym derivative ν is strictly positive almost everywhere on supp(ϱ).


This assumption allows us to write


1/m ∑^m_i=1 μ_i(x) = ν(x) ϱ(x),
 
where the density ν : D → (the Radon–Nikodym derivative) is measurable, positive almost everywhere and satisfies


∫_D ν(x) ϱ(x) = 1.

In what follows it will often be more convenient to work with the reciprocal of this function. We define the weight function w : D → as w(x) = 1/ν(x), x ∈ D.

Given sampling measures μ_1,…,μ_m, we now draw samples x_i ∼μ_i, i = 1,…,m, independently from these measures and consider noisy measurements of an unknown function f : D → the form  


y_i = f(x_i) + e_i ,  i = 1,…,m.

Typically, we will assume that f ∈ L^2_ϱ(D) so that the samples f_meas are almost surely well defined.

We consider a bounded, adversarial noise model, where the e_i's are not assumed to be random, but are assumed to be small in magnitude. Our aim is to derive error bounds in which the noise term scales like the ℓ^2-norm

e_2 = √(∑^m_i=1 |e_i|^2)

of the noise vector e = (e_i)^m_i=1. Random noise models (included unbounded noise)  can also be considered (see <cit.> and <cit.>). 










 §.§ Weighted least-squares approximation


Let ⊂ L^2_ϱ(D) be an arbitrary n-dimensional subspace, where n ≤ m, in which we seek to approximate the unknown f from the measurements f_meas. We term  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 <cit.>. Instead, we consider the weighted least-squares approximation


f̂ ∈p ∈ 1/m ∑^m_i=1 w(x_i) | y_i - p(x_i) |^2.

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


[The scaling factor]
The scaling factors in wls-prob are motivated by noticing that


1/m ∑^m_i=1 w(x_i) | g(x_i) |^2 = 1/m ∑^m_i=1 ∫_D w(x) |g(x)|^2 μ_i(x)
 =  ∫_D |g(x)|^2 ρ(x) = g^2_L^2_ϱ(D),

where the second equality is due to mu_weight_fn. Thus, in the noiseless case, wls-prob can be considered as a empirical approximation to the continuous least-squares approximation


f̂ = p ∈ f - p^2_L^2_ϱ(D),

i.e., the best approximation to f from  in the L^2_ϱ(D)-norm. In particular, if μ_1 = ⋯ = μ_m = μ, then the minimizers of wls-prob converge almost surely to the minimizer of cts-min as m →∞ <cit.>.


The objective of this article is to describe how to choose the measures μ_1,…,μ_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 ϱ. This is equivalent to setting

μ_1 = ⋯= μ_m = ϱ,

which leads, via mu_weight_fn, to ν≡ 1. Thus, wls-prob becomes an unweighted least-squares approximation. 




[Hierarchical approximation]

Often, rather than a fixed subspace , one may wish to construct a sequence of approximations in a nested collection of subspaces

^(1) ⊆^(2) ⊆⋯,

of dimension (^(k)) = n_k. In this case, given numbers 1 ≤ m_1 ≤ m_2 ≤⋯ satisfying m_k ≥ n_k, ∀ k, one aims to design nested collections of sample points

{ x^(1)_i }^m_1_i=1 ⊆{ x^(2)_i }^m_2_i=1 ⊆⋯.

We write f̂^(1),f̂^(2),… for the ensuing (weighted) least-squares approximations, where 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 <cit.> or sequential <cit.>) approximation scheme. 




 §.§ Reformulations of wls-prob


Given a basis {ϕ_i }^n_i=1 for , the problem wls-prob is readily reformulated as an algebraic least-squares problem for the coefficients ĉ = (ĉ_i)^n_i=1∈^n of f̂ = ∑^n_i=1ĉ_i ϕ_i. This takes the form


ĉ ∈c ∈^n A c - b^2_2,

where


A = ( √(w(x_i)/m) ϕ_j(x_i) )^m,n_i,j=1 ∈^m ×n,    b = ( √(w(x_i)/m) (f(x_i) + e_i ) )^m_i=1 ∈^m.

To be precise, every minimizer f̂ satisfying wls-prob has coefficients ĉ that satisfy algls-prob and vice versa. Classical least-squares analysis asserts that any vector ĉ satisfying algls-prob is also a solution of the normal equations


A^* A c = A^* b,

and vice versa. Rewriting the normal equations in terms of functions also leads to the following variational form of wls-prob. 


Find f̂∈ such that f̂p_𝖽𝗂𝗌𝖼,w = fp_𝖽𝗂𝗌𝖼,w + 1/m ∑^m_i=1 w(x_i) e_i p(x_i), ∀p ∈.

This is equivalent to normal-eqns in the same sense as before. Here we wrote


gh_𝖽𝗂𝗌𝖼,w = 1/m ∑^m_i=1 w(x_i) g(x_i) h(x_i),

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_𝖽𝗂𝗌𝖼,w = 1/m ∑^m_i=1 w(x_i) | g(x_i) |^2.

In the noiseless case e =0, the formulation variational-form asserts that f̂ is precisely the orthogonal projection of f onto  with respect to the discrete semi-inner product semi-inner-product. Since semi-inner-product is an empirical approximation to the continuous inner product  ··_L^2_ϱ(D) (recall exp-sum-scaling), this sheds further light on why minimizers of wls-prob converge to cts-min (the orthogonal projection in the L^2_ϱ-inner product).


[Numerical considerations]

Fast numerical computations are not the primary concern of this article. However, we note that algls-prob can be solved using standard linear algebra techniques. Since the matrix A is generally dense and unstructured, each matrix-vector multiplication involves m n floating-point operations (flops). Hence, when using an iterative method such as conjugate gradients, the number of flops that suffice to compute ĉ to an error of η > 0 (in the norm A ·_2 )is roughly cond(A) · m · n ·log(1/η), where cond(A) is the condition number of A. In <ref> we see that the sufficient conditions that ensure accuracy and stability of the approximation f̂ also guarantee that A is well conditioned.





 §.§ Key terminology


We now introduce some key terminology that will be used from now on. First, we say that the approximation f̂ is L^2_ϱ-quasi-optimal or (L^2_ϱ,L^∞_ϱ)-quasi-optimal if, in the absence of noise,

f - f̂_L^2_ϱ(D) ≲inf_p ∈ f - p_L^2_ϱ(D),  or  f - f̂_L^2_ϱ(D) ≲inf_p ∈ f - p_L^∞_ϱ(D),

respectively (note that the term instance optimality also sometimes used <cit.>). 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 - f̂_L^2_ϱ(D) ≲e_(f) + e_2,

where e_(f) is some best approximation error term. Finally, we say that a sampling strategy (i.e., a collection of measures μ_1,…,μ_m) has near-optimal sample complexity or optimal sample complexity if, respectively, m ≥ c n log(n) or m ≥ c n samples suffice for obtaining a quasi-optimal and stable approximation, for some constant c > 0. This is sometimes referred to as rate optimality <cit.>.










§ APPLICATION TO MULTIVARIATE POLYNOMIAL APPROXIMATION


We now introduce an important example considered in this paper, namely, multivariate polynomial approximation in d ≥ 1 dimensions.



 §.§ Spaces of multivariate polynomials


Let D ⊆^d be a domain and ϱ a measure. Let S ⊂^d_0 be a finite set of multi-indices with |S| = n. We consider the polynomial space


= _S  : = { x ↦x^ν : ν∈S } ⊂L^2_ϱ(D).

Here, we use the notation x = (x_i)^d_i=1 for the d-dimensional variable, ν = (ν_1,…,ν_d) for a multi-index and x^ν = x^ν_1_1⋯ x^ν_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^𝖳𝖯_p = { ν∈^d_0 : max_k=1,…,d ν_k ≤p  },    S = S^𝖳𝖣_p = { ν∈^d_0 : ∑^d_k=1 ν_k ≤p }

index sets of order p ∈_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^𝖧𝖢_p  =  { ν∈^d_0 : ∑^d_k=1 (ν_k+1) ≤p +1 }.

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 ≥ 0 (not necessarily an integer), the corresponding anisotropic index sets are defined as

S^𝖳𝖯_p,a = { ν∈^d_0 : max_k=1,…,d a_k ν_k ≤p  },
    S^𝖳𝖣_p,a = { ν∈^d_0 : ∑^d_k=1 a_k ν_k ≤p }

and

S^𝖧𝖢_p,a  =  { ν∈^d_0 : ∑^d_k=1 (ν_k+1)^a_k ≤p +1 }.

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., <cit.> and <cit.>. We remark, however, that all index sets defined above are examples of lower (also known as monotone or downward closed) sets. A set S ⊆^d_0 is lower if whenever ν∈ S and μ≤ν (understood componentwise, once more), one also has μ∈ S.







 §.§ Multivariate orthogonal polynomials on tensor-product domains


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 _S.

Fortunately, when S is lower and D and ϱ 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) ×⋯×(a_d , b_d),    ϱ= ρ_1 ×⋯×ρ_d,

where, for each k = 1,…,d,  -∞≤ a_k < b_k ≤∞ and ρ_k is a probability measure on (a_k,b_d). Then, under mild conditions on ρ_k (see, e.g., <cit.>), there exists a unique sequence of orthonormal polynomials

{ ψ^(k)_i }^∞_i=0 ⊂L^2_ρ_k(a_k,b_k),

where, for each i, ψ^(k)_i is a polynomial of degree i. 
Using this, one immediately obtains an orthonormal basis of L^2_ϱ(D) via tensor products. Specifically,

{ Ψ_ν }_ν∈^d_0 ⊂L^2_ϱ(D),  where Ψ_ν = ψ^(1)_ν_1 ⊗⋯⊗ψ^(d)_ν_d, ∀ν= (ν_k)^d_k=1 ∈^d_0.

What about the subspace _S introduced in PS-def? Fortunately, whenever S is a lower set, the functions Ψ_ν with indices ν∈ S also form an orthonormal basis for this space. In other words,

S lower ⟹ { Ψ_ν : ν∈S } = _S.

See, e.g., <cit.>.
This property, combined with the tensor-product structure of the basis functions, makes optimal sampling and least-squares approximation in the subspaces _S computationally feasible and, in many cases, straightforward, for tensor-product domains and measures. See <ref> 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 α,β > -1, the Jacobi (probability) measure is given by

ρ(x) = c_α,β (1-x)^α (1+x)^β y, y ∈(-1,1),    where c_α,β = ( ∫^1_-1 (1-x)^α (1+x)^β x )^-1.

This measure generates the Jacobi polynomials for general α,β and the ultraspherical polynomials when α = β. Of particular interest are the following cases.


  * α = β = -1/2, which corresponds to the arcsine measure ρ(x) = ( π√(1-x^2))^-1 y. This yields the Chebyshev polynomials of the first kind.

  * α = β = 0,  which corresponds to the uniform measure ρ(x) = 1/2 x. This yields the Legendre polynomials. 

  * α = β = 1/2,  which corresponds to the measure ρ(x) = (2 / π) √(1-x^2) 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.


  * ρ(x) = (2 π)^-1/2^-x^2/2 x over , which yields the Hermite polynomials.

  * ρ(x) = ^-x x over [0,∞), which yields the Laguerre polynomials.








§ THEORY OF WEIGHTED LEAST-SQUARES APPROXIMATION FROM RANDOM SAMPLES




 §.§ Basic accuracy and stability guarantee


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

α_w = inf{ p_𝖽𝗂𝗌𝖼,w : p ∈, p_L^2_ϱ(D) = 1 },
  β_w =  sup{ p_𝖽𝗂𝗌𝖼,w : p ∈, p_L^2_ϱ(D) = 1 }.

In other words, α_w and β_w are the optimal constants in the norm equivalence


α_w p_L^2_ϱ(D) ≤p_𝖽𝗂𝗌𝖼,w ≤β_w p_L^2_ϱ(D),  ∀p ∈.

In approximation theory, this is known as a sampling discretization <cit.> or a (weighted) Marcinkiewicz–Zygmund inequality <cit.>. Squaring and writing out the discrete semi-norm, MZ-inequality is equivalent to


α^2_w p^2_L^2_ϱ(D) ≤1/m ∑^m_i=1 w(x_i) | p(x_i) |^2 ≤β^2_w p^2_L^2_ϱ(D),  ∀p ∈.

Hence, the existence of finite, positive constants 0 < α_w ≤β_w < ∞ implies that ·_𝖽𝗂𝗌𝖼,w is a norm over the n-dimensional space , with α_w,β_w being the constants of the norm equivalence. 

Note that if {ϕ_i }^n_i=1 is an orthonormal basis for , then it is straightforward to show that


α_w = σ_min(A) = √(λ_min(A^*A)),  β_w = σ_max(A) = √(λ_max(A^*A)),

where A is the least-squares matrix (<ref>). 







[Accuracy and stability of weighted least squares]

Let ⊂ L^2_ϱ(D), f ∈  L^2_ϱ(D), x_1,…,x_m ∈ D be sample points at which both f and any p ∈ are finite, e ∈^m and w : D → be such that w(x_i) > 0, ∀ i ∈{1,…,m}. Suppose that α_w > 0. Then the weighted least-squares problem wls-prob has a unique solution f̂. Moreover, this solution satisfies


f - f̂_L^2_ϱ(D) ≤inf_p ∈ { f - p_L^2_ϱ(D) + 1/α_w f - p_𝖽𝗂𝗌𝖼,w } + 1/α_w e_2,w,

where e_2,w = √(1/m∑^m_i=1 w(x_i) | e_i |^2 ).
Also, if {ϕ_i }^n_i=1 is an orthonormal basis of  then the condition number of the least-squares matrix ls-Ab satisfies cond(A) = β_w / α_w.


This result is a standard exercise. In particular, the condition number statement follows immediately from alpha-beta-sigma. 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,…,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.


[The noise bound]

On the face of it, the noise term e_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 <ref> we construct sampling measures that lead to log-linear sample complexity and for which w(x) ≤ 2. Hence, the noise term e_2,w≤√(2)e_2 in this case.















 §.§ The (reciprocal) Christoffel function


We now return to the main setting of this paper, where the samples points are drawn randomly and independently with x_i ∼μ_i, i=1,…,m, for measures μ_i satisfying Assumption <ref>. Our aim is to analyze the sample complexity of weighted least-squares approximation. In view of Lemma <ref>, this involves first bounding the discrete stability constants α_w and β_w.

A key tool in this analysis is the Christoffel function of . Christoffel functions are well-known objects in approximation theory <cit.>, 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.

[Christoffel function]

Let ⊆ L^2_ϱ(D). The (reciprocal) Christoffel function of  is the function = () : D → defined by


(x) = ()(x) : = sup{ | p(x) |^2/p^2_L^2_ϱ(D) : p ∈, p ≠0 },  ∀x ∈D.



In other words, (x) measures how large in magnitude an element of  can be at x ∈ D in relation to its L^2_ϱ-norm.
This function also admits an explicit expression. Let {ϕ_i }^n_i=1 be an arbitrary orthonormal basis of . Then it is a short exercise to show that


(x) = ∑^n_i=1 | ϕ_i(x) |^2.

Often taken as the definition of , this formulation will be useful in our subsequent analysis. It also emphasizes the fact that  is precisely the diagonal of the reproducing kernel of  in L^2_ϱ(D) <cit.>. 

For reasons that will become clear soon, we are particularly interested in the maximal behaviour of the function w(x) (x), where w = 1/ν is the weight function defined by mu_weight_fn. We therefore let


κ_w = κ_w() : = _x ∼ϱ w(x) ()(x).

To continue the connection with approximation theory, it is worth noting that κ_w is the optimal constant in the (weighted) Nikolskii-type inequality (see, e.g., <cit.> and references therein),


√(w(·)) p(·) _L^∞_ϱ(D) ≤√(κ_w) p_L^2_ϱ(D),  ∀p ∈.

Thus, κ_w measures how large the scaled element √(w(·)) p(·) can be uniformly in relation to the L^2_ϱ-norm of p.
It is important to observe that


κ_w() ≥n,

for any weight function w and n-dimensional subspace . This bound follows by integrating both sides with respect to the measure ϱ and noticing that ∫_D(x) ϱ(x) = n, the latter being an immediate consequence of Kappa-def-alt. This gives

n = ∫_D (x) ϱ(x) =  ∫_D w(x) (x) 1/w(x) ϱ(x) ≤κ_w ∫_D 1/w(x) ϱ(x) = κ_w,

where in the last equality we used w_normalization and the fact that w = 1/ν.







 §.§ Bounding the discrete stability constants


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



[Estimates for α_w and β_w in probability]

Let 0 < δ,ϵ < 1, ⊂ L^2_ϱ(D) be a finite-dimensional subspace with () = n and μ_1,…,μ_m be probability measures satisfying Assumption <ref>. Consider sample points drawn randomly and independently with x_i ∼μ_i, i = 1,…,m. Then


√(1-δ) < α_w ≤β_w < √(1+δ),

with probability at least 1-ϵ, provided


m ≥c_δ ·κ_w() ·log(2 n/ϵ),  where c_δ = ((1+δ) log(1+δ) - δ))^-1

and w = 1/ν and κ_w are as in mu_weight_fn and kappa-w-def, respectively.


This result is well known. In view of alpha-beta-sigma, 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 <cit.> (see also Theorem <ref>).  See the appendix for a proof.


[One-sided estimates]
The conclusions of Lemma <ref> only rely on bounding the lower discrete stability constant α_w from below. This can be done with a slightly smaller sampling condition than m-bound-alpha-beta. It follows readily from the proof of Theorem <ref> that α_w > √(1-δ)
with probability at least 1-ϵ, whenever

m ≥c'_δ ·κ_w() ·log(n/ϵ),  where c'_δ = ((1-δ) log(1-δ) + δ)^-1.

However, bounding β_w from above yields a bound on the condition number of A (see Lemma <ref>), which, as discussed in Remark <ref>, is important for numerical purposes.







 §.§ Error bounds in probability




We next combine Lemma <ref> and Theorem <ref> 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 < δ < 1 appearing in Theorem <ref>. However, it is generally informative to think of this as a fixed scalar, e.g., δ = 1/2.



[First uniform error bound error bound in probability]

Let 0 < δ,ϵ < 1, ⊂ L^2_ϱ(D) be a finite-dimensional subspace with () = n and μ_1,…,μ_m be probability measures satisfying Assumption <ref>. Consider sample points drawn randomly and independently with x_i ∼μ_i, i = 1,…,m, where 

m ≥c_δ ·κ_w() ·log(4 n /ϵ),    c_δ = ((1+δ) log(1+δ) - δ))^-1

and w= 1/ν and κ_w are as in mu_weight_fn and kappa-w-def, respectively. Then the following hold with probability at least 1-ϵ. For any e ∈^m and f ∈ L^2_ϱ(D) that is defined everywhere in D, the weighted least-squares approximation f̂ is unique and satisfies


f - f̂_L^2_ϱ(D) ≤( 1 +  1/√(1-δ) c_w )  inf_p ∈ f - p_L^∞_ϱ(D) + 1/√(1-δ) e_2,w,

where c_w = √(1/m∑^m_i=1 w(x_i)). Moreover, if {ϕ_i }^n_i=1 is an orthonormal basis of , then the condition number of the least-squares matrix ls-Ab satisfies cond(A) ≤√(1+δ/1-δ).



This result follows immediately from Lemma <ref> and Theorem <ref> via the estimate

f-p_𝖽𝗂𝗌𝖼,w ≤c_w f - p_L^∞_ϱ(D).

Now suppose that δ = 1/2 (for concreteness) and assume further that w(x) ≲ 1, a.e. x ∼ϱ. This will be the case in Section <ref> when we construct near-optimal sampling measures. Then ls-err-bd-prob-unif-1 reads

f - f̂_L^2_ϱ(D) ≲inf_p ∈ f - p_L^∞_ϱ(D) + e_2.

Hence the approximation is stable and (L^2_ϱ,L^∞_ϱ)-quasi-optimal. 
In some problems, the difference between the L^2_ϱ- and L^∞_ϱ-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., <cit.>). 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^∞_ϱ-norm error decaying at least √(n) slower than the L^2_ϱ-norm error (see, e.g., <cit.>). Thus, the crude bound ls-err-bd-prob-unif-1 may underestimate the convergence rate of the least-squares approximation. Motivated by these considerations, we next discuss how to establish L^2_ϱ-quasi-optimality results.


[Uniform versus nonuniform]

Corollary <ref> 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 ⊂ L^2_ϱ(D) be a set of functions that are defined everywhere and for which

E_() = sup_f ∈ inf_p ∈ f - p_L^∞_ϱ(D) < ∞.

Typically,  may be a unit ball of some Banach space – for example, the space of Sobolev functions H^k_ϱ(D). Then, ignoring noise for simplicity and assuming as before that δ = 1/2 and w(x) ≲ 1, a.e. x ∼ϱ, Corollary <ref> implies the following uniform bound with high probability:

sup_f ∈ f - f̂_L^2_ϱ(D) ≲E_().
 
As we discuss in <ref>, this has implications in the study of sampling numbers in information-based complexity and the efficacy of pointwise samples (standard information).


[The term c_w]

As an alternative to assuming that w(x) ≲ 1, a.e. x ∼ϱ, one may also bound the term c_w by assuming that  contains a function h with h_L^2_ϱ(D) = 1 and h(x) ≳ 1, a.e. x ∼ϱ. This holds, for example, whenever the constant function 1 ∈. In this case,

c_w ≲h_𝖽𝗂𝗌𝖼,w ≤β_w h_L^2_ϱ(D) ≤√(1+δ).

However, as noted in Remark <ref>, we can always construct w so that the former assumption holds.


We now present a nonuniform `in probability' bound that provides L^2_ϱ-quasi-optimality, at the expense of a poor scaling with respect to the failure probability ϵ. The proof is based on Markov inequality, which, roughly speaking, is used to bound the discrete error term arising in LSerrbd.


[First nonuniform error bound in probability]

Let 0 < δ,ϵ < 1, f ∈ L^2_ϱ(D), ⊂ L^2_ϱ(D) be a finite-dimensional subspace with () = n and μ_1,…,μ_m be probability measures satisfying Assumption <ref>. Consider sample points drawn randomly and independently with x_i ∼μ_i, i = 1,…,m, where 


m ≥c_δ ·κ_w() ·log(4 n /ϵ),    c_δ = ((1+δ) log(1+δ) - δ))^-1,

and w = 1/ν and κ_w are as in mu_weight_fn and kappa-w-def, respectively. Then the following hold with probability at least 1-ϵ. For any e ∈^m, the weighted least-squares approximation f̂ is unique and satisfies


f - f̂_L^2_ϱ(D) ≤( 1 + √(2 κ_w()/m ϵ) 1/1-δ )inf_p ∈ f - p_L^2_ϱ(D) + 1/√(1-δ) e_2,w.

Moreover, if {ϕ_i }^n_i=1 is an orthonormal basis of , then cond(A) ≤√(1+δ/1-δ).


Suppose again that δ = 1/2 and w(x) ≲ 1, ∀ x. Then this bound implies that

f - f̂_L^2_ϱ(D) ≲( 1 + 1/√(ϵlog(4 n / ϵ) ) ) inf_p ∈ f - p_L^2_ϱ(D) + e_2.

While stable and L^2_ϱ-quasi-optimal, the scaling with respect to ϵ is undesirable. To obtain an ϵ-independent bound, this suggests we either need to n to be exponentially large in 1/ϵ, or impose an additional constraint on m that m ≳κ_w() / ϵ. 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 LSerrbd is a sum of independent random variables with bounded variance. It therefore concentrates exponentially fast (in m) around its mean f - p^2_𝖽𝗂𝗌𝖼,w = f-p^2_L^2_ϱ(D) (recall exp-sum-scaling). This leads to the following result, which is also nonuniform.


[Second nonuniform error bound in probability]

Consider the setup of Corollary <ref> with m-cond-one replaced by 


m ≥c_δ ·κ_w() ·log(4n /ϵ)  and  m ≥2 ·k ·log(4/ϵ)

for some k > 0. Then the following hold with probability at least 1-ϵ. For any e ∈^m, the weighted least-squares approximation f̂ is unique and satisfies


f - f̂_L^2_ϱ(D) ≤( 1 + √(2/1-δ)  ) inf_p ∈ { f - p_L^2_ϱ(D) + √(w)(f-p)_L^∞_ϱ(D)/√(k) }  + 1/√(1-δ) e_2,w.

Moreover, if {ϕ_i }^n_i=1 is an orthonormal basis of , then cond(A) ≤√(1+δ/1-δ).


This result asserts a mixed type of quasi-optimality, involving the L^2_ϱ-norm and a (weighted) L^∞-norm divided by the factor √(k). Notice that the factor √(w) can be removed whenever w(x) ≲ 1, a.e. x ∼ϱ, as will be the case later. Therefore, consider, as in Remark <ref>, a case where the L^2_ϱ-norm best approximation error decays algebraically fast in n = ().
As we noted therein, the L^∞_ϱ-norm best approximation error often decays √(n) slower than the former. Hence, one may choose k = n in m-conds-in-prob-2 to show that f̂ achieves the same algebraic convergence rate in L^2_ϱ-norm as the L^2_ϱ-norm best approximation in . This approach has also been used in the related context of function approximation via compressed sensing in <cit.> and <cit.>.






 §.§ Error bounds in expectation


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 <cit.>, which is based on <cit.>. 

Let {ϕ_i }^n_i=1 be an orthonormal basis of . First, we notice that alpha-beta-delta holds whenever

G - I_2 ≤δ,

where G is the discrete Gram matrix

G = ( ϕ_jϕ_k_𝖽𝗂𝗌𝖼,w )^n_j,k=1 ∈^n ×n.

It is then possible to show the following bound.



Let 0 < δ < 1, f ∈ L^2_ϱ(D), ⊂ L^2_ϱ(D) be a finite-dimensional subspace with dim() = n and μ_1,…,μ_m be probability measures satisfying Assumption <ref>. Consider sample points drawn randomly and independently with x_i ∼μ_i, i = 1,…,m, where

m ≥c_δ ·κ_w() ·log(2 n /ϵ),     c_δ = ((1+δ) log(1+δ) - δ))^-1,

and w = 1/ν and κ_w are as in mu_weight_fn and kappa-w-def, respectively.
 Then

( f - f̂^2_L^2_ϱ(D) χ_G-I_2 ≤δ ) ≤( 1 + 2/(1-δ)^2  κ_w()/m ) inf_p ∈ f - p^2_L^2_ϱ(D) + 2/1-δ e^2_2,w,

where χ_E denotes the indicator function of an event E.


This lemma can be used to construct two estimators with error bounds in expectation. The first is the conditioned estimator <cit.>, which is defined as

f̂^𝖼𝖾 = f̂ χ_G - I_2 ≤δ.

Computing this estimator requires one to evaluate

G-I_2 = max{ | 1 - σ^2_max(A) | , | 1 - σ^2_min(A) | }  = max{ | 1 - α^2_w | , | 1- β^2_w | }.

Having done this, one simply sets f̂^𝖼𝖾 = f̂ when G - I_2 ≤δ and f̂^𝖼𝖾 = 0 otherwise. 

The conditioned estimator has the disadvantage that it requires an orthonormal basis for  to be known – a property that may not hold in practice (see <ref>). This can be avoided by using a truncated estimator. This approach assumes an a priori bound for f of the form

f_L^2_ϱ(D) ≤σ,

for some known σ≥ 0. Now define the truncation operator _σ : L^2_ϱ(D) → L^2_ϱ(D) by

    _σ(g) 
    := min{1, L/g_L_ϱ^2()} g
    =
    
    g,    g_L^2_ϱ()≤σ, 
    σ g / g_L_ϱ^2(),    g_L^2_ϱ() > σ,
      ∀ g ∈ L^2_ϱ().

Then this estimator is defined as

f̂^𝗍𝖾 = _σ(f̂).

Note that one can also construct a truncated estimator with respect to the other norms. The L^∞_ϱ-norm has also been commonly used for this purpose <cit.>.



[Nonuniform error bound in expectation]

Let 0 < δ,ϵ < 1, f ∈ L^2_ϱ(D), ⊂ L^2_ϱ(D) be a finite-dimensional subspace with dim() = n and μ_1,…,μ_m be probability measures satisfying Assumption <ref>. Consider sample points drawn randomly and independently with x_i ∼μ_i, i = 1,…,m, where


m ≥c_δ ·κ_w() ·log(2n/ϵ),     c_δ = ((1+δ) log(1+δ) - δ))^-1,

and w = 1/ν and κ_w are as in mu_weight_fn and kappa-w-def, respectively.
Then

f - f̂^𝖼𝖾^2_L^2_ϱ(D)  ≤( 1 + 2/(1-δ)^2  κ_w()/m ) inf_p ∈ f - p^2_L^2_ϱ(D) + 2/1-δ e^2_2,w + f^2_L^2_ϱ(D) ϵ.

The same bound holds for the f̂^𝗍𝖾, except with the final term replaced by 4 σ^2 ϵ.


Observe that the factor

2/(1-δ)^2 κ_w()/m ≤2/(1-δ)^2 c_δ 1/log(2n) ≤2/(1-δ)^2 c_δ 1/log(2n/ϵ) →0

as n →∞. Hence, this bound asserts L^2_ϱ-quasi-optimality of the two estimators (with constant approaching 1) up to the ϵ term. 



[Removing the ϵ term]

As discussed in <cit.>, the additional ϵ term in the error bound may cause problems when striving to achieve rate-optimal bounds in the absence of noise. In particular, if the best approximation error inf_p ∈f-p_L^2_ϱ(D)≲ρ^-n decays exponentially fast (or faster) in n = () for some ρ > 1 then achieving the same rate for the least-squares approximation would require setting ϵ = ρ^-n. This adds an additional multiplicative factor of n in the sample complexity bound sampl-com-exp-err, thus prohibiting rate optimality (recall that κ_w() ≥ n). One way to remove this term, which was introduced in <cit.>, is to repeatedly redraw the sample points {x_1,…,x_m} until the condition G - I_2 ≤δ is met, and then use the resulting points to construct the weighted-least squares estimator f̂^⋆. By doing this, one can achieve a similar error bound in expectation, except without the term ϵ (the term ϵ now only influences the expected number of redraws needed to achieve G - I_2 ≤δ). See <cit.> for further information.





[The noise bound]

Due to Assumption <ref>, the noise term in Theorem <ref> satisfies, in general, the bound

e^2_2,w = 1/m ∑^m_i=1 w(x_i) | e_i|^2 = 1/m ∑^m_i=1 ∫_D w(x) μ_i(x) | e_i|^2 ≤e^2_∞ ∫_D w(x) μ(x) = e^2_∞.

If μ_1 = ⋯ = μ_m = μ, then one also has e^2_2,w = 1/me^2_2.
Moreover, as in Remark <ref>, this can also be achieved (up to a constant) in the general case if w(x) ≲ 1, a.e. x ∼ϱ. 













§ CHRISTOFFEL SAMPLING


We now come to the crux of this article, which is devise random sampling schemes that achieve near-optimal sample complexity bounds.







 §.§ Optimal choice of weight function via the Christoffel function


The results shown in the previous section relate the number of measurements m to the constant κ_w(). Hence, our goal is to choose a weight function w that minimizes this constant. Recall that  

κ_w() = _x ∈ϱ w(x) ()(x).

A natural first choice involves selecting 

w(x) ∝1/()(x).

Applying the normalization condition w_normalization (recall that ν = 1/w) and the fact that ∫_D()(x) ϱ(x) = n (recall Kappa-def-alt), we obtain


w(x) = n/()(x).

This choice is quite popular in the literature. 
However, it requires the additional assumption that ()(x) > 0 almost everywhere. This is a rather mild assumption, which is equivalent to requiring that for almost every x ∼ϱ there exists a p = p_x ∈ for which p_x(x) ≠ 0 (in particular, it holds under the assumption made in Remark <ref>). If this holds, then w-opt-orig is the optimal for choice of w, since κ_w() = n achieves the optimal lower bound (recall kappa-w-lb). 

However, this choice may not be desirable for the reasons considered in  Remark <ref>. If ()(x) is small at x, then w(x) becomes large and potentially blows up the noise. Fortunately, this issue can be resolved and the above assumption removed, by slightly changing w(x) to


w(x) = ( 1/2 + 1/2 ()(x)/n )^-1.

This avoids the assumption on () and leads to a bounded weight function satisfying w(x) ≤ 2. The only cost is suboptimality by a factor of 2, i.e., κ_w() ≤ 2n.
The reader will likely notice that one could replace the 1/2 in w-opt with a weighted combination θ + (1-θ) ()(x)/n for any 0 < θ < 1, giving κ_w() ≤ (1-θ)^-1 n and w(x) ≤θ^-1. For simplicity, we consider the factor 1/2 throughout. 

This aside, having chosen w as in w-opt-orig (one could also consider w-opt), to achieve near-optimal sampling we wish to select measures μ_1,…,μ_m such that mu_weight_fn holds, i.e.,


1/m ∑^m_i=1 μ_i(x) = ()(x)/n  ϱ(x)=∑^n_i=1 | ϕ_i(x) |^2 /n  ϱ(x).

If the measures are chosen so that this holds, then the various sample complexity estimates of the previous section are near-optimal in n. Indeed, letting δ = 1/2, we see that the condition

m ≳n ·log(2n/ϵ)

suffices to ensure the various `in probability' or `in expectation' bounds given in the previous section.




 §.§ Christoffel sampling


There are several ways to achieve mu_weight_fn_opts. Arguably the simplest involves setting


μ_1 = ⋯= μ_m = μ,    where μ(x) =   ()(x)/n ϱ(x) = ∑^n_i=1 | ϕ_i(x) |^2 /n ϱ(x).

However, this strategy is not well suited in the case of hierarchical approximation (Remark <ref>). Indeed, let μ^(k) be the optimal sampling measure for ^(k), as defined in mu-opt-1 with = ^(k). Now suppose that the first m_1 points x^(1)_i ∼_i.i.d.μ^(1), i = 1,…,m_1, where μ^(1) is as in mu-opt-1 for = ^(1). Then we would like to recycle these m_1 points { x^(1)_i }^m_1_i=1 when constructing the second set of sample points { x^(2)_i }^m_2_i=1. However, since μ^(1)≠μ^(2) in general, these m_1 points are drawn with respect to the wrong measure for near-optimal sampling in the subspace ^(2). Thus, it is not possible to achieve near-optimal sampling simply by augmenting the set { x^(1)_i }^m_1_i=1 with m_2 - m_1 new points.
 
One strategy to overcome this limitation involves interpreting μ^(k+1) as an additive mixture involving μ^(k) and a certain update measure σ^(k). One can then use this construct a sampling procedure that recycles `most' of the first m_1 points, while ensuring that the overall sample is drawn i.i.d. from μ^(k+1) <cit.>.
An alternative approach, introduced in <cit.>, involves choosing measures μ_i according to the individual basis functions. Let {ϕ_i }^n_i=1 be an orthonormal basis of  and suppose that m = k n for some k ∈. Then we define


μ_i(x) =  | ϕ_j(x) |^2 ϱ(x),  (j-1) k < i ≤j k, j = 1,…,n.

Observe that

1/m ∑^m_i=1 μ_i(x) =  k/m ∑^n_j=1 | ϕ_j(x) |^2 ϱ(x) = ()(x)/n ϱ(x),

due to Kappa-def-alt. Hence mu_weight_fn_opts also holds for this choice, guaranteeing quasi-optimal sampling. Moreover, as each sampling measure corresponds to a single basis function, rather than an additive mixture of basis functions as in mu-opt-1, this approach readily lends itself to hierarchical approximation. See <cit.> for further details. Note that the distributions corresponding to these measures known as induced distributions <cit.>, as they are induced by the orthonormal basis {ϕ_i }^n_i=1. 

We will, henceforth, refer to either procedure – or, indeed, any selection of measures μ_i for which mu_weight_fn holds for w = 1/ν as in w-opt or w-opt-orig – as Christoffel sampling.




 §.§ Uniform error bounds in probability




To conclude this section, we now describe how a further modification of the near-optimal sampling measure can lead to uniform bounds in probability that improve on the somewhat crude bounds shown in Corollary <ref> and achieve something close to L^2_ϱ-quasi-optimality. This section is based on techniques developed in <cit.> to estimate sampling numbers. See also <ref> for additional discussion.


For this, we assume that there is an orthonormal basis {ϕ_i }^∞_i=1⊂ L^2_ϱ(D) and that


= _n = { ϕ_1,…,ϕ_n}.

For convenience, given f ∈ L^2_ϱ(D) let


e_n(f) = inf_p ∈_n f - p_L^2_ϱ = √(∑^∞_i > n | c_i |^2),

where c_i = fϕ_i_L^2_ϱ(D) is the ith coefficient of f. The second equality is due to Parseval's identity.
We now construct the measure. Define sets

I_l = { n 2^l + 1,…, n 2^l+1 },  l = 0,1,2,…
and consider a sequence (v_l)^∞_l=0 with ∑^∞_l=0 v^2_l = 1. Then let


μ(x) = ( 1/2 + 1/4 ∑^n_i=1 | ϕ_i (x) |^2/n + 1/4 ∑^∞_l=0 v^2_l/|I_l| ∑_i ∈I_l | ϕ_i(x) |^2 ) ϱ(x).





Let 0 < δ,ϵ < 1, n ∈, {ϕ_i }^∞_i=1⊂ L^2_ϱ(D) be an orthonormal basis and = _n be as in Pn-def-unif. Let 0 < p < 2 and v_l = c_θ 2^-θ l for 0 < θ < 1/p-1/2, where c_θ is such that ∑^∞_l=0 v^2_l = 1. Consider sample points x_i ∼_i.i.d.μ, i=1,…,m, where μ is as in opt-meas-krieg and


m ≥4 ·c_δ ·n ·log(4n/ϵ),    c_δ = ((1+δ) log(1+δ)-δ))^-1.

Then the following holds with probability at least 1-ϵ. For any f that is defined everywhere in D and for which (e_n(f) )^∞_n=1∈ℓ^p() and any noise e ∈^m, the weighted least-squares approximation f̂ is unique and satisfies


 f - f̂_L^2_ϱ(D) ≤c_p,θ/√(1-δ) ( 1/n ∑_k ≥n (e_k(f))^p )^1/p  + √(2/1-δ) e_2,

where c_p,θ > 0 is a constant depending on p and θ only. Moreover, the condition number of the least-squares matrix ls-Ab satisfies cond(A) ≤√(1+δ/1-δ).


To put this result into context, consider a class of functions  that are defined everywhere on D and for which 

sup_f ∈ e_n(f) ≍n^-α log^β(n+1)
 
for some α > 1/2 and β∈. This holds, for instance, in the case of polynomial approximation when  is a unit ball of functions of finite regularity. Then the right-hand side of f-hatf-bound satisfies

( 1/n ∑_k ≥n (e_k(f))^p )^1/p ≤c_α,β,p n^-α log^β(n+1).

Hence, with probability at least 1-ϵ, one obtains a matching error bound for the least-squares estimator (up to constants), uniformly for functions f ∈, with near-optimal sample complexity.








§ IMPROVEMENT OVER MONTE CARLO SAMPLING


Having derived near-optimal sampling strategies, in this section we discuss how this compares against standard Monte Carlo sampling, i.e., i.i.d. random sampling from the measure ϱ. Recall that in this case, the weight function w is precisely w ≡ 1, meaning that Monte Carlo sampling corresponds to an unweighted least-squares approximation. Theorem <ref> asserts that the sample complexity of Monte Carlo sampling therefore depends on the unweighted quantity


κ() = κ_1() = ()_L^∞_ϱ(D).

In other words, the maximal behaviour of the Christoffel function (), or equivalently, the optimal constant in the unweighted Nikolskii-type inequality nikolskii.




 §.§ Bounded orthonormal systems




Recall that κ_w() ≥ n for any w, and therefore, κ() ≥ n. On the other hand, if κ() ≤ c n for some c ≥ 1, then Monte Carlo sampling is already a near-optimal strategy, and there may be little need to optimize the sampling measure (besides reducing the constant c).
This situation occurs whenever  has an orthonormal basis {ϕ_i }^n_i=1 that is uniformly bounded. Such a basis is sometimes referred to as a bounded orthonormal system (see, e.g., <cit.> or <cit.>). Specifically, if


ϕ_i^2_L^∞_ϱ(D) ≤c,  ∀i = 1,…,n,

then it follows immediately from Kappa-def-alt that κ() ≤ c n. Subspaces of trigonometric polynomials are a standard example for which this property holds (with c = 1). Closely related are the Chebyshev polynomials of the first kind (see <ref>). The orthonormal Chebyshev polynomials on (-1,1) are defined by

ψ_0(x) = 1,  ψ_i(x) = √(2) cos(i arccos(x)), i ∈.

Therefore they satisfy BOS-c c = 2. In d dimensions, the tensor-product Chebyshev polynomial basis on (-1,1)^d satisfies BOS-c with c = 2^d. Hence it is also a bounded orthonormal system, albeit with a constant that grows exponentially fast as d →∞. 



 §.§ Arbitrarily-bad sample complexity bounds


Unfortunately, bounded orthonormal bases are quite rare in practice. It is also straightforward to construct subspaces  for which κ() can be arbitrarily large in comparison to n. One way to do this involves the Legendre polynomials. Unlike the Chebyshev polynomials, the univariate Legendre polynomials grow with their degree, and satisfy (see, e.g., <cit.>)


ψ_i _L^∞(-1,1) = | ψ_i(±1) | = √(i+1/2).

Therefore, for any set S ⊂_0, the space _S = {ψ_i : i ∈ S } has constant κ(_S) given by κ(_S) = ∑_i ∈ S (i+1/2) (this follows from Kappa-def-alt and MC-kappa). The following result is now immediate.




There exists a probability space (D,,ϱ) such that the following holds. For every n ∈ and C > 0, there exists a subspace ⊂ L^2_ϱ(D) of dimension n such that κ() ≥ C.


The reason for this poor behaviour is that Legendre polynomials are highly peaked near the endpoints x = ± 1. Specifically, the ith such polynomial has absolute value √(i+1/2) at x = ± 1, but is 1 uniformly within compact subsets of (-1,1). This points towards a general observation. We expect poor scaling of κ(), and therefore poor performance of Monte Carlo sampling, whenever  has an orthonormal basis that is highly localized around a common shared point. 



 §.§ Sample complexity bounds for polynomial spaces


Following <ref>, we now describe a series of further bounds for κ() in the case of multivariate polynomial spaces = _S, focusing on the case where S is a lower set.


Chebyshev polynomials. As noted previously, any subspace

_S = { Ψ_ν : ν∈S },    S ⊂^d_0, |S| = n,

of multivariate Chebyshev polynomials satisfies the exponentially-large (in d) bound κ(_S) ≤ 2^d n. Fortunately, if S is a lower set, then one has the d-independent bound (see, e.g., <cit.>)

κ(_S) = κ(_S) ≤n^log(3)/log(2),

This bound is sharp up to a constant, i.e., there is a lower set with κ(_S) ≳ n^log(3)/log(2). For d ≥log_2(n/2), this bound is sharper than the previous bounds  κ(_S) ≤ 2^d n. 

Improvements such as this are typical when lower set structure is imposed – we will see another instance of this in a moment. It is one of the features that makes lower sets desirable for high-dimensional approximation. The underlying reason for it is because a lower set cannot contain too many `large' (or even nonzero) indices. 

Legendre polynomials. Because of Kappa-def-alt, MC-kappa and leg-poly-unif-bd, for any subspace _S = {Ψ_ν : ν∈ S } of Legendre polynomials, one has

κ(_S) = ∑_ν∈S ∏^d_k=1 (ν_k+1),

which can be arbitrarily large in comparison to n = |S|. However, lower set structure leads to a dramatic improvement over the general case. If S is lower, then the following sharp upper bound holds (see, e.g., <cit.>)


κ(_S) ≤n^2.


Ultraspherical and Jacobi polynomials. The situation for Jacobi polynomials with max{α,β} > -1/2 is similar to that of Legendre polynomials. The quantity κ(_S) can be made arbitrarily large for general S. However, if S is lower then one has the finite bound


κ(_S) ≤n^2max{ α, β} + 2 ,    when α,β∈_0

for Jacobi polynomials (see <cit.>) and


κ(_S) ≤n^2 α+ 2 ,    when 2α+1 ∈,

for ultraspherical polynomials  (see <cit.>).



[Sharpness of the rates]
The bounds kappa-legendre–kappa-ultraspherical imply that a superlinear sample complexity suffices for stable and accurate polynomial approximation with random samples drawn from Jacobi measures. These rates are also necessary. This was recently shown in <cit.> in the d = 1 case, based on earlier work <cit.> on deterministic points that are equidistributed with respect to such measures. Specifically, <cit.> shows that choosing m ≍ n^τ, where τ < 2(max{α,β}+1), necessarily implies exponential instability of the least-squares approximation (or, indeed, any other approximation method that achieves similar accuracy).



Bounded, non-tensor product domains. Several of these bounds generalize to bounded non-tensor product domains <cit.>. If ϱ is the uniform measure and D ⊂^d is bounded and satisfies the inner cone condition (see, e.g., <cit.>), then

κ(_S) ≤c_D ·n^2  if S = S^𝖳𝖣_p and n = | S^𝖳𝖣_p|.

See <cit.>.
Thus, the same quadratic growth occurs for the total degree index set. Further, if D has C^2 boundary, then <cit.>

κ(_S) ≤c_D ·n^d+1/d  if S = S^𝖳𝖣_p and n = | S^𝖳𝖣_p|.

See <cit.> and references therein for further results in this direction. These bounds apply only to the total degree index set TP-TD. For arbitrary lower sets, one has (see <cit.>)

κ(_S) ≤n^2 /λ  if S is lower and n = |S|,

whenever D has the λ-rectangle property: there is a 0 < λ < 1 such that for any x ∈ D there exists an axis-aligned rectangle R ⊆ D containing x for which |R| ≥λ |D|.


Hermite and Laguerre polynomials. Unfortunately, this analysis of Monte Carlo sampling says nothing about Hermite and Laguerre polynomial approximations, for the simple reason that such polynomials are not uniformly bounded, and therefore the constant MC-kappa satisfies κ() = +∞. The sample complexity of Hermite and Laguerre polynomial approximation is poorly understood in comparison to that of Jacobi polynomials. A number of empirical studies have suggested a super-polynomial or exponential sample complexity (in n, for fixed d). But relatively few theoretical estimates exist. See <cit.> and <cit.>. Suffice to say, though, Hermite and Laguerre polynomial approximations are examples where one stands to gain substantially from Christoffel sampling, and as such, these has often been used as examples to illustrate its efficacy <cit.>.





 §.§ Summary


We summarize this section as follows. In general, the sample complexity of Monte Carlo sampling depends on the L^∞-norm of the Christoffel function, and it is easy to construct examples where this is arbitrarily large in comparison to n = (). Moreover, these scenarios arise in various polynomial approximation problems, especially when considering unbounded domains.





§ VARIATIONS, EXTENSIONS AND FURTHER TOPICS


To conclude this first part of the article, we now discuss a number of issues not considered so far, along with some variations and extensions.








  §.§.§ Sampling from the near-optimal measures


A key practical issue is drawing samples from the measure mu-opt-1 or mu-opt-2 in a computationally efficient manner. As observed in <cit.>, mu-opt-1 is an additive mixture of measures of the form mu-opt-2. Hence it is enough to be able to draw samples from the induced distributions | ϕ_i(x) |^2 ϱ(x), where {ϕ_i }^n_i=1 is any orthonormal basis for . 
A prerequisite for this is, of course, the existence of an explicit orthonormal basis. This is often the case – for example, the problem considered in <ref>, which involve tensor products of orthogonal polynomials, for which fast algorithms for sampling from the induced distributions exist <cit.> – but not always. One such case is polynomial approximation on nontensorial domains. Orthogonal polynomials can be defined explicitly for certain non-tensorial domains, e.g., simplices, balls and spheres <cit.>. Yet, this is impossible in general. 

A simple way to address this problem is to replace D by a fine grid Z = { z_i }^K_i=1 of points drawn i.i.d. from ϱ and replace the continuous measure ϱ by the discrete uniform measure ϱ̅ = K^-1∑^K_i=1δ_z_i supported on Z <cit.>. Given a nonorthogonal basis {ψ_i }^n_i=1 for , now considered a subspace of L^2_ϱ̅(D), one may construct an orthonormal basis via numerical linear algebra tools such as classical QR factorizations or SVDs, or through more recent approaches such as V+A (Vandermonde with Arnoldi) <cit.>. Sampling from the Christoffel sampling measure is then straightforward, since it is also a discrete measure supported on Z.

This approach, which is a form of leverage score sampling (see next), ensures accurate and stable recovery in the L^2_ϱ̅(D)-norm. To guarantee the same in the original L^2_ϱ(D)-norm one needs to that ensure K is large enough. Since Z is obtained by Monte Carlo sampling, this is governed, as in <ref> and specifically MC-kappa, by the constant MC-kappa. This is another reason why analyzing the maximal behaviour of the Christoffel function is important, since it can inform the size of grid needed in such a construction. See <cit.> for further details.  

Note that the computation of an orthonormal basis over Z is a purely offline cost. It does not involve additional evaluations of the target function f, which are often the main computational bottleneck in practice.  Of course, theoretical bounds for the Christoffel may not be available or, if available, may result value of K that is too large for practical computations. This has spurred several recent works <cit.> which develop more refined strategies for constructing the grid K than Monte Carlo sampling.

An alternative to this approach involves using a structured (e.g. Quasi-Monte Carlo) grid for Z <cit.>. This has the additional benefit of ensuring that the algebraic least-squares problem can be solved efficiently via Fast Fourier Transforms (FFT).







  §.§.§ Connections to matrix sketching and leverage score sampling


Let X ∈^N × n, where N ≥ n, and x ∈^N. In applications in data science, it may be infeasible to solve the `full' least squares problem w ∈z ∈^nX z - x_2. Matrix sketching involves using a sketching matrix S ∈^m × N (a matrix with one nonzero per row) and solving the sketched problem ŵ∈z ∈^nS X z - S x_2. The objective is to find S with a minimal number of rows m such that

X ŵ - x_ℓ^2 ≲X w - x_2.

In random sketching, one considers a discrete probability distribution π = {π_1,…,π_N } on {1,…,N} with π_i > 0, ∀ i, draws j_1,…,j_m ∼_i.i.d.π and then sets S_i,j_i = 1/√(π_j_i) and S_ij = 0 otherwise. A near-optimal choice of distribution π is provided by the leverage scores ℓ_i(X), i = 1,…,N, of the matrix X. These are given by


ℓ_i(X) = max_z ∈^n 
 X z ≠0 |(X z)_i|^2/X z^2_2 ≡Q(i,:)^2_2, 

where Q ∈^N × n is any matrix whose columns form an orthonormal basis for Ran(X).
The resulting procedure is the well-known technique of leverage score sampling <cit.>.

Leverage score sampling can be considered as a special case of Christoffel sampling involving the discrete set D = {1,…,N}. Note that any vector x = (x_i)^N_i=1∈^N can be viewed as a function f ∈ L^2_ϱ(D) via the relation f(i) = x_i and vice versa. One now defines the subspace = { X z : z ∈^n } to cast the above problem into the form introduced in <ref>. In particular, the values of the Christoffel function () over the set D are precisely the leverage scores leverage-scores of the matrix X. We refer to <cit.> and <cit.> for details.










  §.§.§ n sampling strategies, frame subsampling and Kadison–Singer


The Christoffel sampling schemes described in this paper are near-optimal, in the sense the sample complexity is log-linear in n. Due to the coupon collector's problem, this is the best achievable when using i.i.d. samples (see, e.g., <cit.>).

This limitation has spurred a recent line of work on methods that have linear sample complexity in n, i.e., optimal up to possible constants, and even schemes that can achieve interpolation, i.e., m = n. This is based on techniques introduced in <cit.> on frame subsampling. Here, given a frame of m ≥ n vectors in ^n, one asks whether it is possible to extract a subset of size n that, after potential reweighting, still constitutes a frame. This problem is closely related to Weaver's KS_2 conjecture, see <cit.> or <cit.>). Weaver's conjecture (now theorem) is equivalent to the Kadison–Singer problem, and form the basis of the proof in <cit.> for the latter.
The connection to sampling for least-squares approximation comes from the Marcinkiewicz–Zygmund inequality MZ-inequality, which can be recast as a frame condition for the vectors x^(1),…,x^(n) defined by x^(i) = ( √(w(x_i))ϕ_j(x_i)  )^n_j=1, i = 1,…,n.

The work of <cit.> (as well as extensions due to <cit.>) has been used to show the existence of sampling points with n sample complexity for an arbitrary subspace . See <cit.> and references therein. Unfortunately, these works are impractical, as the computational cost in constructing the subsample is (at least) exponential in n.
Fortunately, recent progress has been made by using the approach of <cit.>, as well as techniques from <cit.>, leading to practical algorithms that run in polynomial time. See <cit.> and <cit.> for two such approaches, as well as <cit.> for related work in the discrete setting. A significant result of <cit.> is to provides polynomial-time algorithms that also work down to the interpolation regime m = n = (), albeit with constants in the error bounds that grow algebraically with n.





  §.§.§ Sampling numbers and information-based complexity


Another major area of focus in the last several years has been the use of (weighted) least-squares approximation, Christoffel sampling and subsampling techniques to provide new results in information-based complexity <cit.>. In this line of work, one considers a compact subset  of a Banach space, then studies object such as the (linear) sampling numbers for . These measure how well (in a worst-case sense over ) one can approximate functions in  from m samples using arbitrary (linear) reconstruction maps. New results have emerged that bound the sampling numbers (for arbitrary ) in terms of its Kolmogorov n-width, i.e., the best approximation error (uniformly over ) that can be achieved by any n-dimensional subspace. These results show that pointwise samples (known as standard information) can constitute near-optimal information for recovery. Some of the core ideas of this work can be found in Theorem <ref>, including the construction of the measure opt-meas-krieg which is due to  <cit.>.
Note that sampling number are often formulated with respect to the L^2-norm (as in this article), but recent works also consider other L^p-norms – in particular, the uniform norm. For a selection of the many recent results in this direction, see <cit.> and references therein.




  §.§.§ Sampling discretizations


In tandem with these efforts, there has also been a focus on the development and systematic study of sampling discretizations using these ideas, both in the L^2-norm such as in MZ-inequality and in other L^p-norms. We refer to <cit.> for reviews, as well as <cit.>, and references therein. Note that L^∞-norm sampling discretizations are related to the construction of weakly admissible meshes. See <cit.> for recent work on weakly admissible meshes that employs similar techniques.

As we have seen, sampling discretization are sufficient conditions for accurate and stable recovery via (weighted) least squares. However, they are also necessary conditions for stable recovery by any method. Modifying <cit.>, let : ^m → L^2_ϱ(D) be an arbitrary reconstruction map and suppose that  is δ-accurate over , i.e.,


p - ({ p(x_i) }^m_i=1)_L^2_ϱ(D) ≤δp_L^2_ϱ(D),  ∀p ∈.

Note that his holds with δ = 0 for least squares whenever α_w > 0, due to Lemma <ref>. Now let  be any set of functions that are defined at the { x_i }^m_i=1 and suppose that ⊆. Then it is a short argument to show that the ϵ-Lipschitz constant

L_ϵ(; ) = sup_f ∈ sup_0 < e_2,w ≤ϵ ({f(x_i)}^m_i=1) + e) - ({ f(x_i) }^m_i=1) _L^2_ϱ(D)

of the map  restricted to the domain { (f(x_i))^m_i=1 : f ∈}⊆^m satisfies

L_ϵ(; )  ≥(1-δ)/√(α_w),

where α_w is the lower constant in the sampling discretization MZ-inequality. It follows that a reconstruction map cannot be both accurate (even over , as in R-accurate) and stable without a sampling discretization.


















  §.§.§ Alternative sampling methods


Many other sampling methods have been proposed over the last decade, especially in the context of high-dimensional polynomial approximation. However, these are generally lacking near-optimal sample complexity guarantees. See <cit.> and <cit.> and references therein for an overview of the many different approaches considered.

A limitation of Christoffel sampling is that i.i.d. points may cluster, thereby reducing the practical efficiency of the scheme. Most recently, a number of works have explored ideas such as volume sampling <cit.> using determinantal point processes to overcome this limitation. These are well-known concepts in machine learning <cit.>, in which non-independent samples are drawn from a measure that promotes repulsion between the points. This transpires to be closely related to Christoffel sampling, since the marginals of the sample points follow the same distribution. The application of volume sampling to least-squares approximation in arbitrary subspaces has been considered in <cit.> for reproducing kernel Hilbert spaces and <cit.> for general spaces, along with its theoretical analysis and comparison with Christoffel sampling. Despite practical benefits, however, it is as of yet unclear whether they offer theoretical advantages over Christoffel sampling <cit.>.






  §.§.§ Adaptive methods


Finally, we briefly mention the prospect of adaptive methods. While these methods typically lack full theoretical guarantees, they can prove extremely effective in practice. In a variation of Remark <ref>, in an adaptive scheme one also chooses each subspace ^(k+1) adaptively based on the previous approximation f̂^(k). In this case, we term this procedure an adaptive approximation scheme. This can be done using greedy methods <cit.>, as in <cit.> (which are themselves based on adaptive quadrature routines <cit.>), when given a dictionary of candidate basis functions to use to build the spaces ^(k). This aside, adaptive methods can also be used when constructing approximations in complex, nonlinear approximation spaces. See <ref> for further discussion.  














§ BEYOND LINEAR SPACES AND POINTWISE SAMPLES


Up to now, we have considered approximating an unknown function f : D → from a collection of m pointwise samples in an n-dimensional approximation space ⊂ L^2_ϱ(D). In this final section, we introduce a general framework that significantly extends this setup. This section is primarily based on the framework introduced in <cit.>, which was then further extended in <cit.>. Unlike in previous sections, our presentation will now be less thorough: we aim to convey the main ideas without the full details or variations. See <cit.> for in-depth treatments, and <cit.> for related work.

There are four main extensions we now address. First, the target object f need not be a scalar-valued function, but simply an element of a Hilbert space . Second, the measurements arise as evaluations of arbitrary linear operators, which may be scalar-, vector- or function space-valued. Third, there may be C ≥ 1 different distinct processes generating the measurements. And finally, the approximation space  need not be finite-dimensional subspace of . Examples that motivate these generalizations are discussed in <ref>.






 §.§ The general framework


Let (Ω,,) be a probability space,  be a separable Hilbert space and consider a normed vector subspace of _0 ⊆, termed the object space. Let C ≥ 1 be the number of measurement processes. For each c = 1,…,C, let (D_c,_c,ϱ_c) be a probability space, which we term the measurement domain, _c be a Hilbert space, which we term the measurement space, and

L_c : D_c →(_0,_c),

be a mapping from D_c to the space of (_0,_c) of bounded linear operators _0 →_c, which we term the sampling operator.

 For each c, let μ_c be such that (D_c,_c,μ_c) is a probability space. We assume that μ_c is absolutely continuous with respect to ϱ_c and its Radon–Nikodym derivative ν_c : D_c → is strictly positive almost everywhere on supp(ϱ_c). By definition

∫_D_c ν_c(θ) ϱ_c(θ) = 1.

Note that, in what follows, we use θ or θ_c to denote the variable in D_c, rather than x. For convenience, we also define the weight function w_c : D_c → by w_c(θ) = 1/ν_c(θ), θ∈ D_c.

Let m_1,…,m_C ∈, where m_c is the number of measurements in the cth measurement process.  
With this setup in hand, we now draw samples independently with θ_ic∼μ_c, i = 1,…,m_c, c = 1,…,C and consider the measurements


y_ic = L_c(θ_ic)(f) + e_ic ∈_c,  i =1,…,m_c, c = 1,…,C,

where e_ic∈_c is a noise term.
Finally, we let ⊂ be the approximation space (which now need not be a linear space) and we consider the empirical least-squares fit


f̂ ∈p ∈ ∑^C_c=1 1/m_c ∑^m_c_i=1 w_c(θ_ic) y_ic - L_c(θ_ic)(f) ^2__c.

Note that computing solutions to general-least-squares may be very challenging numerically when  is a nonlinear set. However, this is highly dependent on the choice of , and therefore not the focus of this article.







 §.§ Examples


As shown in <cit.>, this framework includes many problems of practical relevance. We now summarize several such examples. We start by showing that it generalizes the setup of <ref>.

(i) Scalar-valued function approximation from pointwise samples.
The setup of  <ref> be considered in the general framework as follows. Let C = 1, = L^2_ϱ(D), _0 = C(D), D_c = D, ϱ_1 = ϱ and _1 = (with the Euclidean inner product). We then define the sampling operator L_1 = L as the pointwise sampling operator L(x)(f) = f(x) for x ∈ D and f ∈_0. Note that the measurements general-meas and least-squares fit general-least-squares reduce to f_meas and wls-prob, respectively


(ii) Function approximation from gradient-augmented samples.
A simple modification of (i) involves sampling both the function and its gradient. This arises in various applications, including parametric PDEs and UQ <cit.>, seismology, Physics-Informed Neural Networks (PINNs) for PDEs <cit.> and deep learning <cit.>. This problem can be cast into the general framework. Suppose that D ⊆^d. We then let C = 1, = H^1_ϱ(D) be the Sobolev space of order one, _0 = C(D), _1 = ^d+1 and L_1 = L be defined by L(θ)(f) = (f(θ) , ∇ f(θ)^⊤ )^⊤.

The main difference with (i) is that the samples are vector-valued. Further generalizations are also possible. For example, in some cases it may be too expensive to evaluate the gradient at every point. Let m_1 be the number of function samples and m_2 be the number of function-plus-gradient samples. As shown in <cit.>, we can consider this as a multimodal sampling problem with C = 2, _1 =, _2 = ^d+1 and sampling operators L_1(θ)(f) = f(θ) and L_2(θ)(f) = (f(θ) , ∇ f(θ)^⊤ )^⊤.


(iii) Hilbert-valued function approximation from pointwise samples. In some applications, the unknown is a function taking values in a Hilbert space . This arises in, for instance, parametric PDEs and UQ, where f is the parametric solution map of a PDE whose (weak) solutions take values in . Approximating such a function from pointwise samples is easily considered in this framework. We use the setup of the scalar-valued case described above, except with _1 = replaced by _1 =. One can also consider gradient-augmented samples, much as in the previous example. A further extension to this problem is that of operator learning <cit.>, in which the unknown is an operator between two Hilbert spaces. This also fits within this framework.




(iv) Image reconstruction. We now briefly describe a seemingly quite different problem, based on <cit.>. This example highlights that the general framework can handle both discrete and continuous settings, and measurements that do not arise as pointwise samples.
Consider a discrete d-dimensional image of size n ×⋯× n, which we may vectorize and express as a vector f ∈^N, where N = n^d. Let F ∈^N × N be the matrix of the d-dimensional discrete Fourier transform. In Fourier imaging <cit.>, the goal is to recover f from a subset of its frequencies. If Ω⊆{1,…,N}, | Ω | = m, is the set of frequencies sampled, then the measurements of f are

P_Ω F f + e ∈^m,

where e ∈^m is noise and P_Ω∈^m × N is a matrix that selects the rows of F corresponding to the indices in Ω.
Fourier imaging arises in various applications, including Magnetic Resonance Imaging (MRI), Nuclear Magnetic Resonance (NMR) and radio interferometry <cit.>. A key question is how to best choose Ω within any relevant physical constraints.
As described in <cit.>, this problem can be cast into the general framework. The framework can also handle various practical constraints – for instance, the fact that MR imaging devices cannot sample individual frequencies, but may only sample along piecewise smooth curves in frequency space, which can be handled by considering vector-valued measurements. Moreover, the  framework can also handle the more advanced scenario of parallel MRI, where C ≥ 1 coils simultaneously acquire measurements.

(v) Other examples. Various other families of problems can be considered within this framework. For instance, as discussed in <cit.>, many standard measurement constructions in compressed sensing <cit.> become special cases of this approach. One can also readily consider related problems such as matrix completion and, more generally, matrix recovery for linear measurements <cit.>. Various other problems in sampling theory and signal processing, such as so-called mobile sampling <cit.>, also fit into this framework. It can also incorporate recovery problems involving function averages <cit.>, as well as techniques such as stratification and antithetics, which are commonly variance reduction techniques in Monte Carlo integration <cit.>.


(vi) Nonlinear approximation spaces. Many recovery problems call for nonlinear approximation spaces. A standard example is the sparse regression problem. Here, one may consider the setup of <ref> (or a more general setup, as in (ii)–(v)), but, rather than a linear subspace, one defines

= { ∑_i ∈S c_i ϕ_i : c_i ∈, ∀i, S ⊆{1,…,N}, | S | = n },

where N ≥ n ≥ 1 and {ϕ_i }^N_i=1⊂ L^2_ϱ(D) is some known dictionary of functions. The sparse regression problem has been studied extensively (see <cit.> and references therein), especially in the context of dictionaries of polynomials, where it is termed sparse polynomial approximation <cit.>. However, many other nonlinear approximations are widely used in practice. A partial list includes various `structured' sparse models, such as joint, group or block sparsity or sparsity in levels <cit.>, low-rank matrix or tensor models <cit.>, single- <cit.> or multi-layer neural networks, tensor networks <cit.>, Fourier sparse functions <cit.> and spaces defined by generative models <cit.>.










 §.§ The generalized Christoffel function


The key tool in our analysis is a generalization of the Christoffel function (Definition <ref>).


[Generalized Christoffel function]
Let _0 be a normed vector space,  be a Hilbert space, (D,,ϱ) be a measure space, L : D →(_0,) be such that the function θ∈ D ↦ L(θ)(f) ∈ is measurable for every f ∈_0 and suppose that ⊆_0, ≠{ 0 }. The generalized Christoffel function of  with respect to L is the function = (P,L) : D → defined by

(θ) = (,L)(θ) = sup{ L(θ)(f)^2_/f^2_ : f ∈, f ≠0 },  ∀θ∈D.


Notice that  reduces to the standard Christoffel function Kappa-def in the case of (i) above. In general,  measures how large the measurement L(θ)(f) of an arbitrary f ∈ can be (in norm) at an index θ∈ D in relation to the norm of f. For instance, in the Fourier imaging problem (iv) it measures how large the Fourier transform can be at a given frequency for an element of the approximation space in relation to its norm. We remark in passing that  inherits some of the properties of the standard Christoffel function. See, e.g., <cit.>.

Much like in kappa-w-def, given a nonnegative weight function w : D →, we also define


κ_w = κ_w(, L) = _θ∼ϱ w(θ) (,L)(θ)






 §.§ Theoretical guarantee


We now present a theoretical guarantee for this framework. We first require several assumptions.


[Nondegeneracy of the sampling operators]

For each c = 1,…,C and f ∈_0, the map θ∈ D_c ↦ L_c(θ)(f) ∈_c is measurable and θ∈ D_c ↦L_c(θ)(f)^2__c∈ is integrable. Moreover, there are constants 0 < a ≤ b <∞ such that


a f^2_ ≤∑^C_c=1 ∫_D_c L_c(θ_c)(f)^2__c ϱ_c(θ_c) ≤b f^2_,  ∀f ∈_0.


We remark in passing that the lower bound nondegeneracy is, in fact, only required to hold for f ∈', where ' =  -  = { p_1 - p_2 : p_1,p_2 ∈} is the difference set (see the proof of Theorem <ref> below, as well as that of <cit.>). 
This assumption says that the action of the sampling operators preserves the norm of any f ∈_0, up to constants. Note that it holds trivially with a = b = 1 in the standard problem (i), since in that case we have C = 1 and

∫_D_1 L_1(θ_1)(f)^2__1 ϱ_1(θ_1) = ∫_D |f(x)|^2 ϱ(x) = f^2_L^2_ϱ(D) = f^2_.

All other examples discussed in <ref> can also be formulated so that this assumption also holds. 

Recall that in the standard case studied earlier, stability and accuracy is ensured by the sampling discretization samp-disc. The middle term in this inequality is an empirical approximation to the L^2_ϱ-norm. An analogous construct is crucial to the analysis of this general setting, and arises by making an empirical approximation to the integrals in nondegeneracy. Given {θ_ic : i = 1,…,m_c, c = 1,…,C}, we say that empirical nondegeneracy holds for  with constants 0 < α_w ≤β_w < ∞ if


α_w q_ ≤√(∑^C_c=1 1/m_c ∑^m_c_c=1 w(θ_ic) L_c(θ_ic)(q)^2__c ) ≤β_w q_,  ∀q ∈' = - .

This can be seen as a generalization of the well-known Restricted Isometry Property (RIP) in compressed sensing <cit.>. In the case of sparse regression (see (vi) above), it is sometimes termed a universal sampling discretization <cit.>.

[Union of subspaces model]

The set ' =  - satisfies the following.


  (i) ' is a cone, i.e., t p ∈' for any t ≥ 0 and p ∈').

  (ii) ' ⊆_1 ∪⋯∪_d = :, where each _i ⊆_0 is a subspace of dimension n.




This trivially holds with d = 1 and _1 = ' = when  is an n-dimensional subspace. In general, Assumption <ref> is a extension of the union-of-subspaces model, which is well-known in the context of compressed sensing <cit.>. It includes many of the examples of nonlinear model classes used in practice, such as sparse regression problem and many of its generalizations (see (vi) above).

Finally, for the result below we also need a truncation operator. Similar to <ref>, this is the operator _σ : → defined by _σ(g) = min{ 1 , σ / g_} g, where σ≥ 0. Given a minimizer f̂ of general-least-squares, we define f̂^𝗍𝖾 = _σ(f̂). Note that we do not consider the other estimator f̂^𝖼𝖾 introduced in <ref>, although it could be readily formulated in this setting. The reason is that when  is a nonlinear space, it is generally not possible to certify (in polynomial time) that emp-nondegen holds, which, like samp-disc is the key ingredient for accuracy and stable recovery. In sparse recovery, for instance, this is equivalent to certifying that a given matrix has the RIP – a well-known NP-hard problem.




Consider the setup of <ref> and suppose that Assumptions <ref> and <ref> hold. Let 0 < ϵ < 1, κ_w_c be as in kappa-w-def-general and suppose that


m_c ≳a^-1 ·κ_w_c('-' , L_c) ·( log(2d/ϵ) + n ),   c = 1,…,C,

or, if  is as in Assumption <ref>(ii), 


m_c ≳a^-1 ·κ_w_c(, L_c) ·log(2nd/ϵ),    c = 1,…,C.

Then, for any f ∈_0, σ≥f_ and noise { e_ic}, the estimator f̂^𝗍𝗋 satisfies 

f - f̂^𝗍𝗋^2_ ≲(b/a) ·inf_p ∈ f - p_ + e^2_2 / a  + σ^2 ϵ,

where e^2_2 = ∑^C_c=11/m_c∑^m_c_i=1e_ic^2__c.



Note that mc-cond-1 involves the term κ_c evaluated over ' - '. One can replace this with just ', at the cost of a more complicated log term (see the proof of Theorem <ref> and <cit.>).




 §.§ Christoffel sampling


Much as in <ref>, we can use Theorem <ref> to optimize the sampling measures μ_c. Let = ' - ' in the case of mc-cond-1 or = in the case of mc-cond-2. Then, using kappa-w-def-general, we now choose

w_c(θ) = ( 1/2 + 1/2 (, L_c)(θ) /∫_D_c (, L_c)(θ) ϱ_c(θ)  )^-1,

which gives the sampling measures

μ_c(θ) = ( 1/2 + 1/2 (, L_c)(θ) /∫_D_c (, L_c)(θ) ϱ_c(θ)  ) ϱ_c(θ),  c = 1,…,C.

We term this Christoffel sampling.[Note we may assume without loss of generality that ∫_D_c( , L_c)(θ) ϱ_c(θ)  > 0. If not, the sampling operator L_c simply yields zero measurements over the space  almost everywhere, and can therefore be excluded. Nondegeneracy nondegeneracy implies that there is at least one sampling operator yielding nonzero measurements over .]
Substituting this into mc-cond-1 yields the measurement condition


m_c ≳( ∫_D_c (' - ' , L_c)(θ) ϱ_c(θ) ) ·(log(2d/ϵ) + n) 

or, in the case of mc-cond-2,


m_c ≳( ∫_D_c (, L_c)(θ) ϱ_c(θ) ) ·log(2nd/ϵ) .

This approach is `optimal' in the sense that it minimizes (up to a factor of 2) the bound mc-cond-1 over all possible sampling measures μ_c. When  is an n-dimensional subspace – in which case ' - ' = ' = and  in Assumption <ref> can be chosen as = – it is a short argument using the nondegeneracy condition nondegeneracy to see that

∑^C_c=1 ∫_D_c (, L_c)(θ) ϱ_c(θ) ≤b n.

(see <cit.>). Hence, if each m_c is chosen proportional to the right-hand side in mc-cond-2-opt, then the total number of measurements satisfies the near-optimal log-linear scaling

m = m_1 + ⋯+ m_C ≲(b/a) ·n ·log(2n/ϵ).

Unfortunately, in the general case of a nonlinear set , there is no clear way to relate the integral in mc-cond-1-opt to explicit quantities such as n and d. However, it is possible to show the bound

m = m_1 + ⋯+ m_C ≲(b/a) ·n ·d ·log(2n/ϵ).

(see <cit.> once more),
where we recall that d is the number of subspaces in Assumption <ref>(ii). For fixed and small d, this is near-optimal. However, in cases such as sparse regression, d ≫  1. Fortunately, a more refined analysis is possible in these cases. See <cit.> for discussion.

While it is difficult to provide explicit measurement conditions in the general case, it is possible to gain some insight over why Christoffel sampling, in general, improves over Monte Carlo sampling, i.e., the case where μ_c = ϱ_c, ∀ c. Since w_c ≡ 1 in this case, kappa-w-def-general and  Theorem <ref> provides the measurement conditions for Monte Carlo sampling of the form

m_c ≳( _θ∼ϱ_c (' - ' , L_c)(θ) )·(log(2d/ϵ) + n ),  c = 1,…,C,

and likewise in the case of mc-cond-2.
Therefore, comparing with mc-cond-1-opt, the improvement of Christoffel sampling, in general terms, can be equated to the difference between the supremum of the (generalized) Christoffel function and its integral (mean). In particular, if this function is sharply peaked, then we expect significant improvements, while if it is approximately flat, then we expect less improvement. Such observations are witnessed in numerical experiments <cit.> .


















§ CONCLUSIONS AND OUTLOOK


In this article, we have surveyed recent advances in optimal sampling for least-squares approximation, focusing on the key role that the Christoffel function plays in both understanding the sample complexity in general, and in constructing near-optimal sampling strategies. We have also seen in <ref> how these ideas naturally extend to more general types of measurements and play a key role in the sample complexity even for nonlinear spaces. We now offer some concluding thoughts.

First, although the picture for pointwise sampling in linear spaces is increasingly mature, there remain various open questions. While optimal (i.e., n) sampling strategies that are practical (i.e., implementable in polynomial time) are now known (recall <ref>), future investigations are needed on their practical efficacy, especially in the interpolation regime. There are also open questions about uniform recovery, which are especially relevant to sampling numbers and questions in information-based complexity. Finally, the question of optimal sampling with hierarchical or adaptive schemes has not yet been addressed.

By contrast, nonlinear approximation spaces pose many more open problems. First, even in relatively well-studied settings such as sparse regression, it is unknown whether Christoffel sampling generally leads to near-optimal sample complexity bounds. This issue is discussed in detail in <cit.>. Less is known about more complicated nonlinear spaces. See <cit.> for discussion. Second, there is also the practical matter of drawing samples from the resulting Christoffel sampling measure. This is very dependent on the particular nonlinear space and samples under consideration, and may well be highly nontrivial. Finding practical surrogate sampling measures is an interesting open problem. 

In the nonlinear setting, it is worth noting that Christoffel sampling is well-suited only when the approximation space  has low intrinsic complexity that is comparable to the number of samples (e.g., m ≍ n log(n) in the case where  is a linear space with () = n). It is not well suited for `complex' approximation spaces, such as spaces of deep neural networks <cit.> or low-rank tensor networks <cit.>. Christoffel sampling can be implemented in an adaptive manner in such cases. Here, one alternates between adding samples and learning an approximation, and at each stages uses a linearization to obtain an intermediate linear approximation space over which Christoffel sampling can be performed. This idea was developed for approximating functions and solving PDEs via deep learning in <cit.> and <cit.>, and later extended to more general approximation spaces in <cit.>. 

Finally, it is worth noting that Christoffel sampling, in any of its guises, is not a panacea. Depending on the problem or target function, there may be little or no benefit over standard Monte Carlo sampling. This is relevant in various applications, as Monte Carlo samples are often commonly encountered in practice <cit.>. It leads to another interesting line of research, which is understanding function classes where Monte Carlo sampling is near-optimal. One such case is
holomorphic function approximation in high dimensions. In a series of works <cit.>, it has been shown that Monte Carlo sampling is near-optimal information for classes of infinite-dimensional holomorphic functions arising in parametric DE problems. In particular, Monte Carlo sampling yields error bounds that are within a polylogarithmic factor of the adaptive m-width for such classes. For other work in this direction involving functions Sobolev spaces, see <cit.>.




§ ACKNOWLEDGEMENTS


The idea for this article came from a talk given at the 2023 Foundations of Computational Mathematics conference. The author would like to thank the conference organizers, as well as the various participants for insightful comments and questions. Parts of this article were written during a visit to the Isaac Newton Institute for Mathematical Sciences, Cambridge for the programme Discretization and recovery in high-dimensional spaces. The author would like to thank the institute for support and hospitality, and the programme participants for providing a stimulating environment. Finally, he would like to thank Simone Brugiapaglia, Matthew J. Colbrook, Maksym Neyra–Nesterenko and Daniel Fassler for helpful feedback.

This work was supported by EPSRC grant number EP/R014604/1 and NSERC grant number RGPIN-2021-611675. 



plain








§ PROOFS



In this appendix we give proofs for various results in the paper. We commence with Lemma <ref>. This is a standard result. We include a short proof for completeness.


[Proof of Lemma <ref>]

Recall that α_w = σ_min(A), where A is the matrix defined in (<ref>). This matrix is full rank since m ≥ n and α_w > 0, and therefore the LS problem has a unique solution. 

Now let p ∈ be arbitrary and consider the variational form variational-form applied to the element f̂ - p ∈. Subtracting pf̂ - p_𝖽𝗂𝗌𝖼,w from both sides gives

f̂ - p^2_𝖽𝗂𝗌𝖼,w = f - pf̂ - p_𝖽𝗂𝗌𝖼,w +  1/m ∑^m_i=1 w(x_i) e_i f̂(x_i) - p(x_i).

Several applications of the Cauchy-Schwarz inequality now yield


f̂ - p_𝖽𝗂𝗌𝖼,w ≤f - p_𝖽𝗂𝗌𝖼,w + e_2,w. 

We now use the triangle inequality and MZ-inequality to get

f - f̂_L^2_ϱ() ≤f - p_L^2_ϱ() + f̂ - p_L^2_ϱ() ≤f - p_L^2_ϱ() + 1/α_w f̂-p_𝖽𝗂𝗌𝖼,w.

The result now follows from disc-err-bd.


We next prove Theorem <ref>. This has now become a standard exercise involving the matrix Chernoff bound (see, e.g., <cit.>), which we repeat here for convenience. This bound was first used in the context of least-squares approximation from i.i.d. samples in <cit.>.



[Matrix Chernoff bound]

Let X_1,…,X_m be independent, self-adjoint random matrices of dimension n.  Assume that 

X_i ≽0
  and  λ_max(X_i) ≤R

almost surely for each i=1,…,m, and define

μ_min = λ_min ( ∑^m_i=1 X_i )   and  μ_max = λ_max ( ∑^m_i=1 X_i ).

Then, for 0 ≤δ≤ 1,

( λ_min  ( ∑^m_i=1 X_i ) ≤( 1 - δ) μ_min ) ≤n exp( -μ_min ( (1-δ) log(1-δ) + δ)/R ),

and, for δ≥ 0,

( λ_max  ( ∑^m_i=1 X_i ) ≥( 1 + δ) μ_max ) ≤n exp( -μ_max ( (1+δ) log(1+δ)-δ)/R ).





[Proof of Theorem <ref>]

Let {ϕ_i }^n_i=1 be an orthonormal basis for  and A be as in ls-Ab. Observe that

A^* A = ∑^m_i=1 X_i,    X_i : =  1/m ( w(x_i) ϕ_j(x_i) ϕ_k(x_i) )^n_j,k=1 ,

is a sum of independent random matrices. Also, orthonormality of the basis functions implies that alpha-beta-sigma holds, i.e.,

α^2_w = λ_min(A^*A),    β^2_w = λ_max(A^*A).

We now wish to apply Theorem <ref>. Due to mu_weight_fn, the fact that w = 1/ν and orthonormality, we have

((A^*A))_jk = ∑^m_i=1 ((X_i))_jk = 1/m ∑^m_i=1 ∫_D w(x) ϕ_j(x) ϕ_k(x) μ_i(x) = ∫_D ϕ_j(x) ϕ_k(x) w(x) ν(x) ϱ(x) = δ_jk.

Hence (A^*A) = I and therefore μ_min = μ_max = 1. Next, let c ∈^n be arbitrary. Then

c^* X_i c = w(x_i)/m| ∑^n_j=1 c_j ϕ_j(x_i) |^2 = w(x_i)/m | p(x_i) |^2,    where p = ∑^n_j=1 c_j ϕ_j.

We immediately deduce that X_i≽ 0. Moreover, Parseval's identity implies that c^2_2 = p^2_L^2_ϱ(D). Hence, using this and Kappa-def and kappa-w-def, we obtain 

c^* X_i c ≤w(x_i)/m K()(x_i) p^2_L^2_ϱ(D) ≤κ_w()/m c^2_2.

Since c was arbitrary and X_i ≽ 0, we conclude that

λ_max(X_i) ≤R : = κ_w()/m.

We are now ready to apply Theorem <ref>. Using the union bound, we have

( α_w ≤√(1-δ) or β_w ≥√(1+δ) )     ≤(λ_min(A^*A) ≤1-δ) + ( λ_max(A^*A) ≥1-δ)


    ≤n exp( - m a_δ/ κ_w() ) + n exp( - m b_δ/ κ_w() ),

where a_δ =  ( (1-δ) log(1-δ) + δ ) and b_δ = ( (1+δ) log(1+δ)-δ). Notice that a_δ≥ b_δ and b_δ = 1/c_δ, where c_δ is as in m-bound-alpha-beta. Therefore

( α_w ≤√(1-δ) or β_w ≥√(1+δ) ) ≤2 n exp( -m/c_δ κ_w() ) ≤ϵ,

where in the last step we used m-bound-alpha-beta. This completes the proof.


We next prove Corollary <ref>. For this, we require the following lemma, which is based on <cit.> (which is, in turn, based on <cit.>).



Let ⊂ L^2_ϱ(D) be an n-dimensional subspace with orthonormal basis {ϕ_i }^n_i=1 and μ_1,…,μ_m be probability measures satisfying Assumption <ref>. Consider sample points drawn randomly and independently with x_i ∼μ_i, i = 1,…,m. Then

∑^n_i=1 | gϕ_i_𝖽𝗂𝗌𝖼,w |^2 ≤κ_w()/m g^2_L^2_ϱ(D),  ∀g ∈^⊥,

where  w = 1/ν, ··_𝖽𝗂𝗌𝖼,w and κ_w are as in mu_weight_fn,  semi-inner-product and kappa-w-def, respectively.


Let g ∈^⊥ and ł∈{1,… , n } be arbitrary. Then

|gϕ_l_𝖽𝗂𝗌𝖼,w|^2  =      1/m^2 ∑^m_i,j = 1 [ w(x_i) w(x_j) g(x_i) g(x_j) ϕ_l(x_i) ϕ_l(x_j) ]


 =      1/m^2  ∑^m_i,j = 1 
 i ≠j ( w(x_i) g(x_i) ϕ_l(x_i) ) ( w(x_j) g(x_j) ϕ_l(x_j) ) 
 + 1/m^2 ∑^m_i=1 ( | w(x_i) g(x_i) ϕ_l(x_i) |^2 ) .

Now mu_weight_fn implies that

1/m ∑^m_i=1 (w(x_i) g(x_i) ϕ_l(x_i) ) = 1/m ∑^m_i=1 ∫_D w(x) g(x) ϕ_l(x) μ_i(x) = gϕ_l_L^2_ϱ(D) = 0,

since g ∈^⊥. Therefore

|gϕ_l_𝖽𝗂𝗌𝖼,w|^2    =  1/m^2 ∑^m_i=1 [ ( | w(x_i) g(x_i) ϕ_l(x_i) |^2 ) - ( (w(x_i) g(x_i) ϕ_l(x_i) ) )^2 ]


    ≤1/m^2 ∑^m_i=1 ∫_D (w(x))^2 | g(x) |^2 | ϕ_l(x) |^2 μ_i(x)
= 1/m ∫_D w(x) |g(x)|^2 | ϕ_l(x) |^2 ϱ(x),

where in the final step we used mu_weight_fn once more. 
We now apply Kappa-def-alt and kappa-w-def to obtain

∑^n_l=1 |gϕ_l_𝖽𝗂𝗌𝖼,w|^2 ≤1/m ∫_D w(x) |g(x)|^2 ∑^n_l=1 | ϕ_l(x) |^2 ϱ(x)  ≤κ_w()/m g^2_L^2_ϱ(D),

as required.



[Proof of Corollary <ref>]

The condition on m and Theorem <ref> imply that alpha-beta-delta holds with probability at least 1-ϵ/2. Therefore Lemma <ref> asserts that the weighted least-squares problem has a unique solution for any function that is defined at the sample points x_i and any noise vector. Let f ∈ L^2_ϱ(D). Then f is defined at the x_i with probability one. Now consider arbitrary noise e = (e_i)^m_i=1∈^m and write f̂_e ∈ for the corresponding weighted least-squares approximation from noisy samples f(x_i) + e_i. Let p^* ∈ be the (unique) element such that f - p^*_L^2_ϱ(D) = inf_p ∈f - p_L^2_ϱ(D) and write g = f - p^*. Let ĝ = ĝ_0 be the least-squares approximation to g from noiseless samples g(x_i) and 0̂_e ∈ be the least-squares approximation to the zero function from the noisy samples e_i. Notice that p^* = p^*_0 = p^*, since the weighted least-squares approximation is a projection in the discrete inner product. Since the least-squares approximation is also linear, we have f - f̂_e = f - f̂ - 0̂_e = g - ĝ - 0̂_e.
This gives


f-f̂_e_L^2_ϱ(D) ≤g_L^2_ϱ(D) + ĝ_L^2_ϱ(D) + 0̂_e_L^2_ϱ(D).

Lemma <ref> and the fact that alpha-beta-delta holds imply that


0̂_e_L^2_ϱ(D) ≤1/√(1-δ) e_2,w.

Now consider the term ĝ_L^2_ϱ(D). Writing ĝ = ∑^n_i=1ĉ_i ϕ_i and ĉ = (ĉ_i)^n_i=1 and using the variational form variational-form gives

ĝ^2_𝖽𝗂𝗌𝖼,w = gĝ_𝖽𝗂𝗌𝖼,w= ∑^n_j=1 c_j gϕ_j_𝖽𝗂𝗌𝖼,𝗐 ≤c_2 √(∑^n_j=1 |gϕ_j_𝖽𝗂𝗌𝖼,𝗐|^2 ).

Now Parseval's identity and the fact that α_w ≥√(1-δ) gives

ĝ^2_L^2_ϱ(D) ≤1/1-δ ĝ^2_𝖽𝗂𝗌𝖼,w ≤ĝ_L^2_ϱ(D) √(∑^n_j=1 |gϕ_j_𝖽𝗂𝗌𝖼,𝗐|^2 ),

i.e.,


ĝ_L^2_ϱ(D) ≤1/1-δ √(∑^n_j=1 |gϕ_j_𝖽𝗂𝗌𝖼,𝗐|^2 ).

We now bound this term in probability. Consider the random variable X = ∑^n_j=1 |gϕ_j_𝖽𝗂𝗌𝖼,𝗐|^2. Since p^* is the orthogonal projection of f onto , g = f - p^* ∈^⊥ and Lemma <ref> implies that

(X) ≤κ_w()/m g^2_L^2_ϱ(D).

Hence, by Markov's inequality

( ĝ_L^2_ϱ(D) ≥1/1-δ √(2 κ_w()/m ϵ) g_L^2_ϱ(D) ) ≤( X ≥2 (X)/ϵ  ) ≤ϵ/2.

We deduce that 

ĝ_L^2_ϱ(D) ≤1/1-δ  √(2 κ_w()/m ϵ) g_L^2_ϱ(D),

with probability at least 1-ϵ/2. Substituting this and q-bd into main-err-split we deduce, after an application of the union bound, that

f-f̂_e_L^2_ϱ(D) ≤( 1 +  1/1-δ  √(2 κ_w()/m ϵ) g_L^2_ϱ(D) ) g_L^2_ϱ(D) + 1/√(1-δ) e_2,w,

with probability at least 1-ϵ. This completes the proof.



We next prove Corollary <ref>. This follows <cit.> and employs Bernstein's inequality.



[Proof of Corollary <ref>]
Let E be the event that alpha-beta-delta holds Let p = p^* be a polynomial attaining the infimum in ls-err-bd-prob-2. Now let F be the event that

f - p_𝖽𝗂𝗌𝖼,w ≤√(2) ( f - p_L^2_ϱ(D) + √(w)(f-p)_L^∞_ϱ(D)/√(k)  ).

Suppose that E and F occur. Then the definitions of E, F and Lemma <ref> imply that

f - f̂_L^2_ϱ(D)     ≤f - p_L^2_ϱ(D) + 1/√(1-δ) f - p_𝖽𝗂𝗌𝖼,w + 1/√(1-δ) e_2,w


    ≤( 1 + √(2/1-δ) ) f - p_L^2_ϱ(D) + √(2/1-δ)  f-p_L^∞_ϱ(D)/√(k) + 1/√(1-δ) e_2,w.

This yields the desired bound ls-err-bd-prob-2. Hence, by the union bound, it suffices to show that (E^c),(F^c) ≤ϵ /2.

The fact that (E^c) ≤ϵ / 2 follows immediately from the first condition in on m in m-conds-in-prob-2 and Theorem <ref>. We now consider (F^c). Define the random variables

Z_i = w(x_i)| f(x_i) - p(x_i) |^2  and  X_i = Z_i - (Z_i).

Notice that

1/m ∑^m_i=1 (Z_i) = f-p^2_L^2_ϱ(D) : = a

due to exp-sum-scaling, and therefore

f -p^2_𝖽𝗂𝗌𝖼,w = 1/m ∑^m_i=1 Z_i = 1/m ∑^m_i=1 X_i + a.

The idea now is to use Bernstein's inequality to estimate the random variable ∑^m_i=1 X_i. Let

b = _x ∼ϱ w(x) | f(x) - p(x) |^2 ≡√(w)(f - p)^2_L^∞_ϱ(D)

and notice that X_i ≤ Z_i ≤ b and -X_i ≤(Z_i) ≤ b
almost surely. Hence |X_i| ≤ b almost surely. We also have 0 ≤ Z_i ≤ b almost surely, and therefore

∑^m_i=1 ((X_i)^2) ≤∑^m_i=1 ((Z_i)^2) ≤b ∑^m_i=1 (Z_i) = a b m.

Since (X_i) = 0, we may apply Bernstein's inequality (see, e.g., <cit.>) to get

( | 1/m ∑^m_i=1 X_i | ≥t ) ≤2 exp( -t^2 m/2/a b + b t /3 ),  ∀t > 0.

Set t = a + b/k. Then it is a short argument involving the second condition in m-conds-in-prob-2 to show that t^2 m/2/a b + b t /3≥log(4/ϵ). Therefore,

f - p^2_𝖽𝗂𝗌𝖼,w ≤| 1/m ∑^m_i=1 X_i | + a <  2 a + b/k,

with probability at least 1-ϵ/2. Substituting the values for a, b and using the inequality √(s+t)≤√(s) + √(t), we see that 

f - p_𝖽𝗂𝗌𝖼,w ≤√(2) ( f - p_L^2_ϱ(D) + √(w)(f-p)_L^∞_ϱ(D)/√(k) ),

with probability at least 1-ϵ/2. Therefore, (F^c) ≤ϵ/2, as required.



We now prove Lemma <ref> and Theorem <ref>. These ideas go back to <cit.>, but the specific arguments are based on <cit.>.




[Proof of Lemma <ref>]

The setup is the same as the proof of Corollary <ref>. However, since we now square the error terms, we use the fact that g ∈^⊥ to replace main-err-split by

f-f̂_e^2_L^2_ϱ(D) = g^2_L^2_ϱ(D) + ĝ + 0̂_e^2_L^2_ϱ(D) ≤g^2_L^2_ϱ(D) + 2 ĝ^2_L^2_ϱ(D) + 20̂_e^2_L^2_ϱ(D).

Whenever G - I_2≤δ holds we have that alpha-beta-delta holds, and therefore q-bd and g-bd-1 also hold. This implies that 

( f-f̂^2_L^2_ϱ(D) χ_G-I_2 ≤δ ) ≤g^2_L^2_ϱ(D) +  2/1-δ e^2_2,w + 2/(1-δ)^2 ∑^n_l=1 |gϕ_l_𝖽𝗂𝗌𝖼,w|^2.

The result now follows from Lemma <ref>.




[Proof of Theorem <ref>]
Let E be the event that G - I_2 ≤δ. Suppose that E occurs. Then

f - f̂^𝖼𝖾_L^2_ϱ(D) = f - f̂_L^2_ϱ(D),    f - f̂^𝗍𝖾_L^2_ϱ(D) ≤f - f̂_L^2_ϱ(D).

Here, the second bound follows from the facts that f = _σ(f) and _σ is a contraction in the L^2_ϱ(D) norm. On the other hand, if E does not occur then we have

f - f̂^𝖼𝖾_L^2_ϱ(D) = f_L^2_ϱ(D),    f - f̂^𝗍𝖾_L^2_ϱ(D) ≤f +  f̂^𝗍𝖾_L^2_ϱ(D) ≤2 σ.

Now observe that (E^c) ≤ϵ, due to the assumption on m and Theorem <ref>. The result now follows by the law of total expectation and Lemma <ref>.





[Proof of Theorem <ref>]
The weight function w corresponding to opt-meas-krieg is given by


w(x) = ( 1/2 + 1/4 ∑^n_i=1 | ϕ_i (x) |^2/n + 1/4 ∑^∞_l=0 v^2_l/|I_l| ∑_i ∈I_l | ϕ_i(x) |^2 )^-1.

Therefore, by Kappa-def-alt, we have w(x) ≤ 4 / ()(x), which gives
 
κ_w = _y ∼ϱ w(x) ()(x)  ≤4 n.

We deduce from Theorem <ref> and m-est-krieg that √(1-δ) < α_w ≤β_w ≤√(1+δ) with probability at least 1-ϵ/2. Using this and the fact that w(x) ≤ 2 by construction, we deduce from Lemma <ref> that

f - f̂_L^2_ϱ(D) ≤inf_p ∈ { f - p_L^2_ϱ(D) + 1/√(1-δ) f - p_𝖽𝗂𝗌𝖼,w } + √(2/1-δ) e_2
 
for any f and e. Let c_i = fϕ_i_L^2_ϱ(D) be the coefficients of f and p = ∑^n_i=1 c_i ϕ_i be the best approximation to f from . Then this gives

 
f - f̂_L^2_ϱ(D) ≤e_n(f) + 1/√(1-δ) f - p_𝖽𝗂𝗌𝖼,w + √(2/1-δ) e_2.

For l = 0,1,2,…, define the matrices

A^(l) = ( √(w(x_i)/m) ϕ_j(x_i) )^m,2^l+1 n_i = 1,j=2^l n + 1 ∈^m ×2^l n,
    c^(l) = ( c_i )^2^l+1 n_i=2^l n + 1 ∈^2^l n.

Then

f - p_𝖽𝗂𝗌𝖼,𝗐 ≤∑^∞_l=0 A^(l) c^(l)_2 ≤∑^∞_l=0 A^(l) _2 c^(l)_2 ≤∑^∞_l=0 A^(l)_2 e_2^l n(f),

where e_2^l n(f) is as in e-def.

We now wish to use the matrix Chernoff bound to estimate A^(l)_2. Observe that 

( A^(l)^2_2 ≥1+t ) = ( λ_max( 
(A^(l))^* A^(l) ) ≥1+t ).

Now, as in the proof of Theorem <ref>, note that ( (A^(l))^* A^(l)  ) = I
and

(A^(l))^* A^(l) = ∑^m_i=1 X_i,  where X_i = ( w(x_i)/m ϕ_j(x_i) ϕ_k(x_i) ) _j,k ∈I_l.

The matrices X_i are nonnegative definite and satisfy, for any c = (c_j)_j ∈ I_l,

c^* X_i c = w(x_i)/m | ∑_j ∈I_l c_j ϕ_j(x_i) |^2 ≤w(x_i)/m ∑_j ∈I_l | ϕ_j(x_i) |^2 c^2_2 .

Using w-krieg-def and taking the supremum over all such c with c_2 = 1, we deduce that

λ_max(X_i) ≤4 |I_l|/m v^2_l.

Hence, the matrix Chernoff bound (Theorem <ref>) gives that

( A^(l)^2_2 ≥1+t ) ≤|I_l| exp( - m v^2_l h(t)/4 |I_l| ),

where h(t) = (1+t) log(1+t) - t.
Now let ϵ_l = (3/π^2) ϵ / (l+1)^2, so that ∑^∞_l=0ϵ_l = ϵ / 2. We want to choose t = t_l so that ( A^(l)^2_2 ≥ 1+t) ≤ϵ_l. Using the bound m-est-krieg, we see that

( A^(l)^2_2 ≥1+t_t ) ≤| I_l| exp( - c_δ n log(4n/ϵ) h(t_t) v^2_l/|I_l| ) ≤ϵ_l,

provided

h(t_l) ≥| I_l| log(|I_l| / ϵ_l) /c_δ n log(4n/ϵ) v^2_l .

Now h is increasing and h(t) ≥ h(1) t for t ≥ 1. Therefore it suffices to take

t_l =  | I_l| log(|I_l| / ϵ_l) /h(1) c_δ n log(4n/ϵ) v^2_l .

Now we recall that |I_l| = n 2^l, the definition of ϵ_l and c_δ≥ c_1 
≳ 1, to get, after some algebra,

1+t_l ≤c 2^l log(2^l (l+1)^2)/v^2_l

for some numerical constant c> 0. To summarize, we have shown that

( A^(l)^2_2 ≥c 2^l log(2^l (l+1)^2)/v^2_l ) ≤ϵ_l.

Taking the union bound and recalling that ∑^∞_l=0ϵ_l = ϵ/2, we deduce that

f - p_𝖽𝗂𝗌𝖼,𝗐 ≤∑^∞_l=0 A^(l)_2 c^(l)_2 ≤c ∑^∞_l=0 2^l/2 √(log(2^l(l+1)^2))/v_l  e_2^l n(f) ≤c ∑^∞_l=0 2^l/2 (l+1)^3/2/v_l e_2^l n(f),

with probability at least 1-ϵ/2 and a potentially different numerical constant c. Now consider e_2^l n(f). The terms e_k(f) are monotonically nonincreasing in k. Therefore, for l = 1,2,… we have

n (2^l-1) (e_2^l n(f))^p ≤e_n+1(f))^p + ⋯+ (e_2^l n(f))^p ≤∑_k > n (e_k(f))^p ,

Hence

e_2^l n(f) ≤c_p 2^-l/p ( 1/n ∑_k > n (e_k(f))^p )^1/p,  l = 1,2,….

We deduce that 

f - p_𝖽𝗂𝗌𝖼,𝗐 ≤c_p [ e_n(f)/v_0 + ( ∑_k ≥n (e_k(f))^p )^1/p  ( ∑^∞_l=1 2^l(1/2-1/p) l^3/2/v_l  ) ]
 ≤c_p,θ  ( ∑_k > n (e_k(f))^p )^1/p 

where in the final step we used the fact that v_l = 2^-θ l for 0 < θ < 1/p-1/2 to deduce that the final sum converges. Substituting this into f-krieg-bound-1 now gives the result.





[Proof of Theorem <ref>]

Much like in Lemma <ref>, if emp-nondegen holds with α_w > 0 then the estimator f̂ given by general-least-squares satisfies (see <cit.>)

f - f̂_ ≤inf_p ∈ { f-p_ + 2/α_w f-p_𝖽𝗂𝗌𝖼,w } + 2/α_w e_2,w,
 
where, for convenience, we define e^2_2,w = ∑^C_c=11/m_c∑^m_c_i=1 w_c(θ_ic) e_ic^2__c and

g^2_𝖽𝗂𝗌𝖼,w = ∑^C_c=1 1/m_c ∑^m_c_c=1 w(θ_ic) L_c(θ_ic)(g)^2__c ,  g ∈_0.

Now let E be the event that α_w ≥ a/2, where a is as in nondegeneracy (the choice of 1/2 here is arbitrary), and consider the estimator f̂^𝗍𝗋. We argue similarly to the proofs of Lemma <ref> and Theorem <ref>. Since _σ is a contraction, we have

f - f̂^𝗍𝗋_ ≤f - f̂_  and  f - f̂^𝗍𝗋_ ≤2 σ.

Now fix p ∈. Then

f - f̂^𝗍𝗋^2_     = ( f - f̂^𝗍𝗋^2_ | E ) (E) + ( f - f̂^𝗍𝗋^2_ | E^c ) (E^c)


    ≤2 f - p^2_ + 16/a^2 f - p^2_𝖽𝗂𝗌𝖼,w + 16/a^2 e^2_2,w + 4 σ^2 (E^c).

Now observe that

f - p^2_𝖽𝗂𝗌𝖼,w     = ( ∑^C_c=1 1/m_c ∑^m_c_i=1 w_c(θ_ic) L_c(f-p)(θ_ic)^2__c )


    =  ∑^C_c=1 ∫_D_c w_c(θ_c) L_c(f-p)(θ_c)^2__c μ_c(θ_c) 


    = ∑^C_c=1 ∫_D_cL_c(f-p)(θ_c)^2__c ϱ_c(θ_c)  ≤b^2 f-p^2_,

where in the last step we used nondegeneracy (Assumption <ref>). We also have

e^2_2,w = ∑^C_c=1 1/m_c ∑^m_c_i=1 e_ic^2__c ∫_D_c w_c(θ_c) μ(θ_c) = ∑^C_c=1 1/m_c ∑^m_c_i=1 e_ic^2__c,

Here we used the definition of the weight functions w_c and the fact that each ϱ_c is a probability measure. We deduce that

f - f̂^𝗍𝗋^2_ ≤2 f - p^2_ + 16 b^2/a^2 f - p^2_ + 16/a^2 e^2_2 + 4 σ^2 (E^c),

Hence, the result follows, provided (E^c) ≤ϵ.

To show that (E^c) ≤ϵ, we appeal to <cit.>, using conditions (b) and (c) defined therein. The remainder of the proof involves showing how to recast the setup considered in <ref> as a special case of that considered in <cit.>. Let m = m_1+⋯ + m_C. Now, for c = 1,…,C, let ^(c) be the distribution of operators in (_0,_c) defined by A^(c)∼^(c) if

A^(c)(f) = √(m/m_c) √(w_c(θ_c)) L_c(θ_c)(f),  where θ_c ∼μ_c.

Now, following <cit.>, define _1,…,_m by _i = ^(c) if m_1+⋯ + m_c-1 < i ≤ m_1 + ⋯ + m_c for c = 1,…,C. Doing this, the setup of <ref> is a special case of <cit.> with the family = {_i}^m_i=1. In particular, nondegeneracy in the sense of <cit.> is implied by nondegeneracy. 

We now apply <cit.>, and specifically, parts (b) and (c), with = ' and δ = 1/2. This implies that emp-nondegen holds, provided


m ≳a^-1/2 ·Φ(S('-') ; ) ·(log(2 d / ϵ) + n ),

or, with  as in Assumption <ref>(ii), 


m ≳a^-1/2 ·Φ(S() ; ) ·log(2 d n/ ϵ) . 

where Φ is the so-called variation, as defined in <cit.> and, for ⊆_0, S() = { u / u_ : u ∈, u ≠ 0 }. Consider any such set . Using this and the definition of , we see that

Φ(S() ; ) = max_c=1,…,C Φ(S() ; ^(c) ) =  max_c=1,…,C { m/m_c _θ_c ∼ϱ_c sup_u ∈\{ 0 } w_c(θ_c) L_c(θ_c)(u / u_)^2__c }.

Since L_c(θ_c) is linear, we deduce that

Φ(S() ; ) = max_c=1,…,C { m/m_c κ_w_c(; L_c) },

where κ_w_c is as in kappa-w-def-general. Using this and setting, respectively, = ' - ' or = we see that m-cond-gen-1 is equivalent to mc-cond-1 and m-cond-gen-2 is equivalent to mc-cond-2. The result now follows.










