<figure><img src="image_1.png" /></figure>

HAL
open science

# Randomized least-squares with minimal oversampling and interpolation in general spaces

Matthieu Dolbeault, Moulay Abdellah Chkifa

<figure><img src="image_2.png" /></figure>

## To cite this version:

Matthieu Dolbeault, Moulay Abdellah Chkifa. Randomized least-squares with minimal oversampling and interpolation in general spaces. 2023. hal-04243062v2

HAL Id: hal-04243062

https://hal.science/hal-04243062v2

Preprint submitted on 13 Feb 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Randomized least-squares with minimal oversampling and interpolation in general spaces

Matthieu Dolbeault*

Moulay Abdellah Chkifa†

## Abstract

In approximation of functions based on point values, least-squares methods provide more stability than interpolation, at the expense of increasing the sampling budget. We show that near-optimal approximation error can nevertheless be achieved, in an expected $L^2$ sense, as soon as the sample size $m$ is larger than the dimension $n$ of the approximation space by a constant ratio. On the other hand, for $m = n$, we obtain an interpolation strategy with a stability factor of order $n$. The proposed sampling algorithms are greedy procedures based on [BSS09] and [LS18], with polynomial computational complexity.

**Keywords.** Least squares, Interpolation, Christoffel function

**MSC 2020.** 65D15, 41A65, 65D05, 41A81, 41A05

## 1 Introduction and main results

Let $(X, A, \mu)$ be a probability space. We consider the problem of estimating an unknown function $f: X \to \mathbb{C}$ from evaluations of $f$ at chosen points $x_1, \dots, x_m \in X$. We assess the error between $f$ and its estimator $\tilde{f}$ either in the $L^2(X, \mu)$ norm

$$ \|g\|_{L^2} := \left( \int_X |g(x)|^2 d\mu(x) \right)^{1/2}, $$

or in the uniform norm $\|g\|_{L^\infty(X, \mu)}$. Given a subspace $V_n \subset L^2(X, \mu)$ of dimension $n$, we would like the estimator $\tilde{f}$ to belong to $V_n$ and to perform almost as well as the best approximation of $f$ in $V_n$, that is, its orthogonal projection

$$ P_n f = \arg \min_{v \in V_n} \|f - v\|_{L^2}. $$

As we only have access to point-wise observations, we cannot explicitly compute $P_n f$ in general. In this context, a classical approach consists in considering a solution to the weighted least-squares problem

$$ P_n^m f \in \arg \min_{v \in V_n} \sum_{i=1}^{m} s_i |f(x_i) - v(x_i)|^2, $$

<footnote>*RWTH Aachen, Germany (dolbeault@igpm.rwth-aachen.de). The author acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project number 442047500 through the Collaborative Research Center Sparsity and Singular Structures (SFB 1481).</footnote>

<footnote>†Mohammed VI Polytechnic University, Ben Guerir (abdellah.chkifa@um6p.ma).</footnote>
where we may use some weights $s_1, \dots, s_m > 0$. We would like this problem to admit a unique solution, and therefore require $m \ge n$. Similar to $P_n$ for the $L^2$ norm, the operator $P_n^m$ is the orthogonal projector onto $V_n$ with respect to the empirical norm

$$ \|g\|_m := \left( \sum_{i=1}^{m} s_i |g(x_i)|^2 \right)^{1/2} . $$

The approximation accuracy is inherently related to the points $x_i$ and the weights $s_i$. For the sake of illustration, consider the setting where $V_n = \mathbb{P}_{n-1}$ is the space of algebraic polynomials of degree less than $n$, restricted to the interval $X = [-1, 1]$, and choose $m = n$, so that $P_n^n$ is the Lagrange interpolation operator associated with $\{x_1, \dots, x_n\}$. For equally spaced points $x_i$, this corresponds to interpolation on a uniform grid, which is known to be highly unstable, failing to converge towards $f$, even when $f$ is infinitely smooth. This is the so-called Runge phenomenon, see e.g. [MM08]. The phenomenon persists with equally spaced points even for $m = rn$ with $r > 1$ constant, as observed in [BX09] and theoretically explained in [PTK11]. In fact, the Runge phenomenon occurs for any points that do not cluster quadratically like Chebyshev points, see [APS19]. It is however defeated by interpolation on Chebyshev type points.

As far as interpolation in spaces $V_n$ is concerned, there exists no systematic choice of points that prevents all instabilities. A generic choice is that of Fekete points, which ensure that $\|P_n^n g\|_{\mathcal{L}^\infty} \le n\|g\|_{\mathcal{L}^\infty}$ in the strong uniform norm $\|g\|_{\mathcal{L}^\infty} = \sup_{x \in X} |g(x)|^1$, resulting in the stability inequality

$$ \|f - P_n^n f\|_{\mathcal{L}^\infty} \le (n+1) \inf_{v \in V_n} \|f - v\|_{\mathcal{L}^\infty}, $$

see for instance Proposition 1.2.5 in [Nov88]. This guarantees a good approximation $P_n^n f$ of $f$, provided the latter is sufficiently smooth and $V_n$ is well chosen. However, for general domains $X$ and spaces $V_n$, the computation of Fekete points can be intractable.

For polynomial interpolation over compact domains $X \subset \mathbb{R}$ or $\mathbb{C}$, Leja sequences are greedy alternatives to Fekete points. For instance, for interpolation by polynomials in $\mathbb{P}_{n-1}$, restricted to union of closed intervals of $\mathbb{R}$, they provably yield $\|P_n^n g\|_{\mathcal{L}^\infty} \lesssim n^{13/4}\|g\|_{\mathcal{L}^\infty}$, see the recent paper [AN22]. We note that numerical evidence shows that $n^{13/4}$ can be replaced by $n$ for intervals. Similar polynomial growths, with a factor $n^3$ in the estimate, hold for a certain type of polynomial interpolation over tensor product domains, see [CCS14, CC15].

On the other hand, near optimal approximation error in expected $L^2$ sense

$$ \mathbb{E}(\|f - P_n^m f\|_{L^2}^2) \le C \inf_{v \in V_n} \|f - v\|_{L^2}^2, $$

can be attained by taking $m$ larger than $n$. In the case of uniformly distributed points and equal weights $s_i = 1$, this usually requires $m$ to scale polynomially in $n$, see [CDL13, CCM$^+$15], but a logarithmic oversampling is achievable if one considers a different sample distribution [CM17]. Numerical methods [HNP22] and a theoretical solution [CD22] have been proposed to reduce the sample size $m$ to a constant multiple of $n$. A discrete version of the above bound can also be found in [CP19], in the context of statistical machine learning.

In the last five papers, the sample points are drawn at random according to a prescribed measure, and the error bounds are presented in expectation. This setting does not require any additional assumption: indeed, point evaluations of a function $f \in L^2$ are defined in an almost sure sense, and so is the weighted least-squares projection $P_n^m f$. This will also be the main framework for the present article.

<footnote n="1">¹We assume here that $V_n$ is included in $\mathcal{L}^\infty$, the space of bounded functions.</footnote>
Our main theorem, stated below, provides new bounds on the $L^2$ approximation error, depending on the ratio between $m$ and $n$.

**Theorem 1.1.** Let $m \ge n$.

* The conditional weighted least-squares approximation $\tilde{f} \in V_n$ defined in (15), using $m$ evaluations of $f$ at points selected by Algorithm 3.1, with inputs $\varepsilon = r^{-1/4}$ and $\gamma = r^{1/2} - r^{1/4}$, satisfies

$$ \mathbb{E}(\|f - \tilde{f}\|_{L^2}^2) \le \left(1 + \frac{1}{r(1 - r^{-1/4})^7}\right) \min_{v \in V_n} \|f - v\|_{L^2}^2, \quad (1) $$

where $r = \frac{m+1}{n} > 1$.

* In turn, the weighted least-squares estimator $\tilde{f} = P_n^m f$, based on points selected by Algorithm 4.1, with inputs $\delta = r^{-1/2}$ and $\kappa \in [0, 1]$, simultaneously satisfies

$$ \mathbb{E}(\|f - \tilde{f}\|_{L^2}^2) \le \left(1 + \frac{1}{1 - \kappa(1 - 1/\sqrt{r})^2}\right) \min_{v \in V_n} \|f - v\|_{L^2}^2 \quad (2) $$

if $\kappa < 1$ and

$$ \|f - \tilde{f}\|_{L^2} \le \left(1 + \frac{1}{\sqrt{\kappa}(1 - 1/\sqrt{r})}\right) \min_{v \in V_n} \|f - v\|_{L^\infty} \quad \text{a.s.} \quad (3) $$

if $\kappa > 0$. Here $r = \frac{m}{n-1} > 1$, and we assume that $V_n \subset L^\infty$ in the last bound.

The first two bounds are of the same kind, the main difference being in the constant factor on the right-hand side, which is of optimal order $1 + 1/r + o(1/r)$ in (1) as $r$ tends to infinity, but behaves better in (2) when $r$ gets close to 1. The parameter $\kappa$ is a free input parameter selected by the user. The choice $\kappa = 0$ gives the best constant in (2), while $\kappa = 1$ is best suited for (3). One can also construct an estimator $\tilde{f} = P_n^m f$ satisfying both estimates at the same time, by running Algorithm 4.1 with an intermediate value of $\kappa$.

The uniform bound (3) follows the approach developed in [LT22, Tem21, PU22, DT24, BSU23], slightly improving the constants when compared to Theorem 1.1 in [Tem21] and Theorem 6.3 in [BSU23], and linking it to the approach in expectation (2) by the use of a common algorithm. In a third approach, similar deterministic bounds can be proved by assuming more regularity on $f$ through a nested sequence of approximation spaces $(V_n)_{n \in \mathbb{N}}$, in a Hilbert space setting [KU21a, KUV21, MU21, NSU22, BSU23, GW24] or in more general Banach spaces [KU21b, DKU23, KPUU23].

Taking $m = n$ in the last two estimates, and observing that $1/\sqrt{r} = \sqrt{1 - 1/n} \le 1 - 1/2n$, we immediately obtain the following result on interpolation in $L^2$.

**Corollary 1.2.** For $m = n$, the interpolation $\tilde{f} = P_n^m f \in V_n$ of $f$ at random points $x_1, \dots, x_n$ selected by Algorithm 4.1 with input $\delta = \sqrt{1 - 1/n}$ achieves the accuracy bounds

$$ \mathbb{E}(\|f - \tilde{f}\|_{L^2}^2) \le (4n^2 + 1) \min_{v \in V_n} \|f - v\|_{L^2}^2 \quad \text{if } \kappa = 0, $$

and, assuming $V_n \subset L^\infty$,

$$ \|f - \tilde{f}\|_{L^2} \le (2n + 1) \min_{v \in V_n} \|f - v\|_{L^\infty} \quad \text{a.s.} \quad \text{if } \kappa = 1. $$
**Remark 1.3.** In the case of uniform approximation, if each function in $V_n$ is uniformly bounded on $X$ and not just essentially bounded, one can replace $L^\infty$ by the strong supremum norm $\mathcal{L}^\infty$, and remove the almost sure restriction, by considering deterministic samples satisfying the constraints of our algorithms. A similar observation can be found in Remark 3 of [KPUU23].

Finally, one can combine the second estimate in Corollary 1.2 with an inverse inequality between $L^2$ and $\mathcal{L}^\infty$ in $V_n$.

**Corollary 1.4.** If the functions in $V_n$ are bounded and $X$ is compact, Algorithm 4.1, applied with the optimal measure $\mu^*$ from [KW60], provides $n$ points for which the interpolation $\tilde{f}^* = P_n^nf \in V_n$ of $f$ satisfies

$$ \|f - \tilde{f}^*\|_{\mathcal{L}^\infty} \le (1 + 2n\sqrt{n}) \min_{v \in V_n} \|f - v\|_{\mathcal{L}^\infty}. $$

Note that in the case $m = n$, the values of the weights $s_i$ have no importance, since the minimum in the definition of $P_n^m$ is zero. In fact, if measure $\mu^*$ is known, we exhibit a constructive set of points such that the Lebesgue stability constant

$$ \mathbb{L}_n = \max_{f \in C^0(X) \setminus \{0\}} \frac{\|P_n^nf\|_{\mathcal{L}^\infty}}{\|f\|_{\mathcal{L}^\infty}} $$

