{'fraction_non_alphanumeric': 0.04927057132590729, 'fraction_numerical': 0.013855551563061444, 'mean_word_length': 3.844450633616488, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multivariate approximation in downward closed polynomial spaces

The task of approximating a function of d variables from its evaluations at a given number of points is ubiquitous in numerical analysis and engineering applications. When d is large, this task is challenged by the so-called curse of dimensionality. As a typical example, standard polynomial spaces, such as those of total degree type, are often uneffective to reach a prescribed accuracy unless a prohibitive number of evaluations is invested. In recent years it has been shown that, for certain relevant applications, there are substantial advantages in using certain sparse polynomial spaces having anisotropic features with respect to the different variables. These applications include in particular the numerical approximation of high-dimensional parametric and stochastic partial differential equations. We start by surveying several results in this direction, with an emphasis on the numerical algorithms that are available for the construction of the approximation, in particular through interpolation or discrete least-squares fitting. All such algorithms rely on the assumption that the set of multi-indices associated with the polynomial space is downward closed. In the present paper we introduce some tools for the study of approximation in multivariate spaces under this assumption, and use them in the derivation of error bounds, sometimes independent of the dimension d, and in the development of adaptive strategies.


Introduction
The mathematical modeling of complex physical phenomena often demands for functions that depend on a large number of variables. One typical instance occurs when a quantity of interest u is given as the solution to an equation written in general form as P(u, y) = 0, (1.1) where y = (y j ) j=1,...,d ∈ R d is a vector that concatenates various physical parameters which have an influence on u.
Supposing that we are able to solve the above problem, either exactly or approximately by numerical methods, for any y in some domain of interest U ⊂ R d we thus have access to the parameter-to-solution map y → u(y). (1. 2) The quantity u(y) may be of various forms, namely: 1. a real number, that is, u(y) ∈ R; 2. a function in some Banach space, for example when (1.1) is a partial differential equation (PDE); 3. a vector of eventually large dimension, in particular when (1.1) is a PDE whose solution is numerically approximated using some numerical method with fixed discretization parameters.
As a guiding example which will be further discussed in this paper, consider the elliptic diffusion equation Here, a and ψ j are given functions in L ∞ (D), and the y j range in finite intervals that, up to renormalization, can all be assumed to be [−1, 1]. In this example y = (y j ) j≥1 is countably infinite-dimensional, that is, d = ∞. The standard weak formulation of (1.3) in H 1 0 (D), is ensured to be well-posed for all such a under the so-called uniform ellipticity assumption j≥1 |ψ j | ≤ a − r, a.e. on D, (1.5) for some r > 0. In this case, the map y → u(y) acts from U = [−1, 1] N to H 1 0 (D). However, if we consider the discretization of (1.3) in some finite element space V h ⊂ H 1 0 (D), where h refers to the corresponding mesh size, using for instance the Galerkin method, then the resulting map y → u h (y), acts from U = [−1, 1] N to V h . Likewise, if we consider a quantity of interest such as the flux q(u) = Σ a∇u · σ over a given interface Σ ⊂ D with σ being the outward pointing normal vector, then the resulting map y → q(y) = q(u(y)), acts from U = [−1, 1] N to R. In all three cases, the above maps act from U to some finite-or infinite-dimensional Banach space V , which is either H 1 0 , V h or R. In the previous instances, the functional dependence between the input parameters y and the output u(y) is described in clear mathematical terms by equation (1.1). In other practical instances, the output u(y) can be the outcome of a complex physical experiment or numerical simulation with input parameter y. However the dependence on y might not be given in such clear mathematical terms.
In all the abovementioned cases, we assume that are we are able to query the map (1.2) at any given parameter value y ∈ U , eventually up to some uncertainty. Such uncertainty may arise due to: (i) measurement errors, when y → u(y) is obtained by a physical experiment, or (ii) computational errors, when y → u(y) is obtained by a numerical computation.
The second type of errors may result from the spatial discretization when solving a PDE with a given method, and from the round-off errors when solving the associated discrete systems.
One common way of modeling such errors is by assuming that we observe u at the selected points y up to a an additive noise η which may depend on y, that is, we evaluate y → u(y) + η(y), (1.6) where η satisfies a uniform bound for some ε > 0 representing the noise level. Queries of the exact u(y) or of the noisy u(y) + η(y) are often expensive since they require numerically solving a PDE, or setting up a physical experiment, or running a time-consuming simulation algorithm. A natural objective is therefore to approximate the map (1.2) from some fixed number m of such queries at points {y 1 , . . . , y m } ∈ U . Such approximations y → u(y) are sometimes called surrogate or reduced models.
Let us note that approximation of the map (1.2) is sometimes a preliminary task for solving other eventually more complicated problems, such as: 1. Optimization and Control, i.e. find a y which minimizes a certain criterion depending on u(y). In many situations, the criterion takes the form of a convex functional of u(y), and the minimization is subject to feasibility constraints. See e.g. the monographs [3,30] and references therein for an overview of classical formulations and numerical methods for optimization problems.
2. Inverse Problems, i.e. find an estimate y from some data depending on the output u(y). Typically, we face an ill-posed problem, where the parameter-tosolution map does not admit a global and stable inverse. Nonetheless, developing efficient numerical methods for approximating the parameter-to-solution map, i.e. solving the so-called direct problem, is a first step towards the construction of numerical methods for solving the more complex inverse problem, see e.g. [33]. 3. Uncertainty Quantification, i.e. describe the stochastic properties of the solution u(y) in the case where the parameter y is modeled by a random variable distributed according to a given probability density. We may for instance be interested in computing the expectation or variance of the V -valued random variable u(y). Note that this task amounts in computing multivariate integrals over the domain U with respect to the given probability measure. This area also embraces, among others, optimization and inverse problems whenever affected by uncertainty in the data. We refer to e.g. [24] for the application of polynomial approximation to uncertainty quantification, and to [32] for the Bayesian approach to inverse problems.
There exist many approaches for approximating an unknown function of one or several variables from its evaluations at given points. One of the most classical approach consists in picking the approximant in a given suitable n-dimensional space of elementary functions, such that n ≤ m. (1.8) Here, by "suitable" we mean that the space should have the ability to approximate the target function to some prescribed accuracy, taking for instance advantage of its smoothness properties. By "elementary" we mean that such functions should have simple explicit form which can be efficiently exploited in numerical computations. The simplest type of such functions are obviously polynomials in the variables y j . As a classical example, we may decide to use, for some given k ∈ N 0 , the total degree polynomial space of order k, namely Note that since u(y) is V -valued, this means that we actually use the V -valued polynomial space Another classical example is the polynomial space of degree k in each variable, namely A critical issue encountered by choosing such spaces is the fact that, for a fixed value of k, the dimension of P k grows with d like d k , and that of Q k grows like k d , that is, exponentially in d. Since capturing the fine structure of the map (1.2) typically requires a large polynomial degree k in some coordinates, we expect in view of (1.8) that the number of needed evaluations m becomes prohibitive as the number of variables becomes large. This state of affairs is a manifestation of the so-called curse of dimensionality. From an approximation theoretic or information-based complexity point of view, the curse of dimensionality is expressed by the fact that functions in standard smoothness classes such as C s (U ) cannot be approximated in L ∞ (U ) with better rate then n −s/d by any method using n degrees of freedom or n evaluations, see e.g. [18,31]. Therefore, in high dimension, one is enforced to give up on classical polynomial spaces of the above form, and instead consider more general spaces of the general form P Λ := span{y → y ν : ∈ Λ}, (1.9) where Λ is a subset of N d 0 with a given cardinality n := #(Λ). In the case of infinitely many variables d = ∞, we replace N d 0 by the set of finitely supported sequences of nonnegative integers. For V -valued functions, we thus use the space Note that V Λ = P Λ in the particular case where V = R.
The main objective when approximating the map (1.2) is to maintain a reasonable trade-off between accuracy measured in a given error norm and complexity measured by n, exploiting the different importance of each variable. Intuitively, large polynomial degrees should only be allocated to the most important variables. In this sense, if d is the dimension and k is the largest polynomial degree in any variable appearing in Λ, we view Λ as a very sparse subset of {0, . . . , k} d .
As generally defined by (1.9), the space P Λ does not satisfy some natural properties of usual polynomial spaces such as closure under differentiation in any variable, or invariance by a change of basis when replacing the monomials y ν by other tensorized basis functions of the form where the univariate functions {φ 0 , . . . , φ k } form a basis of P k for any k ≥ 0, for example with the Legendre or Chebyshev polynomials. In order to fulfill these requirements, we ask that the set Λ has the following natural property.
where ν ≤ ν means that ν j ≤ ν j for all j.
Downward closed sets are also called lower sets. We sometimes use the terminology of downward closed polynomial spaces for the corresponding P Λ . To our knowledge, such spaces have been first considered in [23] in the bivariate case d = 2 and referred to as polynômes pleins. Their study in general dimension d has been pursued in [25] and [16]. The objective of the present paper is to give a survey of recent advances on the use of downward closed polynomial spaces for high-dimensional approximation.
The outline is the following. We review in Section 2 several polynomial approximation results obtained in [1,2] in which the use of well-chosen index sets allows one to break the curse of dimensionality for relevant classes of functions defined on U = [−1, 1] d , e.g. such as those occuring when solving the elliptic PDE (1.3) with parametric diffusion coefficient (1.4). Indeed, we obtain an algebraic convergence rate n −s , where s is independent of d in the sense that such a rate may even hold when d = ∞. Here, we consider the error between the map (1.2) and its approximant in We also consider the case of lognormal diffusion coefficients of the form where the y j are i.i.d. standard Gaussians random variables. In this case, we have U = R d and the error is measured in L 2 (U, V, dγ) where is the tensorized Gaussian probability measure. The above approximation results are established by using n-term truncations of polynomial expansions, such as Taylor, Legendre or Hermite, which do not necessarily result in downward closed index sets.
In the present paper we provide a general approach to establish similar convergence rates with downward closed polynomial spaces. The coefficients in the polynomial expansions cannot be computed exactly from a finite number of point evaluations of (1.2). One first numerical procedure that builds a polynomial approximation from point evaluations is interpolation. In this case the number m of samples is exactly equal to the dimension n of the polynomial space. We discuss in Section 3 a general strategy to choose evaluation points and compute the interpolant in arbitrarily high dimension. One of its useful features is that the evaluations and interpolants are updated in a sequential manner as the polynomial space is enriched. We study the stability of this process and its ability to achieve the same convergence rates in L ∞ established in Section 2.
A second numerical procedure for building a polynomial approximation is the least-squares method, which applies to the overdetermined case m > n. To keep the presentation concise, we confine to results obtained in the analysis of this method only for the case of evaluations at random points. In Section 4 we discuss standard least squares, both in the noisy and noiseless cases, and in particular explain under which circumstances the method is stable and compares favorably with the best approximation error in L 2 . Afterwards we discuss the more general method of weighted least squares, which allows one to optimize the relation between the dimension of the polynomial space n and the number of evaluation points m that warrants stability and optimal accuracy.
The success of interpolation and least squares is critically tied to the choice of proper downward closed sets (Λ n ) n≥1 with #(Λ n ) = n. The set Λ * n that gives the best polynomial approximation among all possible downward closed sets of cardinality n is often not accessible. In practice we need to rely on some a-priori analysis to select "suboptimal yet good" sets. An alternative strategy is to select the sequence (Λ n ) n≥1 in an adaptive manner, that is, make use of the computation of the approximation for Λ n−1 in order to choose Λ n . We discuss in Section 5 several adaptive and nonadaptive strategies which make use of the downward closed structure of such sets. While our paper is presented in the framework of polynomial approximation, the concept of downward closed set may serve to define multivariate approximation procedures in other nonpolynomial frameworks. At the end of the paper we give some remarks on this possible extension, including, as a particular example, approximation by sparse piecewise polynomial spaces using hierarchical bases, such as sparse grid spaces.
Let us finally mention that another class of frequently used methods in highdimensional approximation is based on Reproducing Kernel Hilbert Space (RKHS) or equivalently Gaussian process regression, also known as kriging. In such methods, for a given Mercer kernel K(·, ·) the approximant is typically searched by minimizing the associated RKHS norm among all functions agreeing with the data at the evaluation points, or equivalently by computing the expectation of a Gaussian process with covariance function K conditional to the observed data. Albeit natural competitors, these methods do not fall in the category discussed in the present paper, in the sense that the space where the approximation is picked varies with the evaluation points. It is not clear under which circumstances they may also break the curse of dimensionality. where for each ν = (ν j ) j≥1 ∈ F the function φ ν : U → R has the tensor product form and u ν ∈ V . Here we assume that (φ k ) k≥0 is a sequence of univariate polynomials such that φ 0 ≡ 1 and the degree of φ k is equal to k. This implies that {φ 0 , . . . , φ k } is a basis of P k and that the above product only involves a finite number of factors, even in the case where d = ∞. Thus, we obtain polynomial approximations of (1.2) by fixing some sets Λ n ⊂ F with #(Λ n ) = n and defining Before discussing specific examples, let us make some general remarks on the truncation of countable expansions with V -valued coefficients, not necessarily of tensor product or polynomial type.

Definition 2.1
The series (2.1) is said to converge conditionally with limit u in a given norm · if and only if there exists an exhaustion (Λ n ) n≥1 of F (which means that for any ν ∈ F there exists n 0 such that ν ∈ Λ n for all n ≥ n 0 ), with the convergence property lim The series (2.1) is said to converge unconditionally towards u in the same norm, if and only if (2.3) holds for every exhaustion (Λ n ) n≥1 of F .
As already mentioned in the introduction, we confine our attention to the error norms L ∞ (U, V ) or L 2 (U, V ) with respect to the uniform probability measure dµ. We are interested in establishing unconditional convergence, as well as estimates of the error between u and its truncated expansion, for both types of norm.
In the case of the L 2 norm, unconditional convergence can be established when (φ ν ) ν∈F is an orthonormal basis of L 2 (U ). In this case we know from standard Hilbert space theory that if (1.2) belongs to L 2 (U, V ) then the inner products are elements of V , and the series (2.1) converges unconditionally towards u in L 2 (U, V ). In addition, the error is given by for any exhaustion (Λ n ) n≥1 . Let us observe that, since dµ is a probability measure, the L ∞ (U, V ) norm controls the L 2 (U, V ) norm, and thus the above holds whenever the map u is uniformly bounded over U . For the L ∞ norms, consider an expansion (2.1) where the functions φ ν : U → R are normalized such that φ ν L ∞ (U) = 1, for all ν ∈ F . Then ( u ν V ) ν∈F ∈ ℓ 1 (F ), and it is easily checked that, whenever the expansion (2.1) converges conditionally to a function u in L ∞ (U, V ), it also converges unconditionally to u in L ∞ (U, V ). In addition, for any exhaustion (Λ n ) n≥1 , we have the error estimate (2.5) The above estimate is simply obtained by triangle inequality, and therefore generally it is not as sharp as (2.4). One particular situation is when (φ ν ) ν∈F is an orthogonal basis of L 2 (U ) normalized in L ∞ . Then, if u ∈ L 2 (U, V ) and the satisfy ( u ν V ) ν∈F ∈ ℓ 1 (F ), we find on the one hand that (2.1) converges unconditionally to a limit in L ∞ (U, V ) and in turn in L 2 (U, V ). On the other hand, we know that it converges toward u ∈ L 2 (U, V ). Therefore, its limit in L ∞ (U, V ) is also u.
A crucial issue is the choice of the sets Λ n that we decide to use when defining the n-term truncation (2.2). Ideally, we would like to use the set Λ n which minimizes the truncation error in some given norm · among all sets of cardinality n.
In the case of the L 2 error, if (φ ν ) ν∈F is an orthonormal basis of L 2 (U ), the estimate (2.4) shows that the optimal Λ n is the set of indices corresponding to the n largest u ν V . This set is not necessarily unique, in which case any realization of Λ n is optimal.
In the case of the L ∞ error, there is generally no simple description of the optimal Λ n . However, when the φ ν are normalized in L ∞ (U ), the right-hand side in the estimate (2.5) provides an upper bound for the truncation error. This bound is minimized by again taking for Λ n the set of indices corresponding to the n largest u ν V , with the error now bounded by the ℓ 1 tail of the sequence ( u ν V ) ν∈F , in contrast to the ℓ 2 tail which appears in (2.4).
The properties of a given sequence (c ν ) ν∈F which ensure a certain rate of decay n −s of its ℓ q tail after one retains its n largest entries are well understood. Here, we use the following result, see [12], originally invoked by Stechkin in the particular case q = 2. This result says that the rate of decay is governed by the ℓ p summability of the sequence for values of p smaller than q. In view of (2.4) or (2.5), application of the above result shows that ℓ p summability of the sequence ( u ν V ) ν∈F implies a convergence rate n −s when retaining the terms corresponding to the n largest u ν V in (2.1). From (2.4), when (φ ν ) ν∈F is an orthonormal basis, we obtain s = 1 p − 1 2 if p < 2. From (2.5), when the φ ν are normalized in L ∞ (U ), we obtain s = 1 p − 1 if p < 1. In the present setting of polynomial approximation, we mainly consider four types of series corresponding to four different choices of the univariate functions φ k : • Taylor (or power) series of the form ν∈F t ν y ν , t ν := 1 ν! ∂ ν u(y = 0), ν! := j≥1 ν j !, (2.6) with the convention that 0! = 1.
• Legendre series of the form where (L k ) k≥0 is the sequence of Legendre polynomials on [−1, 1] normalized with respect to the uniform measure • Renormalized Legendre series of the form with (H k ) k≥0 being the sequence of Hermite polynomials normalized according to R |H k (t)| 2 g(t)dt = 1, and dγ given by (1.11). In this case U = R d and (H ν ) ν∈F is an orthonormal basis of L 2 (U, dγ).
We may therefore estimate the L 2 error resulting from the truncation of the Legendre series (2.7) by application of (2.4), or the L ∞ error resulting from the truncation of the Taylor series (2.6) or renormalized Legendre series (2.8) by application of (2.5). According to Lemma 2.1, we derive convergence rates that depend on the value of p such that of the coefficient In a series of recent papers such summability results have been obtained for various types of parametric PDEs. We refer in particular to [13,1] for the elliptic PDE (1.3) with affine parameter dependence (1.4), to [21,2] for the lognormal dependence (1.10), and to [9] for more general PDEs and parameter dependence. One specific feature is that these conditions can be fulfilled in the infinite-dimensional framework. We thus obtain convergence rates that are immune to the curse of dimensionality, in the sense that they hold with d = ∞. Here, we mainly discuss the results established in [1,2] which have the specificity of taking into account the support properties of the functions ψ j .
One problem with this approach is that the sets Λ n associated to the n largest values in these sequences are generally not downward closed. In the next sections, we revisit these results in order to establish similar convergence rates for approximation in downward closed polynomial spaces.

Summability Results
The summability results in [1,2] are based on certain weighted ℓ 2 estimates which can be established for the previously defined coefficient sequences under various relevant conditions for the elliptic PDE (1.3). We first report below these weighted estimates. The first one from [1] concerns the affine parametrization (1.4). Here, we have V = H 1 0 (D) and V ′ denotes its dual H −1 (D).
for some fixed number r > 0. Then, one has The constants bounding these sums depend on r, f V ′ , a min and a L ∞ .
A few words are in order concerning the proof of these estimates. The first estimate (2.11) is established by first proving that the uniform ellipticity assumption (1.5) implies the ℓ 2 summability of the Taylor sequence ( t ν V ) ν∈F . Since the assumption (2.10) means that (1.5) holds with the ψ j replaced by ρ j ψ j , this gives the ℓ 2 summability of the Taylor sequence for the renormalized map y → u(ρy), ρy = (ρ j y j ) j≥1 , which is equivalent to (2.11). The second estimate is established by first showing that j≥1 ρ j |ψ j | ≤ a − r implies finiteness of the weighted Sobolev-type norm Then, one uses the Rodrigues formula L k (t) = d dt k √ 2k+1 k! 2 k (t 2 − 1) k in each variable y j to bound the weighted ℓ 2 sum in (2.12) by this norm. Remark 2.1 As shown in [1], the above result remains valid for more general classes of orthogonal polynomials of Jacobi type, such as the Chebyshev polynomials which are associated with the univariate measure The second weighted ℓ 2 estimate from [2] concerns the lognormal parametrization (1.10).
Theorem 2.2 Let r ≥ 0 be an integer. Assume that there exists a positive sequence Then, one has where with the convention that k l = 0 when l > k. The constant bounding this sum depends on f V ′ , j≥1 exp(−ρ 2 j ) and on the difference C r − K.
Similar to the weighted ℓ 2 estimate (2.12) for the Legendre coefficients, the proof of (2.14) follows by first establishing finiteness of a weighed Sobolev-type norm under the assumption (2.13) in the above theorem. Then one uses the Rodrigues g(t) , with g given by (1.11), in each variable y j to bound the weighted ℓ 2 sum in (2.14) by this norm.
In summary, the various estimates expressed in the above theorems all take the form Then, one natural strategy for establishing ℓ p summabilty of the sequence (c ν ) ν∈F is to invoke Hölder's inequality, which gives, for all 0 < p < 2, where the sequence (κ ν ) ν∈F is defined by Therefore ℓ p summability of (c ν ) ν∈F follows from ℓ q summability of (κ ν ) ν∈F with 0 < q < ∞ such that 1 q = 1 p − 1 2 . This ℓ q summability can be related to that of the univariate sequence one readily obtains the following elementary result, see [12] for more details.
In the case ω ν = ρ ν , i.e. κ ν = b ν , this shows that the ℓ p summability of the Taylor coefficients ( t ν V ) ν∈F follows if the assumption (2.10) holds with b = (ρ −1 j ) j≥1 ∈ ℓ q (N) and ρ j > 1 for all j. By a similar factorization, it is also easily checked that for any algebraic factor of the form α(ν) This allows us to reach a similar conclusion in the cases ω ν = β(ν) −1 ρ ν or ω ν = β(ν) −2 ρ ν , which correspond to the Legendre coefficients ( w ν V ) ν∈F and ( w ν V ) ν∈F , in view of (2.12). Likewise, in the case where ω ν = ξ 1/2 ν , using the factorization , it is shown in [2] that the sum on the left converges if b ∈ ℓ q , provided that r was chosen large enough such that q > 2 r . This shows that the ℓ p summability of the Hermite coefficients ( h ν V ) ν∈F follows if the assumption (2.13) holds with b = (ρ −1 j ) j≥1 ∈ ℓ q (N). Note that, since the sequence b can be renormalized, we may replace (2.13) by the condition sup x∈D j≥1 without a specific bound.

Approximation by Downward Closed Polynomials
The above results, combined with Lemma 2.1, allow us to build polynomial approximations u Λn with provable convergence rates n −s in L ∞ or L 2 by n-term trunctation of the various polynomial expansions. However, we would like to obtain such convergence rates with sets Λ n that are in addition downward closed. Notice that if a sequence (κ ν ) ν∈F of nonnegative numbers is monotone nonincreasing, that is then the set Λ n corresponding to the n largest values of κ ν (up to a specific selection in case of equal values) is downward closed. More generally, there exists a sequence (Λ n ) n≥1 of downward closed realizations of such sets which is nested, i.e.
We are now in position to state a general theorem that gives a condition for approximation using downward closed sets in terms of weighted ℓ 2 summability.
Proof 2.1 With ( κ ν ) ν∈F being the monotone majorant of (κ ν ) ν∈F , we observe that We pick a nested sequence (Λ n ) n≥1 of downward closed sets, such that Λ n consists of the indices corresponding to the n largest κ ν . Denoting by ( κ n ) n≥1 the decreasing rearrangement of ( κ ν ) ν∈F , we observe that With p such that 1 where we have used Hölder's inequality and the properties of ( κ n ) n≥1 . This gives (2.16) with C := AB.
We now would like to apply the above result with In the case of the Taylor series, where κ ν = b ν , we readily see that if b j < 1 for all j ≥ 1, then the sequence (κ ν ) ν∈F is monotone nonincreasing, and therefore Lemma 2.2 shows that b ∈ ℓ q implies (κ ν ) ν∈F ∈ ℓ q m (F ). By application of Theorem 2.3 with the value r = 1, this leads to the following result.
where u Λn is the truncated Taylor series and Λ n is any downward closed set corresponding to the n largest κ ν .
In the case of the Legendre series, the weight κ ν = b ν β(ν) is not monotone nonincreasing due to the presence of the algebraic factor β(ν). However, the following result holds. Lemma 2.3 For any 0 < q < ∞ and for any algebraic factor of the form α(ν) :=

Proof 2.2
The implication from right to left is a consequence of Lemma 2.2, and so we concentrate on the implication from left to right. For this it suffices to find a majorant κ ν of κ ν := α(ν)b ν which is monotone nonincreasing and such that ( κ ν ) ν∈F ∈ ℓ q (F ). We notice that for any τ > 1, there For some J ≥ 1 and τ to be fixed further, we may thus write The first product is bounded by (1 − θ q ) −J . Each factor in the second product is a converging series which is bounded by 1 + cb q j for some c > 0 that depends on c 1 , c 2 and b ℓ ∞ . It follows that this second product converges. Therefore ( κ ν ) ν∈F belongs to ℓ q (F ).
Finally, we show that κ ν is monotone nonincreasing provided that J is chosen large enough. It suffices to show that κ ν+ej ≤ κ ν for all ν ∈ F and for all j ≥ 1 where e j := (0, . . . , 0, 1, 0, . . . ), is the Kronecker sequence of index j. When j ≤ J this is obvious since κ ν+ej = τ b j κ ν ≤ θ κ ν ≤ κ ν . When j > J, we have Noticing that the sequence a n := 1+c1(n+1) 1+c1n c2 converges toward 1 and is therefore bounded, and that b j tends to 0 as j → +∞, we find that for J sufficiently large, the right-hand side in the above equation is bounded by 1 for all ν and j > J.
From Lemma 2.3, by applying Theorem 2.3 with r = 1 or r = 2, we obtain the following result.
where u Λn is the truncated Legendre series and Λ n is any downward closed set corresponding to the n largest κ ν where κ ν := b ν β(ν). If q < 2, we also have with Λ n any downward closed set corresponding to the n largest κ ν where κ ν : Finally, in the case of the Hermite coefficients, which corresponds to the weight we can establish a similar summability result.

Proof 2.3
For any ν ∈ F and any k ≥ 1 we have and therefore the sequence (κ ν ) ν∈F is monotone nonincreasing. Now we check that (κ ν ) ν∈F ∈ ℓ q (F ), using the factorization (2.18) where the inequality follows from the fact that the value l = n ∧ r := min{n, r} is contained in the sum. The jth factor F j in the rightmost product in (2.18) may be written as where C r,q := n≥r n r −q/2 = (r!) q/2 n≥0 (n + 1) · · · (n + r) −q/2 < ∞, (2.19) since we have assumed that q > 2/r. This shows that each F j is finite. If b ∈ ℓ q (N), there exists an integer J ≥ 0 such that b j < 1 for all j > J. For such j, we can bound F j by 1 + (C r,q + r − 1)b q j , which shows that the product converges.
From this lemma, and by application of Theorem 2.3 with the value r = 2, we obtain the following result for the Hermite series.
where u Λn is the truncated Hermite series and Λ n is a downward closed set corresponding to the n largest κ ν given by (2.17).
In summary, we have established convergence rates for approximation by downward closed polynomial spaces of the solution map (1.2) associated to the elliptic PDE (1.3) with affine or lognormal parametrization. The conditions are stated in terms of the control on the L ∞ norm of j≥1 ρ j |ψ j |, where the ρ j have a certain growth measured by the ℓ q summability of the sequence b = (b j ) j≥1 = (ρ −1 j ) j≥1 . This is a way to quantify the decay of the size of the ψ j , also taking their support properties into account, and in turn to quantify the anisotropic dependence of u(y) on the various coordinates y j . Other similar results have been obtained with different PDE models, see in particular [12]. In the above results, the polynomial approximants are constructed by truncation of infinite series. The remainder of the paper addresses the construction of downward closed polynomial approximants from evaluations of the solution map at m points {y 1 , . . . , y m } ∈ U , and discusses the accuracy of these approximants.

Sparse Interpolation by Downward Closed Polynomials
Interpolation is one of the most standard processes for constructing polynomial approximations based on pointwise evaluations. Given a downward closed set Λ ⊂ F of finite cardinality, and a set of points we would like to build an interpolation operator I Λ , that is, I Λ u ∈ V Λ is uniquely characterized by I Λ u(y) = u(y), y ∈ Γ, for any V -valued function u defined on U .
In the univariate case, it is well known that such an operator exists if and only if Γ is a set of pairwise distinct points, and that additional conditions are needed in the multivariate case. Moreover, since the set Λ may come from a nested sequence (Λ n ) n≥1 as discussed in Section 2, we are interested in having similar nestedness properties for the corresponding sequence (Γ n ) n≥1 , where Such a nestedness property allows us to recycle the n evaluations of u which have been used in the computation of I Λn u, and use only one additional evaluation for the next computation of I Λn+1 u.
It turns out that such hierarchical interpolants can be constructed in a natural manner by making use of the downward closed structure of the index sets Λ n . This construction is detailed in [8] but its main principles can be traced from [23]. In order to describe it, we assume that the parameter domain is of either form with the convention that d = ∞ in the second case. However, it is easily check that the construction can be generalized in a straightforward manner to any domain with Cartesian product form where the J k are finite or infinite intervals. We start from a sequence of pairwise distinct points We denote by I k the univariate interpolation operator on the space V k := V ⊗ P k associated with the k-section {t 0 , . . . , t k } of this sequence, that is, for any V -valued function u defined on [−1, 1]. We express I k in the Newton form and set I −1 = 0 so that we can also write Obviously the difference operator ∆ k annihilates the elements of V k−1 . In addition, since ∆ k u(t j ) = 0 for j = 0, . . . , k − 1, we have The coefficient α k ∈ V can be computed inductively, since it is given by that is, the interpolation error at t k when using I k−1 . Setting we observe that the system {B 0 , . . . , B k } is a basis for P k . It is sometimes called a hierarchical basis.
In the multivariate setting, we tensorize the grid T , by defining We first introduce the tensorized operator recalling that the application of a tensorized operator ⊗ j≥1 A j to a multivariate function amounts in applying each univariate operator A j by freezing all variables except the jth one, and then applying A j to the unfrozen variable. It is readily seen that I ν is the interpolation operator on the tensor product polynomial space associated to the grid of points This polynomial space corresponds to the particular downward closed index set of rectangular shape Λ = R ν := { ν : ν ≤ ν}.
Defining in a similar manner the tensorized difference operators we observe that The following result from [8] shows that the above formula can be generalized to any downward closed set in order to define an interpolation operator. We recall its proof for sake of completeness.
Theorem 3.1 Let Λ ⊂ F be a finite downward closed set, and define the grid Then, the interpolation operator onto V Λ for this grid is defined by Proof 3.1 From the downward closed set property, V ν ⊂ V Λ for all ν ∈ Λ. Hence the image of I Λ is contained in V Λ . With I Λ defined by (3.2), we may write for any ν ∈ Λ. Since y ν ∈ Γ ν , we know that On the other hand, if ν ν, this means that there exists a j ≥ 1 such that ν j > ν j . For this j we thus have ∆ ν u(y) = 0 for all y ∈ U with the jth coordinate equal to t νj by application of ∆ νj in the jth variable, so that ∆ ν u(y ν ) = 0.
The decomposition (3.2) should be viewed as a generalization of the Newton form (3.1). In a similar way, its terms can be computed inductively: if Λ = Λ ∪ {ν} where Λ is a downward closed set, we have Therefore, if (Λ n ) n≥1 is any nested sequence of downward closed index sets, we can compute I Λn by n iterations of Note that (B ν ) ν∈Λ is a basis of P Λ and that any f ∈ V Λ has the unique decom- where the coefficients α ν = α ν (f ) ∈ V are defined by the above procedure. Also note that α ν (f ) does not depend on the choice of Λ but only on ν and f .

Stability and Error Estimates
The pointwise evaluations of the function u could be affected by errors, as modeled by (1 .6) and (1.7). The stability of the interpolation operator with respect to such perturbations is quantified by the Lebesgue constant, which is defined by where the supremum is taken over the set of all V -valued functions f defined everywhere and uniformly bounded over U . It is easily seen that this supremum is in fact independent of the space V , so that we may also write where the supremum is now taken over real-valued functions. Obviously, we have where ε is the noise level from (1.7). The Lebesgue constant also allows us to estimate the error of interpolation u − I Λ u L ∞ (U,V ) for the noiseless solution map in terms of the best approximation error in the L ∞ norm: for any u ∈ L ∞ (U, V ) and any u ∈ V Λ we have We have seen in Section 2 that for relevant classes of solution maps y → u(y), there exist sequences of downward closed sets (Λ n ) n≥1 with #(Λ n ) = n, such that inf u∈VΛ n u − u L ∞ (U,V ) ≤ Cn −s , n ≥ 1, for some s > 0. For such sets, we thus have This motivates the study of the growth of L Λn as n → +∞. For this purpose, we introduce the univariate Lebesgue constants Note that L 0 = 1. We also define an analog quantity for the difference operator In the particular case of the rectangular downward closed sets Λ = R ν , since I Λ = I ν = ⊗ j≥1 I νj , we have Therefore, if the sequence T = (t k ) k≥0 is such that for some θ ≥ 1, we find that for all ν ∈ F . For arbitrary downward closed sets Λ, the expression of I Λ shows that Therefore, if the sequence T = (t k ) k≥0 is such that we find that The following result from [8] shows that this general estimate is also valid under the assumption (3.4) on the growth of L k .
One noticeable feature of the above result is that the bound on L Λ only depends on #(Λ), independently of the number of variables, which can be infinite, as well as of the shape of Λ.
We are therefore interested in univariate sequences T = (t k ) k≥0 such that L k and D k have moderate growth with k. For Chebyshev or Gauss-Lobatto points, given by C k := cos 2l + 1 2k + 2 π : l = 0, . . . , k and G k := cos l k π : l = 0, . . . , k , it is well known that the Lebesgue constant has logarithmic growth L k ∼ ln(k), thus slower than algebraic. However these points are not the k section of a single sequence T , and therefore they are not convenient for our purposes. Two examples of univariate sequences of interest are the following.
• The Leja points: from an arbitrary t 0 ∈ [−1, 1] (usually taken to be 1 or 0), this sequence is recursively defined by Note that this choice yields hierarchical basis functions B k that are uniformly bounded by 1. Numerical computations of L k for the first 200 values of k indicates that the linear bound |e − e l | : |e| ≤ 1 .
These points have the property of accumulating in a regular manner on the unit circle according to the so-called Van der Corput enumeration [4]. It is proven in [5] that the linear bound (3.6) holds for the Lebesgue constant of the complex interpolation operator on the unit disc associated to these points. The sequence of real parts t k := ℜ(e k ), is defined after eliminating the possible repetitions corresponding to e k = e l for two different values of k = l. These points coincide with the Gauss-Lobatto points for values of k of the form 2 n + 1 for n ≥ 0. A quadratic bound is established in [6].
If we use such sequences, application of Theorem 3.2 gives bounds of the form for example with θ = 2 when using the ℜ-Leja points, or θ = 1 when using the Leja points provided that the conjectured bound (3.6) holds. Combining with (3.3), we obtain the convergence estimate which reveals a serious deterioration of the convergence rate when using interpolation instead of truncated expansions. However, for the parametric PDE models discussed in Section 2, it is possible to show that this deterioration actually does not occur, based on the following lemma which relates the interpolation error to the summability of coefficient sequences in general expansions of u.
Lemma 3.1 Assume that u admits an expansion of the type (2.1), where φ ν L ∞ (U) ≤ 1 which is unconditionally convergent towards u in L ∞ (U, V ). Assume in addition that y → u(y) is continuous from U equipped with the product topology toward V . If the univariate sequence T = (t k ) k≥0 is such that (3.4) or (3.5) holds for some θ ≥ 1, then, for any downward closed set Λ, (3.7) Proof 3.2 The unconditional convergence of (2.1) and the continuity of u with respect to the product topology allow us to say that the equality in (2.1) holds everywhere in U . We may thus write where we have used that I Λ φ ν = φ ν for every ν ∈ Λ since φ ν ∈ P Λ . For the second sum on the right-hand side, we observe that for each ν / ∈ Λ, since ∆ ν annihilates P ν whenever ν ≤ ν. Therefore where I stands for the identity operator. This implies Since (3.4) or (3.5) holds, we obtain from Theorem 3.2 that which yields (3.7).

Discrete Least Squares on V -valued Linear Spaces
Least-squares fitting is an alternative approach to interpolation for building a polynomial approximation of u from V Λ . In this approach we are given m observations u 1 , . . . , u m of u at points y 1 , . . . , y m ∈ U ⊆ R d where m ≥ n = #(Λ). We first discuss the least-squares method in the more general setting of V -valued linear spaces, where Y n is the space of real-valued functions defined everywhere on U such that dim(Y n ) = n. In the next section, we discuss more specifically the case where Y n = P Λ . Here we study the approximation error in the L 2 (U, V, dµ) norm for some given probability measure dµ, when the evaluation points y i are independent and drawn according to this probability measure. For notational simplicity we use the shorthand The least-squares method selects the approximant of u in the space V n as In the noiseless case where u i := u(y i ) for any i = 1, . . . , m, this also writes where the discrete seminorm is defined by .
Note that f 2 m is an unbiased estimator of f 2 since we have Let {φ 1 , . . . , φ n } denote an arbitrary L 2 (U, dµ) orthonormal basis of the space Y n . If we expand the solution to (4.1) as n j=1 c j φ j , with c j ∈ V , the V -valued vector c = (c 1 , . . . , c n ) t is the solution to the normal equations where the matrix G has entries and where the V -valued data vector d = (d 1 , . . . , d n ) t is given by This linear system always has at least one solution, which is unique when G is nonsingular. When G is singular, we may define u L as the unique minimal ℓ 2 (R n , V ) norm solution to (4.2).
In the subsequent analysis, we sometimes work under the assumption of a known uniform bound We introduce the truncation operator and notice that it is a contraction: The truncated least-squares approximation is defined by Note that, in view of (4.3), we have u(y)− u T (y) V ≤ u(y)− u L (y) V for any y ∈ U and therefore u − u T ≤ u − u L .
Note that the random matrix G concentrates toward its expectation which is the identity matrix I as m → ∞. In other words, the probability that G is ill-conditioned becomes very small as m increases. The truncation operator aims at avoiding instabilities which may occur when G is ill-conditioned. As an alternative proposed in [15], we may define for some given A > 1 the conditioned least-squares approximation by where cond(G) := λ max (G)/λ min (G) is the usual condition number. The property that G − I 2 ≤ δ for some 0 < δ < 1 amounts to the norm equivalence It is well known that if m ≥ n is too much close to n, least-squares methods may become unstable and inaccurate for most sampling distributions. For example, if U = [−1, 1] and Y n = P n−1 is the space of algebraic polynomials of degree n − 1, then with m = n the estimator coincides with the Lagrange polynomial interpolation which can be highly unstable and inaccurate, in particular for equispaced points. Therefore, m should be sufficiently large compared to n for the probability that G is ill-conditioned to be small. This trade-off between m and n has been analyzed in [11], using the function which is the diagonal of the integral kernel of the L 2 (U, dµ) projector on Y n . This function depends on dµ, but not on the chosen orthonormal basis. It is strictly positive in U under minimal assumptions on the orthonormal basis, for example if one element of the basis is the constant function over all U . Obviously, the function k n satisfies U k n dµ = n.
We define K n := k n L ∞ (U) ≥ n.
The following results for the least-squares method with noiseless evaluations were obtained in [11,27,7,15] for real-valued functions, however their proof extends in a straightforward manner to the present setting of V -valued functions. They are based on a probabilistic bound for the event G − I 2 > δ using the particular value δ = 1 2 , or equivalently the value A = 1+δ 1−δ = 3 as a bound on the condition number of G. then the following hold.
(iii) The conditioned least-squares estimator satisfies, in the noiseless case, where ζ(m) is as in (ii).
(iv) If u satisfies (4.3), then the estimator u E ∈ {u L , u T , u C } satisfies, in the noiseless case, with probability larger than 1 − 2m −r .
In the case of noisy evaluations modeled by (1.6)-(1.7), the observations are given by u i = u(y i ) + η(y i ). (4.6) The following result from [7] shows that (4.5) holds up to this additional perturbation.
with probability larger than 1 − 2n −r , where ε is the noise level in (1.7).
Similar results, with more general assumptions on the type of noise, are proven in [11,28,15].

Downward Closed Polynomial Spaces and Weighted Least Squares
Condition (4.4) shows that K n gives indications on the number m of observations required to ensure stability and accuracy of the least-squares approximation. In order to understand how demanding this condition is with respect to m, it is important to have sharp upper bounds for K n . Such bounds have been proven when the measure dµ on U = [−1, 1] d has the form where θ 1 , θ 2 > −1 are real shape parameters and C is a normalization constant such that U dµ = 1. Sometimes (4.7) is called the Jacobi measure, because the Jacobi polynomials are orthonormal in L 2 (U, dµ). Remarkable instances of the measure (4.7) are the uniform measure, when θ 1 = θ 2 = 0, and the Chebyshev measure, when θ 1 = θ 2 = − 1 2 . When Y n = P Λ is a multivariate polynomial space and Λ is a downward closed multi-index set with #(Λ) = n, it is proven in [7,26] that K n satisfies an upper bound which only depends on n and on the choice of the measure (4.7) through the values of θ 1 and θ 2 .
Lemma 4.1 Let dµ be the measure defined in (4.7). Then it holds A remarkable property of both algebraic upper bounds in (4.8) is that the exponent of n is independent of the dimension d, and of the shape of the downward closed set Λ. Both upper bounds are sharp in the sense that equality holds for multi-index sets of rectangular type Λ = R ν corresponding to tensor product polynomial spaces.
As an immediate consequence of Theorem 4.1 and Lemma 4.1, we have the next corollary. Other types of results on the accuracy of least squares have been recently established in [14], under conditions of the same type as (4.9).
In some situations, for example when n is very large, the conditions (4.9) might require a prohibitive number of observations m. It is therefore a legitimate question to ask whether there exist alternative approaches with less demanding conditions than (4.9) between m and n. At best, we would like that m is of order only slightly larger than n, for example by a logarithmic factor. In addition, the above analysis does not apply to situations where the basis functions φ k are unbounded, such as when using Hermite polynomials in the expansion (2.9). It is thus desirable to ask for the development of approaches that also cover this case.
These questions have an affirmative answer by considering weighted least-squares methods, as proposed in [19,22,15]. In the following, we survey some results from [15]. For the space V n = V ⊗ Y n , the weighted least-squares approximation is defined as for some given choice of weights w i ≥ 0. This estimator is again computed by solving a linear system of normal equations now with the matrix G with entries Of particular interest to us are weights of the form where w is some nonnegative function defined on U such that U w −1 dµ = 1. (4.10) We then denote by dσ the probability measure and we draw the independent points y 1 , . . . , y m from dσ. The case w ≡ 1 and dσ = dµ corresponds to the previously discussed standard (unweighted) least-squares estimator u L . As previously done for u L , we associate to u W a truncated estimator u T and a conditioned estimator u C , by replacing u L with u W in the corresponding definitions. Let us introduce the function y → k n,w (y) := n j=1 w(y)|φ j (y)| 2 , where once again {φ 1 , . . . , φ n } is an arbitrary L 2 (U, dµ) orthonormal basis of the space Y n . Likewise, we define K n,w := k n,w L ∞ (U) .
The following result, established in [15] for real-valued functions, extends Theorem 4.1 to this setting. Its proof in the V -valued setting is exactly similar. The choice (4.12) also gives K n,w = k n,w L ∞ (U) = n, and leads to the next result, as a consequence of the previous theorem. then the same conclusions of Theorem 4.1 hold true with u L replaced by u W , with w given by (4.12) and the weights taken as w i = w(y i ).
The above theorem ensures stability and accuracy of the weighted least-squares approximation, under the minimal condition that m is linearly proportional to n, up to a logarithmic factor. Clearly this is an advantage of weighted least squares compared to standard least squares, since condition (4.4) is more demanding than (4.14) in terms of the number of observations m.
However, this advantage comes with some drawbacks that we now briefly recall, see [15] for an extensive description. In general (4.11) and (4.13) are not product measures, even if dµ is one. Therefore, the first drawback of using weighted least squares concerns the efficient generation of independent samples from multivariate probability measures, whose computational cost could be prohibitively expensive, above all when the dimension d is large. In some specific settings, for example downward closed polynomial spaces Y n = P Λ with #(Λ) = n, and when dµ is a product measure, this drawback can be overcome. We refer to [15], where efficient sampling algorithms have been proposed and analyzed. For any m and any downward closed set Λ, these algorithms generate m independent samples with proven bounds on the required computational cost. The dependence on the dimension d and m of these bounds is linear. For the general measure (4.11) the efficient generation of the sample is a nontrivial task, and remains a drawback of such an approach.
The second drawback concerns the use of weighted least squares in a hierarchical context, where we are given a nested sequence Λ 1 ⊂ . . . ⊂ Λ n of downward closed sets, instead of a single such set Λ. Since the measure (4.13) depends on n, the sets (Λ n ) n≥1 are associated to different measures (dσ n ) n≥1 . Hence, recycling samples from the previous iterations of the adaptive algorithm is not as straighforward as in the case of standard least squares.
As a final remark, let us stress that the above results of Theorem 4.3 and Theorem 4.4 hold for general approximation spaces Y n other than polynomials.

Selection of Downward Closed Polynomial Spaces
The interpolation and least-squares methods discussed in Section 3 and Section 4 allow us to construct polynomial approximations in V Λ = V ⊗ P Λ of the map (1.2) from its pointwise evaluations, for some given downward closed set Λ. For these methods, we have given several convergence results in terms of error estimates either in L ∞ (U, V ) or L 2 (U, V, dµ). In some cases, these estimates compare favorably with the error of best approximation min u∈VΛ u − u measured in such norms.
A central issue which still needs to be addressed is the choice of the downward closed set Λ, so that this error of best approximation is well behaved, for a given map u. Ideally, for each given n, we would like to use the set where D n is the family of all downward closed sets Λ of cardinality n. However such sets Λ n are not explicitely given to us, and in addition the resulting sequence (Λ n ) n≥1 is generally not nested.
Concrete selection strategies aim to produce "suboptimal yet good" nested sequences (Λ n ) n≥1 different from the above. Here, an important distinction should be made between nonadaptive and adaptive selection strategies.
In nonadaptive strategies, the selection of Λ n is made in an a-priori manner, based on some available information on the given problem. The results from Section 2.3 show that, for relevant instances of solution maps associated to parametric PDEs, there exist nested sequences (Λ n ) n≥1 of downward closed sets such that #(Λ n ) = n and min u∈VΛ n u − u decreases with a given convergence rate n −s as n → ∞. In addition, these results provide constructive strategies for building the sets Λ n , since these sets are defined as the indices associated to the n largest κ ν := max ν≥ν κ ν like in Theorem 2.5, or directly to the n largest κ ν like in Theorem 2.4 and Theorem 2.6, and since the κ ν are explicitely given numbers.
In the case where we build the polynomial approximation by interpolation, Theorem 3.3 shows that a good choice of Λ n is produced by taking κ ν to be either π(ν)b ν or π(ν)β(ν) 2 b ν where b = (ρ −1 j ) j≥1 is such that (2.10) holds. In the case where we build the polynomial approximation by least-squares methods, the various results from Section 4 show that under suitable assumptions, the error is nearly as good as that of best approximation in L 2 (U, V, dµ) with respect to the relevant probability measure. In the affine case, Theorem 2.5 shows that a good choice of Λ n is produced by taking κ ν to be b ν β(ν) where b = (ρ −1 j ) j≥1 is such that (2.10) holds. In the lognormal case Theorem 2.6 shows that a good choice of Λ n is produced by taking κ ν to be given by (2.17) where b = (ρ −1 j ) j≥1 is such that (2.15) holds. Let us briefly discuss the complexity of identifying the downward closed set Λ n associated to the n largest κ ν . For this purpose, we introduce for any downward closed set Λ its set of neighbors defined by We may in principle define Λ n = {ν 1 , . . . , ν n } by the following induction.
In the finite-dimensional case d < ∞, we observe that N (Λ k ) is contained in the union of N (Λ k−1 ) with the set consisting of the indices where e j is the Kroenecker sequence with 1 at position j. As a consequence, since the values of the κ ν have already been computed for ν ∈ N (Λ k−1 ), the step k of the induction requires at most d evaluations of κ ν , and therefore the overall computation of Λ n requires at most nd evaluations.
In the infinite-dimensional case d = ∞, the above procedure cannot be practically implemented, since the set of neighbors has infinite cardinality. This difficulty can be circumvented by introducing a priority order among the variables, as done in the next definitions. Obviously, if (c ν ) ν∈F is anchored, one of the sets Λ n corresponding to its n largest values is anchored. It is also readily seen that all sequences ( κ ν ) ν∈F that are used in Theorems 2.4, 2.5, 2.6 or 3.3 for the construction of Λ n are anchored, provided that the sequence b = (ρ −1 j ) j≥1 is monotone nonincreasing. This is always the case up to a rearrangement of the variables. For any anchored set Λ, we introduce the set of its anchored neighbors defined by where j(Λ) := max{j : ν j > 0 for some ν ∈ Λ}.
We may thus modify in the following way the above induction procedure.
This procedure is now feasible in infinite dimension. At each step k the number of active variables is limited by j(Λ k ) ≤ k − 1, and the total number of evaluations of κ ν needed to construct Λ n does not exceed 1 + 2 + · · · + (n − 1) ≤ n 2 /2. In adaptive strategies the sets Λ n are not a-priori selected, but instead they are built in a recursive way, based on earlier computations. For instance, one uses the previous set Λ n−1 and the computed polynomial approximation u Λn−1 to construct Λ n . If we impose that the sets Λ n are nested, this means that we should select an index ν n / ∈ Λ n−1 such that Λ n := Λ n−1 ∪ {ν n }.
The choice of the new index ν n is further limited to N (Λ n−1 ) if we impose that the constructed sets Λ n are downward closed, or to N (Λ n−1 ) if we impose that these sets are anchored. Adaptive methods are known to sometimes perform significantly better than their nonadaptive counterpart. In the present context, this is due to the fact that the apriori choices of Λ n based on the sequences κ ν may fail to be optimal. In particular, the guaranteed rate n −s based on such choices could be pessimistic, and better rates could be obtained by other choices. However, convergence analysis of adaptive methods is usually more delicate. We next give examples of possible adaptive strategies in the interpolation and least-squares frameworks.

Adaptive Selection for Interpolation
We first consider polynomial approximations obtained by interpolation as discussed in Section 3. The hierarchical form may formally be viewed as a truncation of the expansion of u in the hierarchical basis which however may not always be converging, in contrast to the series discussed in Section 2. Nevertheless, we could in principle take the same view, and use for Λ n the set of indices corresponding to the n largest terms of (5.2) measured in some given metric L p (U, V, dµ). This amounts in choosing the indices of the n largest w ν α ν V , where the weight w ν is given by This weight is easily computable when dµ is a tensor product measure, such as the uniform measure. In the case where p = ∞ and if we use the Leja sequence, we know that B ν L ∞ (U) = 1 and therefore this amounts to choosing the largest α ν V . This selection strategy is not practially feasible since we cannot afford this exhaustive search over F . However, it naturally suggests the following adaptive greedy algorithm, which has been proposed in [8].
In the case where p = ∞ and if we use the Leja sequence, this strategy amounts in picking the index ν n that maximizes the interpolation error u(y ν ) − I Λn−1 u(y ν ) V among all ν in N (Λ n−1 ). By the same considerations as previously discussed for the a-priori selection of Λ n , we find that in the finite-dimensional case, the above greedy algorithm requires at most dn evaluation after n steps. When working with infinitely many variables (y j ) j≥1 , we replace the infinite set N (Λ n ) in the algorithm by the finite set of anchored neighbors N (Λ n ) defined by (5.1). Running n steps of the resulting greedy algorithm requires at most n 2 /2 evaluations. Remark 5.1 A very similar algorithm has been proposed in [20] in the different context of adaptive quadratures, that is, for approximating the integral of u over the domain U rather than u itself. In that case, the natural choice is to pick the new index ν n that maximizes | U ∆ ν u dµ| over N (Λ n ) or N (Λ n ).
The main defect of the above greedy algorithm is that it may fail to converge, even if there exist sequences (Λ n ) n≥1 such that I Λn u converges toward u. Indeed, if ∆ ν u = 0 for a certain ν, then no index ν ≥ ν will ever be selected by the algorithm. As an example, if u is of the form u(y) = u 1 (y 1 )u 2 (y 2 ), where u 1 and u 2 are nonpolynomial smooth functions such that u 2 (t 0 ) = u 2 (t 1 ), then the algorithm could select sets Λ n with indices ν = (k, 0) for k = 0, . . . , n − 1, since the interpolation error at the point (t k , t 1 ) vanishes.
One way to avoid this problem is to adopt a more conservative selection rule which ensures that all of F is explored, by alternatively using the rule (5.3), or picking the multi-index ν ∈ N (Λ n ) which has appeared at the earliest stage in the neighbors of the previous sets Λ k . This is summarized by the following algorithm.
Even with such modifications, the convergence of the interpolation error produced by this algorithm is not generally guaranteed. Understanding which additional assumptions on u ensure convergence at some given rate, for a given univariate sequence T such as Leja points, is an open problem.
Remark 5.2 Another variant to the above algorithms consists in choosing at the iteration k more than one new index at a time within N (Λ k−1 ) or N (Λ k−1 ). In this case, we have n k := #(Λ k ) ≥ k. For example we may choose the smallest subset of indices that retains a fixed portion of the quantity ν∈Λ k−1 w ν α ν V . This type of modification turns out to be particularly relevant in the least-squares setting discussed in the next section.

Adaptive Selection for Least Squares
In this section we describe adaptive selections in polynomial spaces, for the leastsquares methods that have been discussed in Section 4. We focus on adaptive selection algorithms based on the standard (unweighted) least-squares method.
As a preliminary observation, it turns out that the most efficient available algorithms for adaptive selection of multi-indices might require the selection of more than one index at a time. Therefore, we adopt the notation that n k := #(Λ k ) ≥ k, where the index k denotes the iteration in the adaptive algorithm.
As discussed in Section 4, stability and accuracy of the least-squares approximation is ensured under suitable conditions between the number of samples and the dimension of the approximation space, see e.g. condition (4.9). Hence, in the development of reliable iterative algorithms, such conditions need to be satisfied at each iteration. When dµ is the measure (4.7) with shape parameters θ 1 , θ 2 , condition (4.9) takes the form of m k ln m k ≥ κ n s k , (5.5) where m k denotes the number of samples at iteration k, and Since n k increases with k, the minimal number of samples m k that satisfies (5.5) has to increase as well at each iteration. At this point, many different strategies can be envisaged for progressively increasing m k such that (5.5) remains satisfied at each iteration k. We may now present a first adaptive algorithm based on standard least squares.
• Set ν k := argmax Similarly to the previously discussed interpolation algorithms, in the case of infinitely many variables (y j ) j≥1 the set N (Λ k ) is infinite and should be replaced by the finite set of anchored neighbors N (Λ k ) defined by (5.1). As for interpolation, we may define a more conservative version of this algorithm in order to ensure that all of F is explored. For example, when k is even, we define ν k according to (5.6), and when k is odd we pick for ν k the multi-index ν ∈ N (Λ k ) which has appeared at the earliest stage in the neighbors of the previous sets Λ k . The resulting algorithm is very similar to the one presented for interpolation, with obvious modifications due to the use of least squares.
As announced at the beginning, it can be advantageous to select more than one index at a time from N (Λ k−1 ), at each iteration k of the adaptive algorithm. For describing the multiple selection of indices from N (Λ k−1 ), we introduce the so-called bulk chasing procedure. Given a finite set R ⊆ N (Λ k−1 ), a nonnegative function E : R → R and a parameter α ∈ (0, 1], we define the procedure bulk := bulk(R, E, α) that computes a set F ⊆ R of minimal positive cardinality such that ν∈F E(ν) ≥ α ν∈R E(ν).
A possible choice for the function E is where c ν is given from an available least-squares estimator that has been already computed on any downward closed set R ⊂ Λ ⊆ Λ k−1 ∪ N (Λ k−1 ). Another choice for E is where u L is the truncation to Λ k−1 of a least-squares estimator u L = ν∈Λ c ν φ ν that has been already computed on any downward closed set Λ k−1 ⊂ Λ ⊆ Λ k−1 ∪ N (Λ k−1 ), using a number of samples m k−1 that satisfies condition (5.5) with n k = #(Λ). The discrete norm in E M (ν) uses the same m k−1 evaluations of u that have been used to compute the least-squares approximation u L on Λ.
Both E L (ν) and E M (ν) should be viewed as estimators of the coefficient u, φ ν . The estimator E M (ν) is of Monte Carlo type and computationally cheap to calculate. Combined use of the two estimators leads to the next algorithm for greedy selection with bulk chasing, that has been proposed in [29].
• Assuming that Λ k−1 has been selected, set where E M uses the least-squares approximation u L = ν∈Λ c ν φ ν of u in V Λ that has been calculated at iteration k − 1 on a downward closed set Λ k−1 ⊂ Λ ⊆ Λ k−1 ∪ N (Λ k−1 ) using a number of samples m k−1 that satisfies (5.5) with n k = #(Λ).
• Compute the least-squares approximation of u on V Λ k−1 ∪F1 using a number of samples m k that satisfies (5.5) with n k = #(Λ k−1 ∪ F 1 ).
The set N (Λ k−1 ) can be large, and might contain many indices that are associated to small coefficients. Discarding these indices is important in order to avoid unnecessary computational burden in the calculation of the least-squares approximation. The purpose of the bulk procedure (5.7) is to perform a preliminary selection of a set F 1 ⊆ N (Λ k−1 ) of indices, using the cheap estimator E M . At iteration k, E M in (5.7) uses the estimator computed in (5.8) at iteration k − 1 and truncated to Λ k−1 . Afterwards, at iteration k, the least-squares approximation in (5.8) is calculated on Λ k−1 ∪ F 1 , using a number of samples m k which satisfies condition (5.5), with n k = #(Λ k−1 ∪ F 1 ). The second bulk procedure (5.9) selects a set F 2 of indices from F 1 , using the more accurate estimator E L . The convergence rate of the adaptive algorithm depends on the values given to the parameters α 1 and α 2 .
Finally we mention some open issues related to the development of adaptive algorithms using the weighted least-squares methods discussed in Section 4, instead of standard least squares. In principle the same algorithms described above can be used with the weighted least-squares estimator u W replacing the standard least-squares estimator u L , provided that, at each iteration k, the number of samples m k satisfies m k ln m k ≥ κ n k , and that the samples are drawn from the optimal measure, see Theorem 4.4. This ensures that at each iteration k of the adaptive algorithm, the weighted least-squares approximation remains stable and accurate. However, no guarantees on stability and accuracy are ensured if the above conditions are not met, for example when the samples from previous iterations are recycled.

Approximation in Downward Closed Spaces: beyond Polynomials
The concept of downward closed approximation spaces can be generalized beyond the polynomial setting. We start from a countable index set S equipped with a partial order ≤, and assume that there exists a root index 0 S ∈ S such that 0 S ≤ σ for all σ ∈ S. We assume that (B σ ) σ∈S is a basis of functions defined on [−1, 1] such that The set F is equipped with a partial order induced by its univariate counterpart: ν ≤ ν if and only if ν j ≤ ν j for all j ≥ 1. We may then define downward closed sets Λ ⊂ F in the same way as in Definition 1.1 which corresponds to the particular case S = N. We then define the associated downward closed approximation space by V Λ := V ⊗ B Λ , B Λ := span{B ν : ν ∈ Λ}, that is the space of functions of the form ν∈Λ u ν B ν with u ν ∈ V .
Then, if Λ ⊂ F is a downward closed set, we may define an interpolation operator I Λ onto V Λ associated to the grid Γ Λ := {y ν : ν ∈ Λ}.
In a similar manner as in the polynomial case, this operator is defined inductively by I Λ u := I Λ u + α ν B ν , α ν := α ν (u) = u(y ν ) − I Λ u(y ν ), where ν / ∈ Λ and Λ is any downward closed set such that Λ = Λ ∪ {ν}. We initialize this computation with Λ 1 = {0 F }, where 0 F is the null multi-index, by defining I Λ1 u as the constant function with value u(y 0F ).
The set S is thus a binary tree where λ −1 is the root node, (0, 0) is a child of λ 1 which is itself a child of λ −1 , every node (j, k) has two children (j + 1, 2k) and (j + 1, 2k + 1), and the relation λ ≤ λ means that λ is a parent of λ. The index j corresponds to the level of refinement, i.e. the depth of the node in the binary tree. We associate with S the sequence T := {t λ−1 , t λ1 , t (0,0) } ∪ t (j,k) := 2k + 1 2 j : (j, k) ∈ S, j ≥ 1 , where t λ−1 = −1, t λ1 = 1 and t (0,0) = 0. The hierarchical basis of piecewise linear functions defined over [−1, 1] is then given by where H(t) := max{0, 1 − |t|}, is the usual hat function. In dimension d = 1, the hierarchical interpolation amounts in the following steps: start by approximating f with the constant function equal to f (−1), then with the affine function that coincides with f at −1 and 1, then with the piecewise affine function that coincides with f at −1, 0 and −1; afterwards refine the approximation in further steps by interpolating f at the midpoint of an interval between two adjacents interpolation points.
Other relevant examples include piecewise polynomials, hierarchical basis functions, and more general interpolatory wavelets, see [10] for a survey.