is at most $2n\sqrt{n}$. Although Fekete points achieve $\mathbb{L}_n \le n$, their computational complexity is exponential in $n$. In our approach, the main computational challenge is to find an optimal measure $\mu^*$ as defined in [KW60], see also [Bos90], which is in general as difficult as looking for the Fekete points.

However, using a different measure $\mu$, we still obtain a result similar to Corollary 1.4 however with the factor $\sqrt{n}$ replaced by a larger power of $n$. For example, if $X = [-1, 1]^d$ and $V_n = \text{span}\{x^\nu : \nu \in \Lambda\}$ is a space of multivariate polynomials indexed by a lower set $\Lambda \subset \mathbb{N}_0^d$, the uniform measure yields a factor $n$ instead of $\sqrt{n}$ and an estimate $\mathbb{L}_n \le 2n^2$, while the tensor product arcsine measure achieves a factor $\min(2^{d/2}\sqrt{n}, n^{\log 3/\log 4})$ instead of $\sqrt{n}$, resulting in the estimate $\mathbb{L}_n \le 2n \min(2^{d/2}\sqrt{n}, n^{\log 3/\log 4})$.

Note that the Lebesgue constant does not only determine the convergence of the approximation, but also its robustness to numerical errors and noise in the measurements, see Section 4 in [CM17], as well as [APS19, PTK11].

The above discussion demonstrates that good interpolation points can be constructed for multivariate polynomial approximation. The estimate $2n^{1+\log 3/\log 4}$ on the Lebesgue constant is moderate. It outperforms the estimate $n^3$ which was established in [CC15] using constructions based on $\mathbb{R}$-Leja sequences. We note however that the new procedure is not hierarchical.

**Remark 1.5.** One can also use an inverse inequality between $L^2$ and $\mathcal{L}^\infty$ in the least-squares regime $m > n$. More precisely, combining (3) with [KW60] as in Corollary 1.4, we obtain

$$ \|f - P_n^m f\|_{\mathcal{L}^\infty} \le \left(1 + \frac{\sqrt{n}}{1 - 1/\sqrt{r}}\right) \min_{v \in V_n} \|f - v\|_{\mathcal{L}^\infty}, $$

with $r = m/(n-1)$. Taking the supremum for $f$ in a compact class of functions, and optimizing on both sides over $V_n$, this implies that the uniform sampling numbers are bounded by $O(\sqrt{n})$ times the uniform Kolmogorov $n$-widths, improving on the $O(n)$ factor of Fekete points, if we allow for a constant oversampling $r > 1$. This can already be seen by applying Corollary 5.3 in [PU22] together with the optimal density from [KW60]. The same result is obtained in Theorem 3 of [KPUU23], where error bounds are also derived in any $L^p$ norm. We additionally refer to [GW24], where implications in the field of Information Based Complexity are drawn,
in a specific Hilbert space setting. The two very recent papers [KPUU23, GW24] rely on the pioneering works [KU21a, KU21b], and on the infinite-dimensional adaptation [DKU23] of the result from [MSS15], see also [FS19].

The rest of the paper is organized as follows. In Section 2, we inspect the weighted least-squares projection, and point out why different strategies should be used when $r \gg 1$ and $r \approx 1$. Sections 3 and 4 introduce and analyze the sampling Algorithms 3.1 and 4.1. They are independent from one another, up to the shared use of a few formulas, and the various estimates of Theorem 1.1 are proved separately. Finally, we discuss some numerical aspects of the presented algorithms in Section 5, and provide numerical illustrations in Section 6.

## 2 Least-squares

Let $\varphi = (\varphi_1, \dots, \varphi_n)$ be an orthonormal basis of $V_n$ in $L^2(X, \mu)$. For any $x \in X$, we consider $\varphi(x)$ as a vector in $\mathbb{C}^n$, denote $\varphi(x)^*$ its conjugate transpose, and $|\varphi(x)|^2 = \varphi(x)^*\varphi(x)$ its squared euclidian norm. Observe that, given any $n \times n$ matrix $M$,

$$ \int_X \varphi(x)^* M \varphi(x) d\mu(x) = \mathrm{Tr} \left( M \int_X \varphi(x) \varphi(x)^* d\mu(x) \right) = \mathrm{Tr}(M). \quad (4) $$

Adopting the formalism from [KUV21, MU21, NSU22]), we let $\mathbf{f} = (\sqrt{s_i}f(x_i))_{1\le i \le m} \in \mathbb{C}^m$ be the measurement vector, and

$$ \mathbf{L} = (\sqrt{s_i} \varphi_j(x_i))_{1\le i \le m, 1\le j \le n} \in \mathbb{C}^{m \times n} $$

the collocation matrix. Then the weighted least-square estimator writes

$$ P_n^m f = \sum_{j=1}^{n} a_j \varphi_j, \quad a = \mathbf{L}^+ \mathbf{f}, $$

where $\mathbf{L}^+ = (\mathbf{L}^* \mathbf{L})^{-1} \mathbf{L}^*$ stands for the Moore-Penrose pseudo-inverse of $\mathbf{L}$. In the sequel, we will make sure that $\mathbf{L}$ has full column rank $n$, so that the Gram matrix

$$ \mathbf{A}_m := \mathbf{L}^* \mathbf{L} = \sum_{i=1}^{m} s_i \varphi(x_i) \varphi(x_i)^* $$

is indeed invertible. Denoting $g = f - P_n f$ the optimal residual error and $\mathbf{g} = (\sqrt{s_i}g(x_i))_{1\le i \le m}$ the associated vector, the least-squares error decomposes as

$$ \|f - P_n^m f\|_{L^2}^2 = \|f - P_n f\|_{L^2}^2 + \|P_n f - P_n^m f\|_{L^2}^2 = \|g\|_{L^2}^2 + \|P_n^m g\|_{L^2}^2. \quad (5) $$

There are two possible strategies for bounding $\|P_n^m g\|_{L^2}^2$: either we use

$$ \|P_n^m g\|_{L^2}^2 = |\mathbf{A}_m^{-1} \mathbf{L}^* \mathbf{g}|^2 \le \|\mathbf{A}_m^{-1}\|_{2 \to 2}^2 |\mathbf{L}^* \mathbf{g}|^2 = \lambda_{\min}(\mathbf{A}_m)^{-2} |\langle \varphi, g \rangle_m|^2, \quad (6) $$

where $\langle \varphi, g \rangle_m := (\langle \varphi_j, g \rangle_m)_{1\le j \le n} \in \mathbb{C}^n$, or we bound it by

$$ \|P_n^m g\|_{L^2}^2 = |\mathbf{L}^+ \mathbf{g}|^2 \le \|\mathbf{L}^+\|_{2 \to 2}^2 |\mathbf{g}|^2 = \|\mathbf{A}_m^{-1}\|_{2 \to 2}^2 \|g\|_m^2 = \lambda_{\min}(\mathbf{A}_m)^{-1} \|g\|_m^2. \quad (7) $$

The first approach is expected to give better estimates when $m$ is much larger than $n$, since in that case the discrete inner product $\langle \cdot, \cdot \rangle_m$ weakly converges to the continuous inner product $\langle \cdot, \cdot \rangle_{L^2}$, and as $g$ is orthogonal to $V_n$, the vector $\langle \varphi, g \rangle_m$ should be small. On the other hand,
when $m$ is close to $n$, $L$ may be ill-conditioned, leading to small values of $\lambda_{\min}(A_m)$, thus favouring the second approach.

In both situations, one should choose the points $x_i$ and weights $s_i$ in order to control the smallest eigenvalue of $A_m$ from below. This can be ensured through a greedy selection of points, based on the effective resistance [SS11]

$$ \varphi(x)^*(A - \ell I)^{-1}\varphi(x) $$

of a point $x \in X$ with respect to a hermitian matrix $A$, with $\ell$ a lower bound on the eigenvalues of $A$, and $I$ the identity matrix. Intuitively, if the points $x_1, \dots, x_{i-1}$ are already fixed, the effective resistance of $x_i$ with respect to the partial Gram matrix

$$ A_{i-1} = \sum_{i=1}^{i-1} s_i \varphi(x_i) \varphi(x_i)^* $$

quantifies how close $\varphi(x_i)$ is to the eigenvectors of $A_{i-1}$ with small eigenvalues, helping the lower potential $\mathrm{Tr}((A - \ell I)^{-1})$ to decrease when going from $A_{i-1}$ to $A_i$. The key idea in [BSS09] is to increase the lower barrier $\ell$ at each step, without increasing the potential more than it had decreased when adding sample $x_i$. Keeping the potential bounded ensures that all updates of $\ell$ are of the same size $\delta$, and its trace form is chosen to allow rank-one matrices $s_i\varphi(x_i)\varphi(x_i)^*$ to decrease the potential, independently of the distribution of eigenvalues of $A_{i-1}$. In this way, we obtain the desired bound $\lambda_{\min}(A_m) > \ell_m$ at the end of the algorithm.

As in [BSS09], it is also possible to ensure an upper bound $\lambda_{\max}(A_m) < u_m$ on the eigenvalues of $A_m$, by considering the upper potential $\mathrm{Tr}((uI - A)^{-1})$. Although $\lambda_{\max}(A_m)$ does not explicitly appear in the least-squares formulation, it is worth noticing that the ratio $u_m/\ell_m$ controls the condition number $\lambda_{\max}(A_m)/\lambda_{\min}(A_m)$ of $A_m$, which determines the computational cost of solving the least-squares problem via an iterative solver, as well as its robustness to numerical error. This is discussed in, e.g., Section 5.3.4 of [ABW22].

Nevertheless, we chose not to include this upper potential, because it would have prevented us from achieving the critical sampling budget $m \approx n$. Heuristically, when we use both upper and lower potentials, the vector $\varphi(x_i)$ should increase the smallest eigenvalues of $A_{i-1}$, while staying far from the eigenvectors with largest eigenvalues. Due to the second condition, the lower barrier can only increase half as much as it could without it, resulting in the constraint $m \ge 2n$. In addition, we still achieve an a posteriori upper bound, see (21), which is not particularly sharp, but sufficient for our purposes. Lastly, a full inversion of matrix $A_m$, with complexity $O(n^3)$, remains reasonable in view of the moderate values of $n$ imposed by the sampling costs, see Section 5.

### 3 Random sampling by effective resistance

We first consider an approach combining (5) and (6) to obtain the first estimate (1) in Theorem 1.1. To observe the appropriate decay as $r \to \infty$, we impose that the vector

$$ \langle \varphi, g \rangle_m = \sum_{i=1}^{m} s_i \varphi(x_i) \overline{g(x_i)} $$

is an unbiased estimator of $\langle \varphi, g \rangle_L^2 = 0$, by taking weights $s_i$ inversely proportional to the sampling density of $x_i$. If the points $x_i$ are independent, the euclidian norm $|\langle \varphi, g \rangle_m|^2$ is therefore
bounded in expectation by the sum of the variances of each term, of the form

$$ \mathbb{E}(s_i^2 | \varphi(x_i)|^2 | g(x_i)|^2) \propto \int_X s_i |\varphi(x_i)|^2 |g(x_i)|^2 d\mu(x_i) \le \|s_i|\varphi(x_i)|^2\|_{L^\infty} \|g\|_{L^2}^2. $$

As a result, to control $\ell_m$ from below and these variances from above, the appropriate sampling density for point $x_i$ is a combination of the effective resistance mentioned earlier and the so-called Christoffel function $|\varphi(x)|^2 = \sum_{j=1}^n |\varphi_j(x)|^2$.

Algorithm 3.1 is inspired by the deterministic procedure from [BSS09] and its randomization presented in [LS18], which has already been applied to least-squares recovery in [CP19]. It takes as inputs the parameters

$$ \varepsilon \in (0, 1), \quad \eta = \frac{\varepsilon}{1 - \varepsilon} \in (0, \infty), \quad \text{and} \quad \gamma \ge 0, $$

which respectively influence the size of the barrier increments $\delta_i = \ell_i - \ell_{i-1}$, the weights $s_i$, and the balance between effective resistance and Christoffel function in the sampling density.

**Algorithm 3.1:** Random sampling by effective resistance

**Input:** probability measure $\mu$, parameter $\varepsilon \in (0, 1)$, and $\gamma \ge 0$.

$$ \mathbf{A}_0 = \mathbf{0} \in \mathbb{R}^{n \times n}, \quad \ell_0 = -n $$

**for** $i = 1, \dots, m$ **do**

Let $\mathbf{Y}_i = (\mathbf{A}_{i-1} - \ell_{i-1} \mathbf{I})^{-1}$

Update $\ell_i = \ell_{i-1} + \delta_i$, where $\delta_i := \varepsilon/\Psi_i$ with $\Psi_i = \text{Tr}(\mathbf{Y}_i) + \gamma$

Let $\mathbf{Z}_i = (\mathbf{A}_{i-1} - \ell_i \mathbf{I})^{-1}$

Let $\rho_i : x \mapsto \varphi(x)^*(\mathbf{Z}_i + \gamma \mathbf{I}/n)\varphi(x)$

Draw $x_i$ from probability measure $\frac{\rho_i(x)}{\Xi_i} d\mu(x)$, where $\Xi_i = \text{Tr}(\mathbf{Z}_i) + \gamma$

Let $s_i = \eta/\rho_i(x_i)$

Update $\mathbf{A}_i = \mathbf{A}_{i-1} + s_i \varphi(x_i) \varphi(x_i)^*$

**end for**

**Output:** sample $\{x_1, \dots, x_m\}$, weights $\{s_1, \dots, s_m\}$, and matrix $\mathbf{A}_m$.

For the sake of the analysis, we also define

$$ \mathbf{Y}_{m+1} = (\mathbf{A}_m - \ell_m \mathbf{I})^{-1}, \quad \Psi_{m+1} = \text{Tr}(\mathbf{Y}_{m+1}) + \gamma, \\ \delta_{m+1} = \frac{\varepsilon}{\Psi_{m+1}} \quad \text{and} \quad \ell_{m+1} = \ell_m + \delta_{m+1}, $$

as would have been done at iteration $i = m + 1$.

**Lemma 3.1.** The algorithm is almost surely well defined and $\mathbf{A}_i \succ \ell_{i+1} \mathbf{I}$ for all $0 \le i \le m$, where $\succ$ stands for the Loewner order between positive semi-definite matrices.

**Proof.** We proceed by induction on $i$. At initialization $\mathbf{A}_0 \succ -n \mathbf{I} = \ell_0 \mathbf{I}$. For $i \in \{1, \dots, m\}$ fixed, assume that $\mathbf{A}_{i-1} \succ \ell_{i-1} \mathbf{I}$. This implies that $\mathbf{Y}_i$ is well-defined, positive definite, and

$$ \delta_i \mathbf{I} = \frac{\varepsilon \mathbf{I}}{\text{Tr}(\mathbf{Y}_i) + \gamma} \le \frac{\varepsilon \mathbf{I}}{\lambda_{\max}(\mathbf{Y}_i)} = \lambda_{\min}(\mathbf{A}_{i-1} - \ell_{i-1} \mathbf{I}) \varepsilon \mathbf{I} \prec \mathbf{A}_{i-1} - \ell_{i-1} \mathbf{I}. \quad (8) $$

Hence $\mathbf{A}_{i-1} \succ (\ell_{i-1} + \delta_i) \mathbf{I} = \ell_i \mathbf{I}$, so $\mathbf{Z}_i$ is well-defined, positive definite, and by identity (4),

$$ \int_X \rho_i(x) d\mu(x) = \text{Tr}(\mathbf{Z}_i + \gamma \mathbf{I}/n) = \Xi_i, $$

proving that $\rho_i d\mu/\Xi_i$ is indeed a probability density. Finally $s_i > 0$ almost surely, so $\mathbf{A}_i \succ \mathbf{A}_{i-1} \succ \ell_i \mathbf{I}$. This completes the induction, and a last iteration shows that $\mathbf{A}_m \succ \ell_{m+1} \mathbf{I}$, concluding the proof. □
Contrarily to most variations on the algorithm from [BSS09], no upper potential is used to bound the eigenvalues of $\mathbf{A}_i$ from above. Instead, parameter $\gamma$ provides a lower bound on $\rho_i$ in terms of the Christoffel function $\varphi(x)^*\varphi(x)$, which turns into an upper bound on the norm of rank-one terms $s_i\varphi(x_i)\varphi(x_i)^*$, and therefore on the norm of their sum $\mathbf{A}_m$.

To bound the eigenvalues of $\mathbf{A}_i$ from below, we use the lower barrier $\ell_{i+1}$, which should increase at each step, at a speed $\delta_i$ controlled by the lower potential $\mathrm{Tr}(\mathbf{Y}_i)$. As the densities $\rho_i$ are positive, the sample $(x_1, \dots, x_m)$ could be in any part of $X^m$, so there cannot be any positive deterministic lower bound on $\lambda_{\min}(\mathbf{A}_m)$. However, the following lemma proves a monotonicity property in expectation, similar to Lemmas 4.4 in [LS18] and [LS17], themselves inspired by Lemmas 3.3 and 3.4 in [BSS09].

**Lemma 3.2.** The sequence $(\mathbb{E}(\mathrm{Tr}(\mathbf{Y}_i)))_{1\le i\le m+1}$ is non-increasing.

*Proof.* Let $1 \le i \le m$, and denote $u_i = \sqrt{s_i}\varphi(x_i)$. Applying the Sherman-Morrison formula yields

$$ \mathbf{Y}_{i+1} = (\mathbf{Z}_i^{-1} + u_i u_i^*)^{-1} = \mathbf{Z}_i - \frac{\mathbf{Z}_i u_i u_i^* \mathbf{Z}_i}{1 + u_i^* \mathbf{Z}_i u_i}. \quad (9) $$

By design of the sampling density and weights,

$$ u_i^* \mathbf{Z}_i u_i = s_i \varphi(x_i)^* \mathbf{Z}_i \varphi(x_i) \le s_i \rho_i(x_i) = \eta, \quad (10) $$

so by definition of $\eta = \varepsilon/(1-\varepsilon)$, which implies $1+\eta = \eta/\varepsilon$, we obtain

$$ \mathbf{Y}_{i+1} - \mathbf{Z}_i \asymp - \frac{\mathbf{Z}_i u_i u_i^* \mathbf{Z}_i}{1 + \eta} = - \frac{\varepsilon}{\eta} s_i \mathbf{Z}_i \varphi(x_i) \varphi(x_i)^* \mathbf{Z}_i. \quad (11) $$

On the other hand, $\mathbf{Y}_i$ and $\mathbf{Z}_i$ commute because their inverses do, and we can write

$$ \mathbf{Z}_i - \mathbf{Y}_i = (\mathbf{Y}_i^{-1} - \mathbf{Z}_i^{-1}) \mathbf{Y}_i \mathbf{Z}_i = \delta_i \mathbf{Y}_i \mathbf{Z}_i. $$

Combining this identity to the previous estimate and taking traces gives

$$ \mathrm{Tr}(\mathbf{Y}_{i+1}) - \mathrm{Tr}(\mathbf{Y}_i) = \delta_i \mathrm{Tr}(\mathbf{Y}_i \mathbf{Z}_i) - \frac{\varepsilon}{\eta} s_i \varphi(x_i)^* \mathbf{Z}_i^2 \varphi(x_i). $$

Observe that $\mathbf{Y}_i$ and $\mathbf{Z}_i$ only depend on the samples $x_1, \dots, x_{i-1}$. As a result, if these samples are fixed, taking the expectation $\mathbb{E}_{x_i} = \mathbb{E}(\cdot|x_1, \dots, x_{i-1})$ with respect to point $x_i$, we arrive at

$$ \begin{aligned} \mathbb{E}_{x_i}(\mathrm{Tr}(\mathbf{Y}_{i+1})) - \mathrm{Tr}(\mathbf{Y}_i) &= \delta_i \mathrm{Tr}(\mathbf{Y}_i \mathbf{Z}_i) - \frac{\varepsilon}{\eta} \int_X \frac{\eta}{\rho_i(x_i)} \varphi(x_i)^* \mathbf{Z}_i^2 \varphi(x_i) \frac{\rho_i(x_i)}{\Xi_i} d\mu(x_i) \\ &= \varepsilon \left( \frac{\mathrm{Tr}(\mathbf{Y}_i \mathbf{Z}_i)}{\mathrm{Tr}(\mathbf{Y}_i) + \gamma} - \frac{\mathrm{Tr}(\mathbf{Z}_i^2)}{\mathrm{Tr}(\mathbf{Z}_i) + \gamma} \right), \end{aligned} \quad (12) $$

where we have used identity (4). Since

$$ \mathbf{Z}_i = \mathbf{Y}_i + \delta_i \mathbf{Y}_i \mathbf{Z}_i \quad \text{and} \quad \mathbf{Z}_i^2 = \mathbf{Y}_i \mathbf{Z}_i + \delta_i \mathbf{Y}_i \mathbf{Z}_i^2, \quad (13) $$

we have by Cauchy-Schwarz inequality

$$ \mathrm{Tr}(\mathbf{Y}_i \mathbf{Z}_i) \mathrm{Tr}(\mathbf{Z}_i) - \mathrm{Tr}(\mathbf{Y}_i) \mathrm{Tr}(\mathbf{Z}_i^2) = \delta_i (\mathrm{Tr}(\mathbf{Y}_i \mathbf{Z}_i)^2 - \mathrm{Tr}(\mathbf{Y}_i) \mathrm{Tr}(\mathbf{Y}_i \mathbf{Z}_i^2)) \le 0. \quad (14) $$

Moreover, $\mathrm{Tr}(\mathbf{Y}_i \mathbf{Z}_i) - \mathrm{Tr}(\mathbf{Z}_i^2) = -\delta_i \mathrm{Tr}(\mathbf{Y}_i \mathbf{Z}_i^2) = -\delta_i \mathrm{Tr}(\mathbf{Z}_i \mathbf{Y}_i \mathbf{Z}_i) < 0$. Hence the right-hand side in (12) is negative. Taking an expectation over the previous samples $x_1, \dots, x_{i-1}$ in turn implies that $\mathbb{E}(\mathrm{Tr}(\mathbf{Y}_{i+1}))$ is smaller than $\mathbb{E}(\mathrm{Tr}(\mathbf{Y}_i))$, which concludes the proof. $\square$
One main difference with [LS18] and [LS17] is the choice of different matrices $\mathbf{Y}_i$ and $\mathbf{Z}_i$ for updating the lower barrier $\ell_i$ and selecting a new point $x_i$. This comes at the expense of the refined analysis (13), (14) to show that (12) is non-positive. Moreover, in view of (8), this changes the sampling density at most by a factor $1/(1 - \varepsilon)$.

However, if we had used $\mathbf{Y}_i$ instead of $\mathbf{Z}_i$ in the sampling density $\rho_i$, in other words if we had defined $\rho_i(x) = \varphi(x)^*(\mathbf{Y}_i + \gamma\mathbf{I}/n)\varphi(x)$ and $\Xi_i = \text{Tr}(\mathbf{Y}_i) + \gamma$, we would have obtained

$$u_i^* \mathbf{Z}_i u_i \le u_i^* \frac{\mathbf{Y}_i}{1 - \varepsilon} u_i \le \frac{s_i \rho_i(x_i)}{1 - \varepsilon} = \frac{\eta}{1 - \varepsilon}$$

instead of (10), where the first inequality comes from the fact that $\mathbf{Y}_i^{-1} - \mathbf{Z}_i^{-1} \le \varepsilon \mathbf{Y}_i^{-1}$ in view of (8). In order to recover the same right-hand side as before in inequality (11), one would then need to assume

$$\frac{1}{1 + \frac{\eta}{1 - \varepsilon}} \ge \frac{\varepsilon}{\eta} \iff 1 + \frac{\eta}{1 - \varepsilon} \le \frac{\eta}{\varepsilon} \iff \eta \ge \frac{1}{\frac{1}{\varepsilon} - \frac{1}{1 - \varepsilon}} = \frac{\varepsilon(1 - \varepsilon)}{1 - 2\varepsilon}.$$

So this simpler approach requires $\varepsilon < \frac{1}{2}$, which is not possible unless $r > 2$, and even in that case our approach gives slightly better estimates.

The control on the lower potentials in expectation from Lemma 3.2 yields a lower bound on the eigenvalues of $\mathbf{A}_m$ in probability.

**Proposition 3.3.** Let $m \ge n$, $p \in (0, 1)$, and define $r = (m + 1)/n$. The random matrix $\mathbf{A}_m$ generated by Algorithm 3.1 satisfies

$$\mathbf{A}_m \succsim \alpha \mathbf{I}, \quad \alpha = n \left( \frac{\varepsilon r}{1/p + \gamma} - 1 \right),$$

with probability at least $1 - p$.

*Proof.* In view of the monotonicity property established in Lemma 3.2, and thanks to the initialization $\ell_0 = -n$, we have

$$\mathbb{E}(\text{Tr}(\mathbf{Y}_i)) \le \text{Tr}(\mathbf{Y}_1) = 1, \quad 1 \le i \le m + 1.$$

By Markov's inequality

$$\mathbb{P}\left(\sum_{i=1}^{m+1} \text{Tr}(\mathbf{Y}_i) > \frac{m+1}{p}\right) \le p \frac{\mathbb{E}\left(\sum_{i=1}^{m+1} \text{Tr}(\mathbf{Y}_i)\right)}{m+1} \le p.$$

As $\Psi_i = \text{Tr}(\mathbf{Y}_i) + \gamma$, we deduce that with probability at least $1 - p$, it holds

$$\sum_{i=1}^{m+1} \Psi_i \le (m+1) \left(\frac{1}{p} + \gamma\right).$$

Together with the inequality $(\sum_{i=1}^{m+1} \Psi_i) (\sum_{i=1}^{m+1} \Psi_i^{-1}) \ge (m+1)^2$, this implies

$$\sum_{i=1}^{m+1} \Psi_i^{-1} \ge \frac{m+1}{1/p + \gamma} = \frac{rn}{1/p + \gamma},$$

and therefore

$$\ell_{m+1} = \ell_0 + \sum_{i=1}^{m+1} \frac{\varepsilon}{\Psi_i} \ge -n + \frac{rn\varepsilon}{1/p + \gamma} = \alpha,$$

with probability at least $1 - p$. Since $\mathbf{A}_m \succ \ell_{m+1} \mathbf{I}$, see Lemma 3.1, the proof is complete. $\square$
**Remark 3.4.** The lower bound $\mathbf{A}_m \succcurlyeq \alpha I$ can be rewritten as a frame inequality

$$c^* \mathbf{A}_m c = \sum_{i=1}^{m} s_i |\langle c, \varphi(x_i) \rangle|^2 \geq \alpha |c|^2, \quad c \in \mathbb{C}^n,$$

or as a Marcinkiewicz-Zygmund inequality: for all $v \in V_n$, $\|v\|_m^2 \geq \alpha \|v\|_{L^2}^2$. All three versions have been extensively used in the literature for emphasizing the relations with subsampling of frames and discretization of continuous norms, see [FS19, NOU13, NOU16, LT22]. Note that our approach is not performing a subsampling in general, although it can be interpreted as such in the context of Remark 5.1.

We are now ready for the proof of the first statement of Theorem 1.1. Let $E$ be the event of probability at least $1-p$ where $\mathbf{A}_m \succcurlyeq \alpha \mathbf{I}$, for $p \in (0, 1)$ and $\alpha > 0$ given by Proposition 3.3, and define

$$\tilde{f} = P_n^m f |E \qquad (15)$$

as the weighted least-squares projection, conditioned to event $E$. In practice, this amounts to relaunching Algorithm 3.1 until $E$ occurs, and computing the weighted least-squares projection of $f$ with the last sample.

*Proof of Theorem 1.1, equation (1).* As $V_n \subset L^2$, we can assume that $f$ is also in $L^2$, otherwise the right-hand side would be infinite. Similar to (5), by Pythagoras theorem,

$$\begin{aligned} \mathbb{E}(\|f - \tilde{f}\|_{L^2}^2) &= \|f - P_n f\|_{L^2}^2 + \mathbb{E}(\|\tilde{f} - P_n f\|_{L^2}^2) \\ &= \|g\|_{L^2}^2 + \mathbb{E}(\|P_n^m f - P_n f\|_{L^2}^2 | E) \\ &= \|g\|_{L^2}^2 + \mathbb{E}(\|P_n^m g\|_{L^2}^2 | E). \end{aligned}$$

where $g = f - P_n f$. Combining this with equation (6), together with the definition of $E$, yields

$$\mathbb{E}(\|f - \tilde{f}\|_{L^2}^2) \leq \|g\|_{L^2}^2 + \alpha^{-2} \mathbb{E}(|\langle \varphi, g \rangle_m|^2 | E).$$

We now use Proposition 3.3:

$$\mathbb{E}(|\langle \varphi, g \rangle_m|^2 | E) = \frac{\mathbb{E}(|\langle \varphi, g \rangle_m|^2 \chi_E)}{\mathbb{P}(E)} \leq \frac{1}{1-p} \mathbb{E}(|\langle \varphi, g \rangle_m|^2).$$

Develop the discrete inner products $\langle \cdot, \cdot \rangle_m$ as

$$\mathbb{E}(|\langle \varphi, g \rangle_m|^2) = \mathbb{E} \left( \sum_{k=1}^{n} |\langle \varphi_k, g \rangle_m|^2 \right) = \sum_{k=1}^{n} \mathbb{E} \left( \sum_{i=1}^{m} \overline{s_i \varphi_k(x_i)} g(x_i) \sum_{j=1}^{m} s_j \varphi_k(x_j) \overline{g(x_j)} \right).$$

When $i > j$, we have that

$$\mathbb{E}_{x_i} \left( s_i \overline{\varphi_k(x_i)} g(x_i) \right) = \frac{\eta}{\Xi_i} \int_X \overline{\varphi_k(x)} g(x) d\mu(x) = 0,$$

since $g$ is orthogonal to $V_n$, so the corresponding term is zero. The same holds true for $i < j$, by exchanging $i$ and $j$. We are only left with the diagonal terms

$$\mathbb{E}(|\langle \varphi, g \rangle_m|^2) = \mathbb{E} \left( \sum_{i=1}^{m} s_i^2 |\varphi(x_i)|^2 |g(x_i)|^2 \right) \leq \frac{nn\eta}{\gamma} \mathbb{E} \left( \sum_{i=1}^{m} s_i |g(x_i)|^2 \right) = \frac{nn\eta}{\gamma \Xi_i} \|g\|_{L^2}^2,$$
where we have first used the uniform bound $s_i|\varphi(x_i)|^2 \le n\eta/\gamma$, then integrated over $x_i$ using $s_i\rho(x_i) = \eta$, both identities following from the definition of $s_i$ in Algorithm 3.1. Since $\Xi_i \ge \gamma$, $m \le rn$ and $\eta = \varepsilon/(1-\varepsilon)$, we deduce that

$$ \begin{aligned} \mathbb{E}(\|f - \tilde{f}\|_{L^2}^2) &\le \left(1 + \frac{1}{\alpha^2(1-p)} \frac{rn^2\varepsilon^2}{\gamma^2(1-\varepsilon)^2}\right) \|g\|_{L^2}^2 \\ &= \left(1 + \frac{r}{1-p} \frac{1}{\beta(\varepsilon, \gamma)^2}\right) \|g\|_{L^2}^2, \end{aligned} \quad (16) $$

where

$$ \beta(\varepsilon, \gamma) = \frac{\alpha\gamma(1-\varepsilon)}{n\varepsilon} = \gamma \left( \frac{r\varepsilon}{1/p + \gamma} - 1 \right) \left( \frac{1}{\varepsilon} - 1 \right) $$

given the definition of $\alpha$ in Proposition 3.3. For $p$ fixed, maximizing $\beta(\varepsilon, \gamma)$ over $\varepsilon \in (0, 1)$ and $\gamma > 0$ yields conditions

$$ \frac{1}{p} + \gamma = r\varepsilon^2 \quad \text{and} \quad \left(\frac{1}{p} + \gamma\right)^2 = \frac{r\varepsilon}{p}, $$

which are equivalent to $r\varepsilon^3 = 1/p$ and $\gamma = r\varepsilon^2(1-\varepsilon)$, and the maximum value of $\beta(\varepsilon, \gamma)$ is $r(1-\varepsilon)^3$. Letting $p = r^{-1/4} \in (0, 1)$, and thus $\varepsilon = r^{-1/4} \in (0, 1)$ and $\gamma = r^{1/2} - r^{1/4} > 0$, we obtain the upper bound in (1). Note that the first optimality condition above automatically implies $\alpha = n(\varepsilon^{-1} - 1) > 0$, making the proof consistent. $\square$

Notice that, with this last choice $p = \varepsilon = r^{-1/4}$, $\gamma = (1-\varepsilon)/\varepsilon^2$ and $\alpha = n(1-\varepsilon)/\varepsilon$. In particular, as $r$ tends to infinity, the probability of failure $1-p$ in Proposition 3.3 goes to zero, and $\gamma$ goes to infinity, meaning that the sampling density gets closer to the Christoffel function, as in [CM17].

**Remark 3.5.** In the example where $V_n$ consists of piecewise constant functions on a fixed partition of $X$ into pieces of equal measure, consider $f$ such that $f - P_n f$ is a realisation of a white noise of variance 1. It is easily seen that the optimal approximation error is the Monte-Carlo rate $1+1/r+o(1/r)$, with $\tilde{f}$ made of averages over $r$ random samples in each piece. In that sense, Algorithm 3.1 achieves the optimal decay rate of the error, similar to [CM17, HNP22], but with a remainder independent of $n$.

# 4 Refined randomized sampling algorithm

We now seek to optimize the approximation strategy in the regime $r \approx 1$. Although Algorithm 3.1 works as soon as $r \ge 1$, several improvements can be performed for small values of $r$. First, using (7) instead of (6) reduces the factor $\eta^2/\alpha^2$ to $\eta/\alpha$ in bound (16). Secondly, bounding the weights $s_i$ by a multiple of $1/\gamma$ becomes a crude estimate when $r \to 1$, since $\gamma \to 0$ in that case. Lastly, the use of Markov's inequality in Proposition 3.3 comes at the expense of a factor $1/(1-p)$, which grows as $r$ gets close to 1; to avoid the last issue, we would like the lower barrier $\ell$ to grow in a steady, deterministic fashion. This prevents the sample points $x_i$ from being drawn anywhere in the domain $X$.

For all these reasons, we consider Algorithm 4.1, with inputs $\delta \in (0, 1)$, the size of lower barrier increments, and $\kappa \in [0, 1]$, a parameter balancing the $L^2$ and $L^\infty$ settings. Note that we will always fix $m$ beforehand, and set $\delta = 1/\sqrt{r}$, where $r = m/(n-1)$.

Again, we define $\mathbf{Y}_{m+1} = (\mathbf{A}_m - \ell_m \mathbf{I})^{-1}$, as we would have done at iteration $i = m+1$.
**Algorithm 4.1:** Random sampling with fixed increments of the barrier

**Input:** probability measure $\mu$, parameter $\delta \in (0, 1)$, and $\kappa \in [0, 1]$.

$A_0 = 0 \in \mathbb{R}^{n \times n}$, $\ell_0 = -n$

for $i = 1, \dots, m$ do

Let $\mathbf{Y}_i = (\mathbf{A}_{i-1} - \ell_{i-1} \mathbf{I})^{-1}$

Update $\ell_i = \ell_{i-1} + \delta$

Let $\mathbf{Z}_i = (\mathbf{A}_{i-1} - \ell_i \mathbf{I})^{-1}$

Let $\mathbf{W}_i = (\mathrm{Tr}(\mathbf{Z}_i) - \mathrm{Tr}(\mathbf{Y}_i))^{-1} \mathbf{Z}_i^2 - \mathbf{Z}_i$ and $w_i : x \mapsto \varphi(x)^* \mathbf{W}_i \varphi(x)$

Define $R_i(x) = w_i(x) \mathbf{1}_{w_i(x) \ge \kappa \frac{1-\delta}{\delta}}$

Draw $x_i$ from probability density $\frac{R_i(x)}{\Gamma_i} d\mu(x)$, where $\Gamma_i = \int_X R_i(x) d\mu(x)$

Let $s_i = 1/R_i(x_i) = 1/w_i(x_i)$

Update $\mathbf{A}_i = \mathbf{A}_{i-1} + s_i \varphi(x_i) \varphi(x_i)^*$

end for

**Output:** Sample $\{x_1, \dots, x_m\}$, weights $\{s_1, \dots, s_m\}$, and matrix $\mathbf{A}_m$.

**Lemma 4.1.** The algorithm is well-defined and $\mathrm{Tr}(\mathbf{Y}_i) = 1$ for $1 \le i \le m + 1$.

*Proof.* We use an induction on $i$. First, note that $\mathbf{Y}_1 = \mathbf{I}/n$ hence $\mathrm{Tr}(\mathbf{Y}_1) = 1$. Let $i \in \{1, \dots, m\}$, and assume that $\mathbf{Y}_i$ is positive definite with $\mathrm{Tr}(\mathbf{Y}_i) = 1$. We have

$$ \lambda_{\min}(\mathbf{Y}_i^{-1}) = \lambda_{\max}(\mathbf{Y}_i)^{-1} \ge (\mathrm{Tr}(\mathbf{Y}_i))^{-1} = 1 > \delta, $$

Therefore $\mathbf{Z}_i = (\mathbf{Y}_i^{-1} - \delta \mathbf{I})^{-1}$ is well defined and positive definite, and $\mathbf{W}_i$ is also well defined since $\mathbf{Z}_i \succ \mathbf{Y}_i$. We then need to show that $R_i(x)$ is not zero a.e., in other words that $w_i(x) \ge \kappa(1 - \delta)/\delta$ on a set of positive measure. To this end, we will rely on an averaging argument by showing that $\int_X w_i d\mu$, which is equal to $\mathrm{Tr}(\mathbf{W}_i)$ according to (4), is larger than $(1 - \delta)/\delta$. Recalling (13) and (14), we observe as in Claim 3.6 of [BSS09]

$$ \begin{aligned} \mathrm{Tr}(\mathbf{W}_i) &= \frac{\mathrm{Tr}(\mathbf{Z}_i^2)}{\mathrm{Tr}(\mathbf{Z}_i) - \mathrm{Tr}(\mathbf{Y}_i)} - \mathrm{Tr}(\mathbf{Z}_i) \\ &= \frac{1}{\delta} \left( \frac{\mathrm{Tr}(\mathbf{Z}_i^2)}{\mathrm{Tr}(\mathbf{Y}_i \mathbf{Z}_i)} - \frac{\mathrm{Tr}(\mathbf{Z}_i)}{\mathrm{Tr}(\mathbf{Y}_i)} \right) + \frac{1-\delta}{\delta} \mathrm{Tr}(\mathbf{Z}_i) > \frac{1-\delta}{\delta}, \end{aligned} \quad (17) $$

where we used the induction hypothesis $\mathrm{Tr}(\mathbf{Y}_i) = 1$, and its immediate consequence $\mathrm{Tr}(\mathbf{Z}_i) > 1$. Finally, using the Sherman-Morrison formula as in (9), with $u_i = \sqrt{s_i} \varphi(x_i)$,

$$ \mathrm{Tr}(\mathbf{Y}_{i+1}) = \mathrm{Tr}(\mathbf{Z}_i) - \frac{u_i^* \mathbf{Z}_i^2 u_i}{1 + u_i^* \mathbf{Z}_i u_i} = \mathrm{Tr}(\mathbf{Z}_i) - \frac{\varphi(x_i)^* \mathbf{Z}_i^2 \varphi(x_i)}{w_i(x_i) + \varphi(x_i)^* \mathbf{Z}_i \varphi(x_i)} = \mathrm{Tr}(\mathbf{Y}_i) = 1, $$

which concludes the induction. $\square$

**Remark 4.2.** The symmetric matrix $\mathbf{W}_i$ is not necessarily positive semi-definite, but we have a framing on its maximal eigenvalue

$$ 0 < \frac{\mathrm{Tr}(\mathbf{W}_i)}{n} \le \lambda_{\max}(\mathbf{W}_i) \le \frac{\lambda_{\max}(\mathbf{Z}_i^2)}{\mathrm{Tr}(\mathbf{Z}_i) - \mathrm{Tr}(\mathbf{Y}_i)} \le \frac{\mathrm{Tr}(\mathbf{Z}_i^2)}{\delta \mathrm{Tr}(\mathbf{Y}_i \mathbf{Z}_i)} \le \frac{1}{\delta(1-\delta)}, $$

since $\mathbf{Z}_i \succ \mathbf{Y}_i \succ (1 - \delta) \mathbf{Z}_i$. This property will be useful in the discussion of Section 5.

As $\ell_i = i\delta - n$, we easily get a lower bound on the eigenvalues of $\mathbf{A}_m$.
**Proposition 4.3.** Let $m \ge n$ and define $r = m/(n-1)$. The random matrix $\mathbf{A}_m$ generated by Algorithm 4.1 verifies

$$ \mathbf{A}_m \succsim \alpha \mathbf{I} \quad \text{a.s., with } \alpha = (n-1)(\delta r - 1). $$

*Proof.* After the last iteration of the algorithm, we have

$$ \lambda_{\min}(\mathbf{A}_m) = \ell_m + (\lambda_{\max}(\mathbf{Y}_{m+1}))^{-1} \ge m\delta - n + (\mathrm{Tr}(\mathbf{Y}_{m+1}))^{-1} = (\delta r - 1)(n-1). $$

□

We now have all the tools for the proof of the main theorem.

*Proof of Theorem 1.1, equation (2).* Introducing the set

$$ E_i = \{x \in X, w_i(x) \ge \kappa \frac{1-\delta}{\delta}\}, $$

observe that

$$ \begin{aligned} \Gamma_i &= \int_{E_i} R_i d\mu = \int_X w_i(x) d\mu(x) - \int_{E_i^c} w_i(x) d\mu(x) \\ &\ge \mathrm{Tr}(\mathbf{W}_i) - \int_X \kappa \frac{1-\delta}{\delta} d\mu \ge (1-\kappa) \frac{1-\delta}{\delta}. \end{aligned} \quad (18) $$

Taking the expectation $\mathbb{E}_{x_i}$ over point $x_i$, with $x_1, \dots, x_{i-1}$ fixed, we compute

$$ \mathbb{E}_{x_i}(s_i|g(x_i)|^2) = \int_{E_i} \frac{1}{R_i(x)} |g(x)|^2 \frac{R_i(x)}{\Gamma_i} d\mu(x) \le \frac{\|g\|_{L^2}^2}{\Gamma_i} \le \frac{1}{1-\kappa} \frac{\delta}{1-\delta} \|g\|_{L^2}^2, $$

where $g = f - P_n f \in L^2(X, \mu)$. The same holds for all $1 \le i \le m$, and consequently

$$ \mathbb{E} (\|g\|_m^2) \le \frac{m}{1-\kappa} \frac{\delta}{1-\delta} \|g\|_{L^2}^2, $$

which together with (7) and Proposition 4.3 gives

$$ \mathbb{E}(\|P_n^m g\|^2) \le \alpha^{-1} \mathbb{E}(\|P_n^m g\|_m^2) \le \alpha^{-1} \mathbb{E}(\|g\|_m^2) \le \frac{r}{m(\delta r - 1)} \frac{m}{1-\kappa} \frac{\delta}{1-\delta} \|g\|_{L^2}^2. $$

Taking $\delta = 1/\sqrt{r}$ and adding (5), we arrive at (2).

□

**Remark 4.4.** For $\kappa = 0$, the density of point $x_i$ is proportional to $w_i(x)_+ = \max(w_i(x), 0)$. A sampling density proportional to the positive part of a function also occurs in [LS17], equation (6). There, it is used to select points which increase the lower potential without increasing too much the upper potential.

*Proof of Theorem 1.1, equation (3).* As $V_n \subset L^\infty$, we can assume that $f$ is also in $L^\infty$, otherwise the right-hand side would be infinite. We introduce the notation

$$ P_n^\infty f \in \arg \min_{v \in V_n} \|f - v\|_{L^\infty} $$

for a Chebyshev projection of $f$ onto $V_n$ with respect to the Banach norm $\| \cdot \|_{L^\infty}$, and take $g = f - P_n^\infty f$ the associated residual. As $1/s_i = w_i(x_i) \ge \kappa(1-\delta)/\delta$, it holds

$$ \|g\|_m^2 = \sum_{i=1}^m s_i |g(x_i)|^2 \le \sum_{i=1}^m s_i \|g\|_{L^\infty}^2 \le \frac{m\delta}{\kappa(1-\delta)} \|g\|_{L^\infty}^2 \quad \text{a.s.} $$
By a triangular inequality

$$ \|f - P_n^m f\|_{L^2} \le \|f - P_n^\infty f\|_{L^2} + \|P_n^m f - P_n^\infty f\|_{L^2} = \|g\|_{L^2} + \|P_n^m g\|_{L^2}. $$

We conclude by bounding the operator norm $\|P_n^m\|_{L^\infty \to L^2}$ through

$$ \|P_n^m g\|_{L^2}^2 \le \alpha^{-1} \|P_n^m g\|_m^2 \le \alpha^{-1} \|g\|_m^2 \le \frac{r}{\delta r - 1} \frac{\delta}{\kappa(1-\delta)} \|g\|_{L^\infty}^2. \quad (19) $$

As a result, Algorithm 4.1 with $\delta = 1/\sqrt{r}$ implies inequality (3). □

**Remark 4.5.** In the randomized setting, $\mu$ could be any positive measure on $X$. On the contrary, in the uniform setting, we need $\mu$ to be a probability measure to bound $L^2$ norms by $L^\infty$ norms. Similar results would hold if $\mu$ was a positive measure of finite mass, with all $L^\infty$ norms multiplied by a factor $\mu(X)$, see for instance [LT22, Tem21, KPUU23].

We turn to the proof of the last result of the introduction, which is stated in strong uniform norm to avoid technicalities.

*Proof of Corollary 1.4.* Denote $\mathcal{L}^\infty$ the set of bounded functions, and $\|f\|_{\mathcal{L}^\infty} := \sup_{x \in X} |f(x)|$ the associated norm. As $V_n \subset \mathcal{L}^\infty$, we can assume that $f$ is also in $\mathcal{L}^\infty$, otherwise the right-hand side would be infinite. Taking $g = f - P_n^\infty f$ the $\mathcal{L}^\infty$-residual, we write as above

$$ \|f - P_n^m f\|_{\mathcal{L}^\infty} \le \|f - P_n^\infty f\|_{\mathcal{L}^\infty} + \|P_n^m f - P_n^\infty f\|_{\mathcal{L}^\infty} = \|g\|_{\mathcal{L}^\infty} + \|P_n^m g\|_{\mathcal{L}^\infty}. $$

For any probability measure $\mu$ on $X$, consider $(\varphi_1, \dots, \varphi_n)$ an orthonormal basis of $V_n$ in $L^2 = L^2(X, \mu)$. Then the supremum of the Christoffel function is, by a Cauchy-Schwarz inequality,

$$ K_n := \sup_{x \in X} |\varphi(x)|^2 = \sup_{x \in X} \max_{c \in \mathbb{R}^n \setminus \{0\}} \frac{\left| \sum_{j=1}^n c_j \varphi_j(x) \right|^2}{|c|^2} = \max_{v \in V_n \setminus \{0\}} \frac{\|v\|_{L^\infty}^2}{\|v\|_{L^2}^2}. $$

By simply considering (19) in the case $m = n$, $\kappa = 1$, and $\delta = 1/\sqrt{r} \le 1 - 1/2n$, we derive the inequality $\|P_n^m g\|_{L^2} \le (1-\delta)^{-1} \|g\|_{\mathcal{L}^\infty} \le 2n \|g\|_{\mathcal{L}^\infty}$. As $P_n^m g \in V_n$, we obtain

$$ \|P_n^m g\|_{\mathcal{L}^\infty} \le \sqrt{K_n} \|P_n^m g\|_{L^2} \le 2n \sqrt{K_n} \|g\|_{\mathcal{L}^\infty}. $$

We now rely on the main corollary of the paper [KW60]. Since $X$ is compact, there exists a measure $\mu^*$ for which $K_n = n$. Using such a measure in our framing, we obtain by application of the above $\|P_n^m g\|_{\mathcal{L}^\infty} \le 2n\sqrt{n}\|g\|_{\mathcal{L}^\infty}$, which concludes the proof. □

Let us briefly explain the result in [KW60] to our setting. Given $\phi = (\phi_1, \dots, \phi_n)$ a basis of $V_n$ and $\mu$ a probability measure over $X$, we consider the $n \times n$ symmetric positive semi-definite matrix $M(\mu) := \int_X \phi\phi^* d\mu$. Assuming this matrix is non-singular, $\varphi := M(\mu)^{-1/2}\phi$ is also basis of $V_n$ and is orthonormal with respect to $\mu$. In particular $\varphi(x)^*\varphi(x) = \phi^* M(\mu)^{-1}\phi$ is the associated Christoffel function. The optimal probability measures from [KW60] are characterized by

$$ \mu^* \in \arg\min_{\mu \in \text{Prob}(X)} \sup_{x \in X} \phi(x)^* M(\mu)^{-1} \phi(x) = \arg\max_{\mu \in \text{Prob}(X)} \det(M(\mu)). $$

These two equivalent extremum problems are difficult to solve for general compacts $X$ and spaces $V_n$, thus making the approach from Corollary 1.4 impractical. Nevertheless, it gives a theoretical bound on the approximation error achievable by this method. We refer to [Bos90] for an alternative proof of the result in [KW60], with applications to interpolation, in particular
to the so-called Fejer problem. We also note that [KW60] answers Remark 5.5 in [PU22].

By inspection of the proof of Corollary 1.4, we observe that any measure $\mu$ used in Algorithm 4.1 yields a bound of the form

$$ \|f - P_n^n f\|_{\mathcal{L}^\infty} \le (1 + 2n\sqrt{K_n}) \|f - P_n^\infty f\|_{\mathcal{L}^\infty}. \quad (20) $$

Corollary 1.4 can also be stated for $m \ge n$, as done in Remark 1.5. For instance, Algorithm 4.1 for $m = 2n$, hence with parameters $\delta = 1/\sqrt{r} \le 1/\sqrt{2}$, $\kappa = 1$, and any measure $\mu$, yields

$$ \|f - P_n^m f\|_{\mathcal{L}^\infty} \le (1 + 4\sqrt{K_n}) \|f - P_n^\infty f\|_{\mathcal{L}^\infty}. $$

The above remarks has to be combined with upper bounds on $K_n$, for which we give a few examples in polynomial approximation.

We first consider the space $V_n = \mathbb{P}_{n-1}$ of algebraic polynomials of degree less than $n$, restricted to $X = [-1, 1]$. By considering the uniform measure $d\mu(x) = dx/2$, we have $K_n = \sum_{j=0}^{n-1}(2j + 1) = n^2$, while with the arcsine measure $d\mu(x) = \frac{dx}{\pi\sqrt{1-x^2}}$, we obtain $K_n = 1 + \sum_{j=1}^{n-1} 2 = 2n - 1$, where we used uniform bounds on the Legendre and Chebyshev polynomials, respectively.

Similar results are also available for spaces of multivariate polynomials indexed in lower set, see [CCM+15, Lemma 3.1]. We recall that a lower set, or downward closed set, is a set of multi-indices $\Lambda \subset \mathbb{N}_0^d$ such that $\nu \in \Lambda$ implies $\nu' \in \Lambda$ for all $\nu' \le \nu$, the last inequality being understood component-wise. For $X = [-1, 1]^d$, we consider polynomial approximation spaces of the form

$$ V_n = \operatorname{Span}\{x_1^{\nu_1} \dots x_d^{\nu_d}, \nu \in \Lambda\}, $$

with $\Lambda$ a lower set of cardinality $n$. Among the associated spaces $V_n$, one can find polynomials of bounded degree in each variable, of bounded total degree, as well as hyperbolic cross spaces, which are of important use in high-dimensional approximation [DTU18]. The following estimates hold:

* for $\mu$ the uniform measure, $K_n = \sum_{\nu \in \Lambda} \prod_{j=1}^d (1 + 2\nu_j)$ is bounded by $n^2$,

* for $\mu$ the tensorized arcsine measure, $K_n = \sum_{\nu \in \Lambda} 2^{\nu|0}$ is bounded by $n^{\frac{\log 3}{\log 2}}$.

The estimates $n^2$ and $n^{\frac{\log 3}{\log 2}}$ are established in [CCM+15]. We note also the plain linear estimates $C_\Lambda n$ where $C_\Lambda = \max_{\nu \in \Lambda} \prod_{j=1}^d (1 + 2\nu_j)$ and $C_\Lambda = \max_{\nu \in \Lambda} 2^{\nu|0} \le 2^d$, respectively. Such estimates might be much smaller, for instance in the second setting when $d$ is small. Combining these bounds with (20) implies the results given at the end of the introduction.

In the particular case of polynomials of fixed total degree, invariance by affine transformations allows to treat more general cases.

**Theorem 4.6** ([CD21], Theorems 5.4 and 5.6). For $X \subset \mathbb{R}^d$ a compact domain, $d\mu(x) = dx/|X|$ the uniform measure, and $V_n$ a set of polynomials of fixed total degree, it holds $K_n \le C_X n^2$ if $X$ has a Lipschitz boundary, and $K_n \le C_X n^{1+1/d}$ if $X$ has a smooth boundary.

It is worth noting that the bound $K_n \le C_X n^2$ also holds for arbitrary lower sets when the domain is a union of rectangles of fixed volume, see Theorem 6.5 of [AH20]. All these results show that, in many instances, we have $K_n = O(n^2)$ with the uniform measure. Therefore Algorithm 4.1 provides a feasible way to sample points for which the Lebesgue constant is $L_n = O(n^2)$.
# 5 Numerical aspects

Our two algorithms have a polynomial complexity in $n$, since one has to solve linear systems of size $n$ at each iteration. In fact, $\text{Tr}(\mathbf{Y}_{i+1})$ can be computed without the knowledge of $\mathbf{Y}_{i+1}$, by taking a trace in the Shermann-Morrison formula (9), i.e.

$$ \text{Tr}(\mathbf{Y}_{i+1}) = \text{Tr}(\mathbf{Z}_i) - \frac{\varphi(x_i)^* \mathbf{Z}_i^2 \varphi(x_i)}{s_i + \varphi(x_i)^* \mathbf{Z}_i \varphi(x_i)}, $$

so we only need a matrix inversion for $\mathbf{Z}_i$. Moreover, this can be done efficiently by updating the eigenvalues and eigenvectors of $\mathbf{A}_i$ from one iteration to the next, see [LS18] for details.

In [LS18] and [LS17], polynomial and exponential modifications of the lower potential are also proposed, allowing to draw the samples by batches, and therefore reducing the complexity to a linear expression in $n$ up to logarithmic factors, in the particular case of subsampling a graph laplacian. The proofs are more involved, produce larger constants and require a large value of $r$, so we did not adapt them to our setting.

The essential remaining difficulty is to draw each point $x_i$ according to its prescribed density $\rho_i$ (or $R_i$). We start with the following observation.

**Remark 5.1.** An important application is the case where $X = \{x_1, \dots, x_M\}$ is a finite set, and $\mu$ is the uniform measure on $X$. The orthonormality of the basis $\varphi = (\varphi_1, \dots, \varphi_n)$ is encoded by the following decomposition of the identity into a sum of rank-one matrices:

$$ \sum_{i=1}^{M} v_i v_i^* = \mathbf{I}, \quad v_i := \frac{1}{\sqrt{M}} (\varphi_1(x_i), \dots, \varphi_n(x_i))^\top. $$

In this situation, it suffices to evaluate $\rho_i$ at all points of $X$, and to select $x_i = x$ with probability $\rho_i(x) / \sum_{x' \in X} \rho_i(x')$. Our algorithms 3.1 and 4.1 can also work as a sample reduction technique, see for example [CD22], where the evaluation points are extracted from a large initial sample $x_1, \dots, x_M$. Indeed, it suffices to apply our algorithms to the set $\{x_1, \dots, x_M\}$, equipped with the empirical measure $\mu = \frac{1}{M} \sum_{i=1}^{M} \delta_{x_i}$.

When drawing points directly from an infinite set $X$, one may use an acceptance-rejection strategy. Notice that $\rho_i$ in Algorithm 3.1 and $w_i$ in Algorithm 4.1 are of the form $\varphi(x)^* \mathbf{M} \varphi(x)$, with either $\mathbf{M} = \mathbf{Z}_i + \gamma \mathbf{I}/n$ or $\mathbf{M} = \mathbf{W}_i$. Moreover, some routines have been developed [CM17, AC20, ACD23, Mig21, Mig19, CD21] for drawing points from the Christoffel measure $\frac{1}{n}|\varphi(x)|^2 d\mu(x)$ in relevant multivariate settings.

In the case of Algorithm 3.1, as $\mathbf{M} = \mathbf{Z}_i + \gamma \mathbf{I}/n$ is positive semi-definite, it suffices to draw candidate points $x$ from the Christoffel measure, and to accept them with probability

$$ \frac{\varphi(x)^* \mathbf{M} \varphi(x)}{\lambda_{\max}(\mathbf{M}) |\varphi(x)|^2} \le 1. $$

The probability of accepting a point is

$$ \int_{x \in X} \frac{\varphi(x)^* \mathbf{M} \varphi(x)}{\lambda_{\max}(\mathbf{M}) |\varphi(x)|^2} \frac{|\varphi(x)|^2}{n} d\mu(x) = \frac{\text{Tr}(\mathbf{M})}{n \lambda_{\max}(\mathbf{M})} \ge \frac{1}{n}, $$

therefore we need in average $n\lambda_{\max}(\mathbf{M}) / \text{Tr}(\mathbf{M}) \le n$ draws from the Christoffel measure to find $x_i$.

In the case of Algorithm 4.1, $\mathbf{M} = \mathbf{W}_i$ may have negative eigenvalues, so the bound $\lambda_{\max}(\mathbf{M}) \le \text{Tr}(\mathbf{M})$ may not hold. In order to draw a sample from the probability density
$R_i(x)/\Gamma_i$, we can still rely on acceptance-rejection from the Christoffel measure, and accept points with probability

$$ \frac{R_i(x)}{\lambda_{\max}(\mathbf{W}_i)|\varphi(x)|^2} \in [0, 1]. $$

The acceptance probability is then

$$ \int_{x \in X} \frac{R_i(x)}{\lambda_{\max}(\mathbf{W}_i)|\varphi(x)|^2} \frac{|\varphi(x)|^2}{n} d\mu(x) = \frac{\Gamma_i}{n\lambda_{\max}(\mathbf{W}_i)} \ge \frac{(1-\kappa)(1-\delta)^2}{n} $$

in view of Remark 4.2, and we need to at most $n(1-\kappa)^{-1}(1-\delta)^{-2}$ candidate points in average. When $\kappa = 1$, inspection of (18) and (17) respectively shows that

$$ \Gamma_i \ge \mathrm{Tr}(\mathbf{W}_i) - \frac{1-\delta}{\delta} \quad \text{and} \quad \mathrm{Tr}(\mathbf{W}_i) \ge \frac{1-\delta}{\delta} \mathrm{Tr}(\mathbf{Z}_i), $$

so that

$$ \frac{\delta}{1-\delta} \Gamma_i \ge \mathrm{Tr}(\mathbf{Z}_i) - 1 = \delta \mathrm{Tr}(\mathbf{Y}_i \mathbf{Z}_i) \ge \delta \mathrm{Tr}(\mathbf{Y}_i^2) \ge \frac{\delta}{n} \mathrm{Tr}(\mathbf{Y}_i)^2 = \frac{\delta}{n}, $$

hence we can replace the factor $1-\kappa$ in the acceptance probability by a factor $\delta/n$.

**Remark 5.2.** In the first iterations $1 \le i < n$, the maximal eigenvalue of $M$ has multiplicity at least $n-i$, so the bound on the acceptance probability is $n-i$ times higher. This implies for instance that there are in average less than 2 (respectively, less than $2(1-\kappa)^{-1}(1-\delta)^{-2}$) rejections for $i = 1, \dots, n/2$ in Algorithm 3.1 (respectively, Algorithm 4.1), which we observe in practice.

Even with such methods, the computational complexity remains dominated by the sampling step: Algorithm 3.1 runs in $O(mn^3t)$ time, where the factor $m$ comes from the number of iterations, a factor $n^2$ from the matrix-vector multiplication in the evaluation of $\rho_i$, the last factor $n$ from the average number of rejections, and $t$ is the time needed to generate a point with the Christoffel density. The linear algebra takes $O(n^3)$ time per iteration, hence $O(mn^3)$ time in total. Similarly, Algorithm 4.1 runs in $O(mn^3t(1-\kappa)^{-1}(1-\delta)^{-2})$ time, the only difference coming from the bound on the acceptance ratio. Therefore, applying our sampling strategies may prove quite challenging for large values of $n$.

**Remark 5.3.** In the deterministic setting, it is not necessary to draw each point $x_i$ with a density proportional to $\rho_i(x)$: the only requirement is that $w_i(x_i) \ge \kappa \frac{1-\delta}{\delta}$ for all $1 \le i \le m$. The idea of using the algorithm from [BSS09], still in a deterministic context, but with an infinite set $X$, seems to originate in [DPS+21], which investigates Marcinkiewicz-type discretization theorems. However, it seems that searching for $x_i$ by rejection sampling works better in practice than sorting $X$ and looking for the first point that achieves the above condition [BSU23].

We next comment on the relatively limited importance of the weights. In Algorithm 3.1, the sampling density decomposes as

$$ \frac{\rho_i(x)}{\Xi_i} = (1-p) \frac{\varphi(x)^* \mathbf{Z}_i \varphi(x)}{\mathrm{Tr}(\mathbf{Z}_i)} + p \frac{|\varphi(x)|^2}{n}, $$

with $p = \gamma/\Xi_i$. In others words, we are sampling from a mixture of effective resistance and Christoffel density. Here we only know that $s_i|\varphi(x_i)|^2 \le n\eta/\gamma$, which by the way yields

$$ \lambda_{\max}(\mathbf{A}_m) \le \mathrm{Tr}(\mathbf{A}_m) = \sum_{i=1}^{m} s_i |\varphi(x_i)|^2 \le \frac{mn\eta}{\gamma}. \quad (21) $$
But there is no a priori uniform bound on the weights. If we desire such a bound, it suffices to add a third, constant term in the sampling density:

$$ \frac{\rho_i(x)}{\Xi_i} = \frac{\mathrm{Tr}(\mathbf{Z}_i) \varphi(x)^* \mathbf{Z}_i \varphi(x)}{\Xi_i \mathrm{Tr}(\mathbf{Z}_i)} + \frac{\gamma}{\Xi_i} \frac{|\varphi(x)|^2}{n} + \frac{\gamma_\infty}{\Xi_i} $$

for a new parameter $\gamma_\infty > 0$, where we redefine $\Xi_i$ as $\mathrm{Tr}(\mathbf{Z}_i) + \gamma + \gamma_\infty$. To choose a point according to this density, one can simply draw it from $d\mu(x)$ with probability $\gamma_\infty/\Xi_i$, from the Christoffel measure with probability $\gamma/\Xi_i$, and from the effective resistance measure otherwise. This choice immediately yields $s_i = \eta/\rho_i(x_i) \le \eta/\gamma_\infty$, however taking $\gamma_\infty > 0$ deteriorates a bit the other estimates, which is why we did not include it earlier.

Having a uniform bound on the weights $s_i$ is particularly important in case of noisy samples, where one can only observe $f(x_i) + \varepsilon_i$ for some random $\varepsilon_i$. Then the estimates of Theorem 1.1 hold with an additional noise term of the form

$$ \|\mathbf{L}^+(\sqrt{s_i}\varepsilon_i)_{1\le i\le m}\|_2^2 \le \lambda_{\min}(\mathbf{A}_m)^{-1} \frac{\eta}{\gamma_\infty} \|(\varepsilon_i)_{1\le i\le m}\|_2^2, $$

due the definition of the weighted least-squares estimator in Section 2. We refer to Theorems 5.3 and 5.19 of [ABW22] for details.

Concerning Algorithm 4.1, in the deterministic setting, one can replace each weight $s_i$ by its upper bound $\frac{\delta}{\kappa(1-\delta)}$, resulting in an unweighted discrete norm

$$ \|g\|_m^2 = \frac{1}{m} \sum_{i=1}^m |g(x_i)|^2, $$

without changing our estimates. We refer to [BSU23] for earlier results on subsampling of frames with unweighted discrete norms. In any case, in numerical experiments, the weights seem to have lesser importance than the position of the points.

# 6 Numerical experiments

We conclude with a few illustrations in the space of multivariate polynomials, focusing on a function whose best approximation is explicitly known. Our code is available at the following address: github.com/Belloliva/minimal-oversampling-for-multivariate-polynomial-least-squares.

In the univariate setting, the Legendre polynomials $(L_k)_{k\ge 0}$ are orthonormal in $L^2([-1, 1], dx/2)$, and satisfy

$$ \|L_k\|_{L^\infty} = L_k(1) = \sqrt{2k+1}. $$

We recall that their generating function is given by

$$ g_y(x) = \sum_{k=0}^{\infty} \frac{L_k(x)}{L_k(1)} y^k = \frac{1}{\sqrt{1-2xy+y^2}}, \quad x \in [-1, 1], \quad y \in (-1, 1). $$

Now, on the domain $X = [-1, 1]^d$ with $d \in \mathbb{N}$, one can consider tensorized Legendre polynomials

$$ L_k(\boldsymbol{x}) := \prod_{j=1}^{d} L_{k_j}(x_j), \quad \boldsymbol{x} = (x_1, \dots, x_d) \in X, \quad \boldsymbol{k} = (k_1, \dots, k_d) \in \mathbb{N}_0^d, $$
which are orthonormal for the uniform measure $d\mu(\boldsymbol{x}) = d\boldsymbol{x}/2^d$ over $X$.

Fixing $\boldsymbol{y} = (y_1, \dots, y_d) \in (0, 1)^d$, we would like to approximate the multivariate function

$$ \boldsymbol{x} \mapsto g_{\boldsymbol{y}}(\boldsymbol{x}) = \prod_{j=1}^{d} g_{y_j}(x_j) = \sum_{\boldsymbol{k} \in \mathbb{N}_0^d} c_{\boldsymbol{k}} L_{\boldsymbol{k}}(\boldsymbol{x}), \quad c_{\boldsymbol{k}} = \prod_{j=1}^{d} \frac{y_j^{k_j}}{\sqrt{2k_j + 1}}. $$

This is a typical representer of holomorphic functions on $X$ with anisotropic dependence on $\boldsymbol{x}$ parameterized by $\boldsymbol{y}$. Its $L^2(X, \mu)$ norm is analytically given by

$$ \|g_{\boldsymbol{y}}\|_{L^2}^2 = \prod_{j=1}^{d} \int_{-1}^{1} |g_{y_j}(x_j)|^2 \frac{dx_j}{2} = \prod_{j=1}^{d} \frac{1}{2y_j} \ln \left( \frac{1+y_j}{1-y_j} \right). $$

Moreover, the best approximation space of dimension $n$ for $g_{\boldsymbol{y}}$ w.r.t. $\|\cdot\|_{L^2}$ is

$$ V_n = \text{Span}\{L_k, k \in \Lambda_n\}, $$

where $\Lambda_n$ is the set of $n$ multi-indices $k$ with largest values of $c_k$. Note that $\Lambda_n$ is lower and is easily computed, since the coefficients $c_k$ are decreasing for the component-wise order, and that the approximation error is

$$ \min_{v \in V_n} \|g_y - v\|_{L^2}^2 = \left\| g_y - \sum_{k \in \Lambda_n} c_k L_k \right\|_{L^2}^2 = \|g_y\|_{L^2}^2 - \sum_{k \in \Lambda_n} c_k^2. \quad (22) $$

**Remark 6.1.** When looking at approximations in the uniform norm, although $V_n$ and $\sum_{k \in \Lambda_n} c_k L_k$ may not be optimal, one can still compute an upper bound

$$ \min_{v \in V_n} \|g_y - v\|_{L^\infty} \le \left\| g_y - \sum_{k \in \Lambda_n} c_k L_k \right\|_{L^\infty} = \prod_{j=1}^{d} \frac{1}{1 - y_j} - \sum_{k \in \Lambda_n} \prod_{j=1}^{d} y_j^{k_j}, $$

since the maximum is attained at $(1, \dots, 1) \in X$.

In Figure 1, we take parameters $\boldsymbol{y} = (0.9, 0.8, 0.7, 0.6)$ in spatial dimension $d=4$, use Legendre expansions of size $n=128$, and draw $m=2n=256$ points $\boldsymbol{x}_1, \dots, \boldsymbol{x}_m \in X$ according to one of the following strategies:

(a) The points are i.i.d according to the uniform measure $d\mu(\boldsymbol{x})$

(b) The points are i.i.d according to the tensor product arcsine measure

$$ \prod_{j=1}^{d} \frac{dx_j}{\pi \sqrt{1 - x_j^2}} $$

(c) The points are i.i.d according to the Christoffel measure

$$ \frac{1}{n} \sum_{k \in \Lambda_n} |L_k(\boldsymbol{y})|^2 d\mu(\boldsymbol{x}) $$

(d) The points and weights are generated by Algorithm 3.1 with input parameters $\varepsilon = r^{-1/4}$ and $\gamma = r^{1/2} - r^{1/4}$, where $r = (m+1)/n$

(e) The points and weights are generated by Algorithm 4.1 with input parameters $\delta = r^{-1/2}$ and $\kappa = 1/2$, where $r = m/(n-1)$.
<figure><img src="image_4.png" /><figcaption>Figure 1: Histograms of condition number of $\mathbf{A}_m$ for various random sampling strategies</figcaption></figure>

In the first three cases, the weights $s_i$ are taken as the inverse of the sampling density (w.r.t. $d\mu(\mathbf{x})$) at points $\mathbf{x}_i$. The uniform (a) and Christoffel (c) settings are the ones studied in [CDL13] and [CM17], and the arcsine measure (b) is the limit of the Christoffel measure when $d$ is fixed and $n$ tends to infinity, thus yielding similar properties while being slightly simpler to sample from.

Each scheme outputs a random matrix $\mathbf{A}_m = (\sum_{i=1}^m s_i L_k(\mathbf{x}_i) L_{k'}(\mathbf{x}_i))_{\mathbf{k},\mathbf{k}' \in \Lambda_n}$, and Figure 1 displays histograms of their condition number $\lambda_{\max}(\mathbf{A}_m)/\lambda_{\min}(\mathbf{A}_m)$ over 400 runs. Indeed, this condition number is a classical indicator of the accuracy and robustness of least-squares, compared to the best possible approximation.

We immediately observe that the last two methods achieve smaller condition numbers, often between 30 and 40 for Algorithm 3.1, and between 20 and 30 for Algorithm 4.1, compared to i.i.d points, which generally give values larger than 40, and sometimes much larger, especially in the case of uniform points. Heuristically, we expect the condition number of Algorithm 4.1 to be in average a bit larger than the eigenvalues of

$$ \frac{\mathbb{E}(\mathbf{A}_m)}{\lambda_{\min}(\mathbf{A}_m)} \lesssim \frac{(1-\kappa)^{-1}m\delta(1-\delta)^{-1}}{(n-1)(\delta r - 1)} \mathbf{I} = \frac{2}{(1-1/\sqrt{r})^2} \mathbf{I} \lesssim 22.9 \mathbf{I}, $$

where we used arguments from the proof of Theorem 1.1, equation (2), together with the fact that $r = \frac{256}{127}$. Therefore, there is a very good agreement between the theoretical and numerical spectral properties of $\mathbf{A}_m$.

It should be mentioned that in Algorithm 3.1, the random event from Proposition 3.3 was always realized, resulting in usual least-squares instead of the conditioned version (15). Nevertheless, redrawing the whole sample when the final condition number is large, as proposed in [HNP22] for Christoffel points, remains a good option in practice for all sampling schemes, since it amounts to truncating the tails of the above histograms.

In Figure 2, we plot the number of rejections as a function of the iteration index $i = 1, \dots, m$, averaged over the 400 runs, for Algorithms 3.1 and 4.1. In accordance with Remark 5.2, few rejections are observed in the first iterations $i < n$. However, the number of rejections remains moderate in the second half $n \le i \le m$, staying below 15 in the case of Algorithm 3.1, which is much better than the pessimistic upper bound $n = 128$ from Section 5. Algorithm 4.1 incurs about 7 times more rejections, which is again smaller than the factor $(1 - \kappa)^{-1}(1 - \delta)^{-2} \approx 23$ encountered in the analysis.

Finally, in Figure 3, we keep the same parameters $d = 4$, $\mathbf{y} = (0.9, 0.8, 0.7, 0.6)$ and $n = 128$, and plot the normalized error

$$ \mathbb{E}_{\text{emp}}(\|g_{\mathbf{y}} - P_n^m g_{\mathbf{y}}\|_{L^2})/E^* $$

as a function of the number of samples $m$, for $m = n + 2p$ and $p = 0, \dots, 40$. Here $\mathbb{E}_{\text{emp}}$ stands for an empirical average over 100 runs, and $E^* \approx 0.402882$ is the error of best approximation, computed by (22). For the sake of clarity, we only display the results for the Christoffel measure
<figure><img src="image_5.png" /><figcaption>Figure 2: Average over 400 runs of the number of rejections in each iteration</figcaption></figure>

and our algorithms. We again observe that the latter perform better, especially in the regime $m \approx n$, and that they are within a factor 2 of the optimal error as soon as $m \ge n + 30$.

<figure><img src="image_6.png" /><figcaption>Figure 3: Ratio between the least-squares error and the best approximation error for $n = 128$ and $m$ ranging between 128 and 208, when the points and weights are i.i.d according to the Christoffel measure, or generated by Algorithms 3.1 and 4.1.</figcaption></figure>

**Remark 6.2.** In view of the discussion from Section 5, the last three schemes necessitate sampling from Christoffel measure, which can be efficiently implemented by first drawing $\mathbf{k}$ uniformly in $\Lambda_n$, and then drawing $\mathbf{x}$ from the measure

$$|L_{\mathbf{k}}(\mathbf{x})|^2 d\mu(\mathbf{y}) = \prod_{j=1}^{d} |L_{k_j}(x_j)|^2 \frac{dx_j}{2}.$$

For the second step, exploiting the product structure of $L_{\mathbf{k}}$, it suffices to draw each component $x_j$ independently from $|L_{k_j}(x)|^2 \frac{dx_j}{2}$. When $k_j = 0$, this is just the uniform measure over $[-1, 1]$. Otherwise, we rely once more on acceptance/rejection from the arcsine measure, that is from the cosine of uniform points. Thanks to the so-called Berstein inequality for Legendre polynomials

$$\frac{|L_k(x)|^2}{2} \le \frac{2k+1}{k} \frac{1}{\pi\sqrt{1-x^2}},$$

the acceptance probability is at least $k_j/(2k_j + 1) \ge 1/3$ for $k_j \ge 1$, so there are at most $2 \min(k, d)$ rejections in average. As a conclusion, the time needed to generate a
point from the Christoffel probability measure is $t = O(d \min(n, d))$. We refer to Section 5 of [CM17], as well as [CD21] and [Mig21], for an overview of efficient sampling strategies on more general domains.

**Aknowledgement:** The authors would like to thank Albert Cohen for insightful feedback and discussions all along the elaboration of the paper, David Krieg, Mario Ullrich and Tino Ullrich for their enriching questions and comments, and the reviewers for their careful reading and valuable suggestions.

## References

[ABW22] Ben Adcock, Simone Brugiapaglia, and Clayton G. Webster. *Sparse polynomial approximation of high-dimensional functions*, volume 25 of *Computational Science & Engineering*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, [2022] ©2022.

[AC20] Ben Adcock and Juan M. Cardenas. Near-optimal sampling strategies for multivariate function approximation on general domains. *SIAM J. Math. Data Sci.*, 2(3):607–630, 2020.

[ACD23] Ben Adcock, Juan M. Cardenas, and Nick Dexter. An adaptive sampling and domain learning strategy for multivariate function approximation on unknown domains. *SIAM J. Sci. Comput.*, 45(1):A200–A225, 2023.

[AH20] Ben Adcock and Daan Huybrechs. Approximating smooth, multivariate functions on irregular domains. *Forum of Mathematics, Sigma*, 8:e26, 2020.

[AN22] Vladimir Andrievskii and Fedor Nazarov. A simple upper bound for Lebesgue constants associated with Leja points on the real line. *J. Approx. Theory*, 275:Paper No. 105699, 13, 2022.

[APS19] Ben Adcock, Rodrigo B. Platte, and Alexei Shadrin. Optimal sampling rates for approximating analytic functions from pointwise samples. *IMA J. Numer. Anal.*, 39(3):1360–1390, 2019.

[Bos90] Len Bos. Some remarks on the Fejér problem for Lagrange interpolation in several variables. *J. Approx. Theory*, 60(2):133–140, 1990.

[BSS09] Joshua D Batson, Daniel A Spielman, and Nikhil Srivastava. Twice-ramanujan sparsifiers. In *Proceedings of the forty-first annual ACM symposium on Theory of computing*, pages 255–262, 2009.

[BSU23] Felix Bartel, Martin Schäfer, and Tino Ullrich. Constructive subsampling of finite frames with applications in optimal function recovery. *Appl. Comput. Harmon. Anal.*, 65:209–248, 2023.

[BX09] John P. Boyd and Fei Xu. Divergence (Runge phenomenon) for least-squares polynomial approximation on an equispaced grid and Mock-Chebyshev subset interpolation. *Appl. Math. Comput.*, 210(1):158–168, 2009.

[CC15] Abdellah Chkifa and Albert Cohen. On the stability of polynomial interpolation using hierarchical sampling. In *Sampling theory, a renaissance*, Appl. Numer. Harmon. Anal., pages 437–458. Birkhäuser/Springer, Cham, 2015.
[CCM+15] Abdellah Chkifa, Albert Cohen, Giovanni Migliorati, Fabio Nobile, and Raul Tem- pone. Discrete least squares polynomial approximation with random evaluations— application to parametric and stochastic elliptic PDEs. *ESAIM Math. Model. Nu- mer. Anal.*, 49(3):815–837, 2015.

[CCS14] Abdellah Chkifa, Albert Cohen, and Christoph Schwab. High-dimensional adaptive sparse polynomial interpolation and applications to parametric PDEs. *Found. Comput. Math.*, 14(4):601–633, 2014.

[CD21] Albert Cohen and Matthieu Dolbeault. Optimal sampling and christoffel functions on general domains. *Constructive Approximation*, pages 1–43, 2021.

[CD22] Albert Cohen and Matthieu Dolbeault. Optimal pointwise sampling for $L^2$ approximation. *Journal of Complexity*, 68:101602, 2022.

[CDL13] Albert Cohen, Mark A. Davenport, and Dany Leviatan. On the stability and accuracy of least squares approximations. *Found. Comput. Math.*, 13(5):819–834, 2013.

[CM17] Albert Cohen and Giovanni Migliorati. Optimal weighted least-squares methods. *The SMAI journal of computational mathematics*, 3:181–203, 2017.

[CP19] Xue Chen and Eric Price. Active regression via linear-sample sparsification. In *Conference on Learning Theory*, pages 663–695. PMLR, 2019.

[DKU23] Matthieu Dolbeault, David Krieg, and Mario Ullrich. A sharp upper bound for sampling numbers in $L^2$. *Applied and Computational Harmonic Analysis*, 63:113–134, 2023.

[DPS+21] Feng Dai, Andriy Prymak, Alexei Shadrin, Vladimir Temlyakov, and Serguey Tikhonov. Entropy numbers and Marcinkiewicz-type discretization. *J. Funct. Anal.*, 281(6):Paper No. 109090, 25, 2021.

[DT24] Feng Dai and Vladimir Temlyakov. Random points are good for universal discretization. *Journal of Mathematical Analysis and Applications*, 529(1):127570, 2024.

[DTU18] Dinh Dung, Vladimir Temlyakov, and Tino Ullrich. *Hyperbolic cross approximation*. Advanced Courses in Mathematics. CRM Barcelona. Birkhäuser/Springer, Cham, 2018. Edited and with a foreword by Sergey Tikhonov.

[FS19] Daniel Freeman and Darrin Speegle. The discretization problem for continuous frames. *Adv. Math.*, 345:784–813, 2019.

[GW24] Jiaxin Geng and Heping Wang. On the power of standard information for tractability for $L_\infty$ approximation of periodic functions in the worst case setting. *Journal of Complexity*, 80:101790, 2024.

[HNP22] Cécile Haberstich, Anthony Nouy, and Guillaume Perrin. Boosted optimal weighted least-squares. *Math. Comp.*, 91(335):1281–1315, 2022.

[KPUU23] David Krieg, Kateryna Pozharska, Mario Ullrich, and Tino Ullrich. Sampling recovery in the uniform norm. *arXiv preprint arXiv:2305.07539*, 2023.

[KU21a] David Krieg and Mario Ullrich. Function values are enough for $L_2$-approximation. *Found. Comput. Math.*, 21(4):1141–1151, 2021.
[KU21b] David Krieg and Mario Ullrich. Function values are enough for $L_2$-approximation: Part II. *J. Complexity*, 66:Paper No. 101569, 14, 2021.

[KUV21] Lutz Kämmerer, Tino Ullrich, and Toni Volkmer. Worst-case recovery guarantees for least squares approximation using random samples. *Constr. Approx.*, 54(2):295–352, 2021.

[KW60] Jack Kiefer and Jacob Wolfowitz. The equivalence of two extremum problems. *Canadian J. Math.*, 12:363–366, 1960.

[LS17] Yin Tat Lee and He Sun. An SDP-based algorithm for linear-sized spectral sparsification. In *STOC'17—Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing*, pages 678–687. ACM, New York, 2017.

[LS18] Yin Tat Lee and He Sun. Constructing linear-sized spectral sparsification in almost-linear time. *SIAM J. Comput.*, 47(6):2315–2336, 2018.

[LT22] Irina Limonova and Vladimir Temlyakov. On sampling discretization in $L_2$. *J. Math. Anal. Appl.*, 515(2):Paper No. 126457, 14, 2022.

[Mig19] Giovanni Migliorati. Adaptive approximation by optimal weighted least-squares methods. *SIAM J. Numer. Anal.*, 57(5):2217–2245, 2019.

[Mig21] Giovanni Migliorati. Multivariate approximation of functions on irregular domains by weighted least-squares methods. *IMA J. Numer. Anal.*, 41(2):1293–1317, 2021.

[MM08] Giuseppe Mastroianni and Gradimir V. Milovanovic. *Interpolation Processes: Basic Theory and Applications*. Springer Publishing Company, Incorporated, 1 edition, 2008.

[MSS15] Adam W. Marcus, Daniel A. Spielman, and Nikhil Srivastava. Interlacing families II: Mixed characteristic polynomials and the Kadison—Singer problem. *Annals of Mathematics*, pages 327–350, 2015.

[MU21] Moritz Moeller and Tino Ullrich. $L_2$-norm sampling discretization and recovery of functions from RKHS with finite trace. *Sampl. Theory Signal Process. Data Anal.*, 19(2):Paper No. 13, 31, 2021.

[NOU13] Shahaf Nitzan, Alexander Olevskii, and Alexander Ulanovskii. A few remarks on sampling of signals with small spectrum. *Tr. Mat. Inst. Steklova*, 280(Ortogonal'nye Ryady, Teoriya Priblizheniui i Smezhnye Voprosy):247–254, 2013.

[NOU16] Shahaf Nitzan, Alexander Olevskii, and Alexander Ulanovskii. Exponential frames on unbounded sets. *Proc. Amer. Math. Soc.*, 144(1):109–118, 2016.

[Nov88] Erich Novak. *Deterministic and stochastic error bounds in numerical analysis*, volume 1349 of *Lecture Notes in Mathematics*. Springer-Verlag, Berlin, 1988.

[NSU22] Nicolas Nagel, Martin Schäfer, and Tino Ullrich. A new upper bound for sampling numbers. *Foundations of Computational Mathematics*, 22(2):445–468, 2022.

[PTK11] Rodrigo B. Platte, Lloyd N. Trefethen, and Arno B. J. Kuijlaars. Impossibility of fast stable approximation of analytic functions from equispaced samples. *SIAM Rev.*, 53(2):308–318, 2011.
[PU22] Kateryna Pozharska and Tino Ullrich. A note on sampling recovery of multivariate functions in the uniform norm. *SIAM J. Numer. Anal.*, 60(3):1363–1384, 2022.

[SS11] Daniel A. Spielman and Nikhil Srivastava. Graph sparsification by effective resistances. *SIAM J. Comput.*, 40(6):1913–1926, 2011.

[Tem21] Vladimir Temlyakov. On optimal recovery in $L_2$. *J. Complexity*, 65:Paper No. 101545, 11, 2021.