Optimal weighted least-squares methods ∗

arXiv:1608.00512v1 [math.NA] 1 Aug 2016

Albert Cohen† and Giovanni Migliorati‡
August 2, 2016

Abstract
We consider the problem of reconstructing an unknown bounded function u defined on a domain X ⊂
Rd from noiseless or noisy samples of u at n points (xi )i=1,...,n . We measure the reconstruction error in a
norm L2 (X, dρ) for some given probability measure dρ. Given a linear space Vm with dim(Vm ) = m ≤ n,
we study in general terms the weighted least-squares approximations from the spaces Vm based on
independent random samples. It is well known that least-squares approximations can be inaccurate and
unstable when m is too close to n, even in the noiseless case. Recent results from [4, 5] have shown the
interest of using weighted least squares for reducing the number n of samples that is needed to achieve an
accuracy comparable to that of best approximation in Vm , compared to standard least squares as studied
in [3]. The contribution of the present paper is twofold. From the theoretical perspective, we establish
results in expectation and in probability for weighted least squares in general approximation spaces Vm .
These results show that for an optimal choice of sampling measure dµ and weight w, which depends on the
space Vm and on the measure dρ, stability and optimal accuracy are achieved under the mild condition
that n scales linearly with m up to an additional logarithmic factor. In contrast to [3], the present
analysis covers cases where the function u and its approximants from Vm are unbounded, which might
occur for instance in the relevant case where X = Rd and dρ is the Gaussian measure. From the numerical
perspective, we propose a sampling method which allows one to generate independent and identically
distributed samples from the optimal measure dµ. This method becomes of interest in the multivariate
setting where dµ is generally not of tensor product type. We illustrate this for particular examples of
approximation spaces Vm of polynomial type, where the domain X is allowed to be unbounded and high
or even infinite dimensional, motivated by certain applications to parametric and stochastic PDEs.

AMS classification numbers: 41A10, 41A25, 41A65, 62E17, 93E24.
Keywords: multivariate approximation, weighted least squares, error analysis, convergence rates, random matrices, conditional sampling, polynomial approximation.

1

Introduction

Let X be a Borel set of Rd . We consider the problem of estimating an unknown function u : X → R from
pointwise data (y i )i=1,...,n which are either noiseless or noisy observations of u at points (xi )i=1,...,n from X.
In numerous applications of interest, some prior information is either established or assumed on the function
u. Such information may take various forms such as:
∗ This

research is supported by Institut Universitaire de France and the ERC AdV project BREAD.
Universités, UPMC Univ Paris 06, CNRS, UMR 7598, Laboratoire Jacques-Louis Lions, 4, place Jussieu 75005,
Paris, France. email: cohen@ljll.math.upmc.fr
‡ Sorbonne Universités, UPMC Univ Paris 06, CNRS, UMR 7598, Laboratoire Jacques-Louis Lions, 4, place Jussieu 75005,
Paris, France. email: migliorati@ljll.math.upmc.fr
† Sorbonne

1

(i) regularity properties of u, in the sense that it belongs to a given smoothness class;
(ii) decay or sparsity of the expansion of u in some given basis;
(iii) approximability of u with some prescribed error by given finite-dimensional spaces.
Note that the above are often related to one another and sometimes equivalent, since many smoothness
classes can be characterized by prescribed approximation rates when using certain finite-dimensional spaces
or truncated expansions in certain bases.
This paper uses the third type of prior information, taking therefore the view that u can be “well
approximated” in some space Vm of functions defined everywhere on X, such that dim(Vm ) = m. We work
under the following mild assumption:
f or any x ∈ X, there exists v ∈ Vm such that v(x) 6= 0.

(1)

This assumption holds, for example, when Vm contains the constant functions. Typically, the space Vm
comes from a family (Vj )j≥1 of nested spaces with increasing dimension, such as algebraic or trigonometric
polynomials, or piecewise polynomial functions on a hierarchy of meshes.
We are interested in measuring the error in the L2 (X, dρ) norm
Z
1/2
kvk :=
|v|2 dρ
,
X

where dρ is a given probability measure on X. We denote by h·, ·i the associated inner product. One typical
strategy is to pick the estimate from a finite-dimensional space Vm such that dim(Vm ) = m. The ideal
estimator is given by the L2 (X, dρ) orthogonal projection of u onto Vm , namely
Pm u := argmin ku − vk.
v∈Vm

In general, this estimator is not computable from a finite number of observations. The best approximation
error
em (u) := min ku − vk = ku − Pm uk,
v∈Vm

thus serves as a benchmark for a numerical method based on a finite sample. In the subsequent analysis,
we make significant use of an arbitrary L2 (X, dρ) orthonormal basis {L1 , . . . , Lm } of the space Vm . We also
introduce the notation
em (u)∞ := min ku − vkL∞ ,
v∈Vm

∞

where L is meant with respect to dρ, and observe that em (u) ≤ em (u)∞ for any probability measure dρ.
The weighted least-squares method consists in defining the estimator as
n

uW := argmin
v∈Vm

1X i
w |v(xi ) − y i |2 ,
n i=1

(2)

where the weights wi > 0 are given. In the noiseless case y i = u(xi ), this also writes
argmin ku − vkn ,

(3)

v∈Vm

where the discrete seminorm is defined by
n

kvkn :=

1X i
w |v(xi )|2
n i=1
2

!1/2

.

(4)

This seminorm is associated with the semi-inner product h·, ·in . If we expand the solution to (3) as
Pm
j=1 vj Lj , the vector v = (vj )j=1,...,m is the solution to the normal equations
Gv = d,

(5)

where the matrix G has entries Gj,k = hLj , Lk in and where the data vector d = (dj )j=1,...,m is given
Pn
by dj := n1 i=1 wi y i Lj (xi ). This system always has at least one solution, which is unique when G is
nonsingular. When G is singular, we may define uW as the unique minimal ℓ2 norm solution to (5).
Note that G is nonsingular if and only if k · kn is a proper norm on the space Vm . Then, if the data are
noisefree that is, when y i = u(xi ), we may also write
n
uW = Pm
u,
n
where Pm
is the orthogonal projection onto Vm for the norm k · kn .
In practice, for the estimator (2) to be easily computable, it is important that the functions L1 , . . . , Lm
have explicit expressions that can be evaluated at any point in X so that the system (5) can be assembled.
Let us note that computing this estimator by solving (5) only requires that {L1 , . . . , Lm } is a basis of the
space Vm , not necessarily orthonormal in L2 (X, dρ). Yet, since our subsequent analysis of this estimator
makes use of an L2 (X, dρ) orthonormal basis, we simply assume that {L1 , . . . , Lm } is of such type.
In our subsequent analysis, we sometimes work under the assumption of a known uniform bound

kukL∞ ≤ τ.

(6)

We introduce the truncation operator
z 7→ Tτ (z) := sign(z) min{|z|, τ },
and we study the truncated weighted least-squares approximation defined by
uT := Tτ ◦ uW .
Note that, in view of (6), we have |u − uT | ≤ |u − uW | in the pointwise sense and therefore
ku − uT k ≤ ku − uW k.
The truncation operator aims at avoiding unstabilities which may occur when the matrix G is ill-conditioned.
In this paper, we use randomly chosen points xi , and corresponding weights wi = w(xi ), distributed in such
a way that the resulting random matrix G concentrates towards the identity I as n increases. Therefore, if
no L∞ bound is known, an alternative strategy consists in setting to zero the estimator when G deviates
from the identity by more than a given value in the spectral norm. We recall that for m × m matrices X,
this norm is defined as kXk2 := supkvk2 =1 kXvk2 . More precisely, we introduce the conditioned least-squares
approximation, defined by
(
uW , if kG − Ik2 ≤ 12 ,
uC :=
0,
otherwise.
The choice of 12 as a threshold for the distance between G and I in the spectral norm is related to our
subsequent analysis. However, the value 21 could be be replaced by any real number in ]0, 1[ up to some
minor changes in the formulation of our results. Note that
kG − Ik2 ≤

1
=⇒ cond(G) ≤ 3.
2
3

(7)

It is well known that if n ≥ m is too much close to m, weighted least-squares methods may become
unstable and inaccurate for most sampling distributions. For example, if X = [−1, 1] and Vm = Pm−1 is the
space of algebraic polynomials of degree m − 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.
The question that we want to address here in general terms is therefore:
Given a space Vm and a measure dρ, how to best choose the samples y i and weights wi in order to ensure that the L2 (X, dρ) error ku − ũk is comparable to em (u), with n being as close as possible to m, for
ũ ∈ {uW , uT , uC } ?
We address this question in the case where the xi are randomly chosen. More precisely, we draw independently the xi according to a certain probabiity measure dµ defined on X. A natural prescription for the
success of the method is that kvkn approaches kvk as n tends to +∞. Therefore, one first obvious choice is
to use
dµ = dρ and wi = 1, i = 1, . . . , n,
(8)
that is, sample according to the measure in which we plan to evaluate the L2 error and use equal weights.
When using equal weights wi = 1, the weighted least-squares estimator (2) becomes the standard leastsquares estimator, as a particular case. The strategy (8) was analyzed in [3], through the introduction of
the function
m
X
|Lj (x)|2 ,
x 7→ km (x) :=
j=1

which is the diagonal of the integral kernel of the projector Pm . This function only depends on Vm and dρ.
It is strictly positive in X due to Assumption 1. Its reciprocal function is characterized by
1
=
min
kvk2 ,
km (x) v∈Vm ,v(x)=1
and is called Christoffel function in the particular case where Vm is the space of algebraic polynomials of
total degree m − 1, see [10]. Obviously, the function km satisfies
Z
km dρ = m.
(9)
X

We define
Km = Km (Vm , dρ) := kkm kL∞ ,
and recall the following results from [3, 7] for the standard least-squares method with the weights and the
sampling measure chosen as in (8).
Theorem 1 For any r > 0, if m and n are such that the condition
Km ≤ κ

n
1 − ln 2
, with κ := κ(r) =
ln n
2 + 2r

(10)

is satisfied, then the following hold:
(i) The matrix G satisfies the tail bound


1
≤ 2n−r .
Pr kG − Ik2 >
2
4

(11)

(ii) If u ∈ L2 (X, dρ) satisfies a uniform bound (6), then the truncated least-squares estimator satisfies, in
the noiseless case,
E(ku − uT k2 ) ≤ (1 + ε(n))em (u)2 + 8τ 2 n−r ,
(12)
where ε(n) :=

4κ
ln(n)

→ 0 as n → +∞, and κ as in (10).

(iii) If u ∈ L∞ (X, dρ), then the truncated and nontruncated least-squares estimators satisfy, in the noiseless
case,
√
(13)
ku − uT k ≤ ku − uW k ≤ (1 + 2)em (u)∞ ,
with probability larger than 1 − 2n−r .
The second item in the above result shows that the optimal accuracy em (u) is met in expectation, up
to an additional term of order n−r . When em (u) has polynomial decay O(m−s ), we are ensured that this
additional term can be made negligible by taking r strictly larger than s/2, which amounts in taking κ(r)
small enough. Condition (10) imposes a minimal number of samples to ensure stability and accuracy of
standard least squares. Since (9) implies that Km ≥ m, the fulfillment of this condition requires that n is
at least of the order m ln(m). However simple examples show that the restriction can be more severe, for
example if Vm = Pm−1 on X = [−1, 1] and with ρ being the uniform probability measure. In this case, one
√
choice for the Lj are the Legendre polynomials with proper normalization kLj kL∞ = |Lj (1)| = 1 + 2j so
that Km = m2 , and therefore condition (10) imposes that n is at least of order m2 ln(m). Other examples in
the multivariate setting are discussed in [1, 2] which show that for many relevant approximation spaces Vm
and probability measures dρ, the behaviour of Km is superlinear in m, leading to a very demanding regime
in terms of the needed number n of samples. In the case of multivariate downward closed polynomial spaces,
precise upper bounds for Km have been proven in [2, 6] for measures associated to Jacobi polynomials.
In addition, note that the above theory does not cover simple situations such as algebraic polynomials
over unbounded domains, for example X = R equipped with the Gaussian measure, since the orthonormal
polynomials Lj are unbounded for j ≥ 2 and thus Km = ∞ if m ≥ 2.

2

Main results

In the present paper, we show that these limitations can be overcome, by using a proper weighted leastsquares method. We thus return to the general form of the discrete norm (4) used in the definition of the
weighted least-squares estimator. We now use a sampling measure dµ which generally differs from dρ and is
such that
wdµ = dρ,
R
where w is a positive function defined everywhere on X and such that X w−1 dρ = 1, and we then consider
the weighted least-square method with weights given by
wi = w(xi ).
With such a choice, the norm kvkn again approaches kvk as n increases. The particular case dµ = dρ and
w ≡ 1 corresponds to the standard least-squares method analyzed by Theorem 1. Note that changing the
sampling measure is a commonly used strategy for reducing the variance in Monte Carlo methods, where it
is referred to as importance sampling.

5

With Lj again denoting the L2 (X, dρ) orthonormal basis of Vm , we now introduce the function
x 7→ km,w (x) :=

m
X

w(x)|Lj (x)|2 ,

j=1

which only depends on Vm , dρ and w, as well as
Km,w = Km,w (Vm , dρ, w) := kkm,w kL∞ .

R
√
√
Note that, since the wLj are an L2 (X, dµ) orthonormal basis of wVm , we find that X km,w dµ = m and
thus Km,w ≥ m. We prove in this paper the following generalization of Theorem 1.
Theorem 2 For any r > 0, if m and n are such that the condition
1 − ln 2
n
, with κ :=
ln n
2 + 2r

(14)



1
≤ 2n−r .
Pr kG − Ik2 >
2

(15)

Km,w ≤ κ
is satisfied, then the following hold:
(i) The matrix G satisfies the tail bound

(ii) If u ∈ L2 (X, dρ) satisfies a uniform bound (6), then the truncated weighted least-squares estimator
satisfies, in the noiseless case,
E(ku − uT k2 ) ≤ (1 + ε(n))em (u)2 + 8τ 2 n−r ,
where ε(n) :=

4κ
ln(n)

(16)

→ 0 as n → +∞, and κ as in (10).

(iii) If u ∈ L∞ (X, dρ), then the nontruncated weighted least-squares estimators satisfy, in the noiseless case,
√
ku − uW k ≤ (1 + 2)em (u)∞ ,
(17)
with probability larger than 1 − 2n−r .
(iv) If u ∈ L2 (X, dρ), then the conditioned weighted least-squares estimator satisfies, in the noiseless case,
E(ku − uC k2 ) ≤ (1 + ε(n))em (u)2 + 2kuk2n−r ,
where ε(n) :=

4κ
ln(n)

(18)

→ 0 as n → +∞, and κ as in (10).

Let us mention that the quantity Km,w has been considered in [4], where similar stability and approximation results have been formulated in a slightly different form (see in particular Theorem 2.1 therein), in
the specific framework of total degree polynomial spaces.
The interest of Theorem 2 is that it leads us in a natural way to an optimal sampling strategy for the
weighted least-square method. We simply take
w :=

m
m
= Pm
,
2
km
j=1 |Lj |

(19)

and with such a choice for w one readily checks that

dµ :=

km
dρ,
m
6

(20)

R
is a probability measure on X since X km dρ = m.
In addition, we have for this particular choice that
km,w = wkm = m,
and therefore
Km,w = m.
We thus obtain the following result as a consequence of Theorem 2, which shows that the above choice of w
and dµ allows us to obtain near-optimal estimates for the truncated weighted least-squares estimator, under
the minimal condition that n is at least of the order m ln(m).
Corollary 1 For any r > 0, if m and n are such that the condition
m≤κ

n
1 − ln 2
, with κ :=
ln n
2 + 2r

(21)

is satisfied, then the conclusions (i), (ii), (iii) and (iv) of Theorem 2 hold for weighted least squares with the
choice of w and dµ given by (19) and (20).
One of the interests of the above optimal sampling strategy is that it applies to polynomial approximation
on unbounded domains that were not covered by Theorem 1, in particular X = R equipped with the Gaussian
measure. In this case, the relevant target functions u are often nonuniformly bounded and therefore the
results in items (ii) and (iii) of Theorem 2 do not apply. The result in item (iv) for the conditioned estimator
uC remains valid, since it does not require uniform boundedness of u.
Let us remark that all the above results are independent of the dimension d of the domain X. However,
raising d has the unavoidable effect of restricting the classes of functions for which the best approximation
error em (u) or em (u)∞ have some prescribed decay, due to the well-known curse of dimensionality.
Note that the optimal pair (dµ, w) described by (19) and (20) depends on Vm , that is
w = wm

and dµ = dµm .

This raises a difficulty for properly choosing the samples in settings where the choice of Vm is not fixed
a-priori, such as in adaptive methods. In certain particular cases, it is known that wm and dµm admit limits
w∗ and dµ∗ as m → ∞ and are globally equivalent to these limits. One typical example is given by the
univariate polynomial spaces Vm = Pm−1 , when X = [−1, 1] and dρ = ρdx where ρ is a Jacobi weight and
dx is the Lebesgue measure on X. In this case dµ∗ is the pluripotential equilibrium measure
dµ∗ =

dx
√
,
2π 1 − x2

see e.g. [11, 9], and one has
cdµ∗ ≤ dµm ≤ Cdµ∗ ,

m ≥ 1,

for some fixed constants 0 < c < C < ∞. Thus, in such a case, the above corollary also holds for the choice
w = w∗ and dµ = dµ∗ under the condition m ≤ Cc κ lnnn . The development of sampling strategies in cases of
varying values of m without such asymptotic equivalences is the object of current investigation.
A closely related weighted least-squares strategy was recently proposed and analyzed in [5], in the polynomial framework. There, the authors propose to use the renormalized Christoffel function (19) in the
definition of the weights, however sampling from the fixed pluripotential equilibrium measure dµ∗ . Due to
7

the fact that dµm differs from dµ∗ , the main estimate obtained in [5] (see p.3 therein) does not have the
same simple form of a direct comparison between ku − uT k and em (u) as in (ii) of Theorem 2. In particular,
it involves an extra term d(f ) which does not vanish even as n → ∞.
One intrinsic difficulty when using the optimal pair (dµ, w) = (dµm , wm ) described by (19) and (20) is the
effective sample generation, in particular in the multivariate framework since the measure dµm is generally
not of tensor product type. One possible approach is to use Markov Chain Monte Carlo methods such as
the Metropolis-Hastings algorithm, as explored in [4]. In such methods the samples are mutually correlated,
and only asymptotically distributed according to the desired sampling measure. One contribution of the
present paper is to propose a straightforward and effective sampling strategy for generating an arbitrary
finite number n of independent samples identically distributed according to dµm . This strategy requires that
dρ has tensor product structure and that the spaces Vm are spanned by tensor product bases, such as for
multivariate polynomial spaces, in which case dµm is generally not of tensor product type.
The rest of our paper is organized as follows. The proof of Theorem 2 is given in §3 in a concise form
since it follows the same lines as the original results on standard least squares from [3, 7]. We devote §4
to analog results in the case of samples affected by additive noise, proving that the estimates are robust
under condition (14). The proposed method for sampling the optimal measure dµm is discussed in §5, and
we illustrate its effectiveness in §6 by numerical examples.

3

Proof of Theorem 2

The proof is structurally similar to that of Theorem 1 given in [3] for items (i) and (ii) and in [2] for item
Pn
(iii), therefore we only sketch it. We observe that G = n1 i=1 Xi where the Xi are i.i.d. copies of the rank
1 random matrix


X = X(x) := w(x)Lj (x)Lk (x)
,
j,k=1,...,m

with x a random variable distributed over X according to µ. One obviously has E(X) = I. We then invoke
the Chernov bound from [12] to obtain that if kXk2 ≤ R almost surely, then, for any 0 < δ < 1,
Pr {kG − Ik2 > δ} ≤ 2m



 c 
1/R
e−δ
δ
=
2m
exp
−
,
(1 − δ)1−δ
R

(22)

with cδ := δ + (1 − δ) ln(1 − δ) > 0. Taking δ = 12 , and observing that
kX(x)k2 =

m
X
1
Km,w (x)
|Lj (x)|2 =
w(x)
,
n
n
j=1

K

which yields (15) in item (i).
we may thus take R = m,w
n
For the proof of (16) in item (ii), we first consider the event where kG − Ik2 ≤ 12 . In this case we write
n
n
ku − uT k2 = kTτ (u) − Tτ (uW )k2 ≤ ku − uW k2 = ku − Pm
uk2 ≤ kgk2 + kPm
gk2 ,
n
where we have used that Pm
Pm u = Pm u and that g is orthogonal to Vm , and thus

ku − uT k2 ≤ em (u)2 +
where a = (aj )j=1,...,m is solution of the system
Ga = b,
8

m
X
j=1

|aj |2 ,

g := u − Pm u,

and b := (hg, Lk in )k=1,...,m . Since kG−1 k2 ≤ 2, it follows that
ku − uT k2 ≤ em (u)2 + 4
In the event where kG − Ik2 >

1
2,

m
X

k=1

|hg, Lk in |2 .

we simply write ku − uT k ≤ 2τ . It follows that

E(ku − uT k2 ) ≤ em (u)2 + 4

m
X

k=1

E(|hg, Lk in |2 ) + 8τ 2 n−r .

For the second term, we have
n
n
1 XX
E(w(xi )w(xj )g(xi )g(xj )Lk (xi )Lk (xj ))
n2 i=1 j=1

1
= 2 n(n − 1)|E(w(x)g(x)Lk (x))|2 + nE(|w(x)g(x)Lk (x)|2 )
n
Z

1
1
|hg, Lk i|2 +
|w(x)|2 |g(x)|2 |Lk (x)|2 dµ
= 1−
n
n X
Z
1
=
w(x)|g(x)|2 |Lk (x)|2 dρ,
n X

E(|hg, Lk in |2 ) =

where we have used the fact that g is L2 (X, ρ)-orthogonal to Vm and thus to Lk . Summing over k, we obtain
m
X

k=1

E(|hg, Lk in |2 ) ≤

κ
Km,w
kgk2 ≤
em (u)2 ,
n
ln(n)

and we therefore obtain (16).
For the proof of (17) in item (iii) we place ourselves in the event where kG − Ik2 ≤ 12 . This property also
means that
3
1
kvk22 ≤ hGv, vi2 ≤ kvk22 , v ∈ Rm ,
2
2
which can be expressed as a norm equivalence over Vm ,
1
3
kvk2 ≤ kvk2n ≤ kvk2 ,
2
2

v ∈ Vm .

(23)

We then write that for any v ∈ Vm ,
n
ku − Pnm uk ≤ ku − vk + kv − Pm
uk
√
n
ukn
≤ ku − vk + 2kv − Pm
√
≤ ku − vk + 2ku − vkn
√
≤ (1 + 2)ku − vkL∞ ,

n
n
where we have used (23), the Pythagorean identity ku − vk2n = ku − Pm
uk2n + kv − Pm
uk2n , and the fact that
both k · k and k · kn are dominated by k · kL∞ . Since v is arbitrary, we obtain (17).
Finally, (18) in item (iv) is proven in a very similar way as (16) in item (ii), by writing that in the event
kG − Ik2 > 12 , we have ku − uC k = kuk, so that

E(ku − uC k2 ) ≤ em (u)2 + 4

m
X

k=1

and we conclude in the same way.

E(|hg, Lk in |2 ) + 2kuk2 n−r ,
✷

9

4

The noisy case

In a similar way as in [3, 8], we can analyze the case where the observations of u are affected by an additive
noise. In practical situations the noise may come from different sources, such as a discretization error when u
is evaluated by some numerical code, or a measurement error. The first one may be viewed as a perturbation
of u by a deterministic funtion h, that is, we observe
y i = u(xi ) + h(xi ).
The second one is typically modelled as a stochastic fluctuation, that is, we observe
y i = u(xi ) + η i .
where η i are independent realizations of the centered random variable η = y − u(x). Here, we do not
necessarily assume η and x to be independent, however we typically assume that the noise is centered, that
is,
E(η|x) = 0,
(24)
and we also assume uniformly bounded conditional variance
σ 2 := sup E(|η|2 |x) < ∞.

(25)

x∈X

Note that we may also consider consider a noncentered noise, which amounts in adding the two contributions,
that is,
y i = u(xi ) + β i , β i = h(xi ) + η i ,
(26)
with h(x) = E(β|x). The following result shows that the estimates in Theorem 2 are robust under the
presence of such an additive noise.
Theorem 3 For any r > 0, if m and n are such that condition (14) is satisfied, then the following hold for
the noise model (26):
(i) if u ∈ L2 (X, dρ) satisfies a uniform bound (6), then the truncated weighted least-squares estimator
satisfies
K m,w σ 2
+ 8τ 2 n−r ,
(27)
E(ku − uT k2 ) ≤ (1 + 2ε(n))em (u)2 + (8 + 2ε(n))khk2 +
n
(ii) if u ∈ L2 (X, dρ), then the conditioned weighted least-squares estimator satisfies
K m,w σ 2
+ 2kuk2 n−r ,
n
R
→ 0 as n → +∞, with κ as in (10), and K m,w := X km,w dρ.

E(ku − uC k2 ) ≤ (1 + 2ε(n))em (u)2 + (8 + 2ε(n))khk2 +
where in both cases ε(n) :=

4κ
ln(n)

(28)

Proof: We again first consider the event where kG − Ik2 ≤ 12 . In this case we write
ku − uT k ≤ ku − uW k,

n
and use the decomposition u − uW = g − Pm
g − h where g = u + Pm u as in the proof of Theorem 2 and h
stands for the solution to the least-squares problem for the noise data (β i )i=1,...,n . Therefore
n
n
n
ku − uW k2 = kgk2 + kPm
g + hk2 ≤ kgk2 + 2kPm
gk2 + 2khk2 = kgk2 + 2kPm
gk2 + 2

10

m
X
j=1

|nj |2 ,

where n = (nj )j=1,...,m is solution to
Gn = b,

b :=

n
1 X

n

i=1


β i w(xi )Lk (xi )

k=1,...,m

= (bk )k=1,...,m .

Since kG−1 k2 ≤ 2, it follows that
ku − uT k2 ≤ em (u)2 + 8

m
X

k=1

|hg, Lk in |2 + 8

m
X

k=1

|bk |2 .

Compared to the proof of Theorem 2, we need to estimate the expectation of the third term on the right
side. For this we simply write that
E(|bk |2 ) =

n
n
1 XX
E(β i w(xi )Lk (xi )β j w(xj )Lk (xj )).
n2 i=1 j=1

For i 6= j, we have
E(β i w(xi )Lk (xi )β j w(xj )Lk (xj )) = E(βw(x)Lk (x))2 = E(h(x)w(x)Lk (x))2 =

Z

2

hwLk dµ

X

= |hh, Lk i|2 .

Note that the first and second expectations are with respect to the joint density of (x, β) and the third one
with respect to the density of x, that is, µ. For i = j, we have
E(|β i w(xi )Lk (xi )|2 ) = E(|βw(x)Lk (x)|2 )
Z
E(|βw(x)Lk (x)|2 |x)dµ
=
X
Z
E(|β|2 |x)|w(x)Lk (x)|2 dµ
=
ZX
E(|β|2 |x)w(x)|Lk (x)|2 dρ
=
X
Z
(|h(x)|2 + E(|η|2 |x))w(x)|Lk (x)|2 dρ
=
X
Z
≤
(|h(x)|2 + σ 2 )w(x)|Lk (x)|2 dρ.
X

Summing up on i, j and k, and using condition (14), we obtain that
m
X

k=1


K m,w 2 
K m,w σ 2
1
Km,w
κ 
E(|bk |2 ) ≤ 1 − 2 khk2 +
khk2 +
khk2 +
σ ≤ 1+
.
n
n
n
log n
n

(29)

For the rest we proceed as for item (ii) and (iv) in the proof of Theorem 2, using that in the event kG−Ik2 > 12
we have ku − uT k ≤ 2τ and ku − uC k = kuk.
✷
Remark 1 Note that for the standard least-squares method, corresponding to the case where w ≡ 1, we know
2
that K m,w = m. The noise term thus takes the stardard form mσ
n , as seen for example in Theorem 3 2of [3]
κσ
or in Theorem 1 of [8]. Note that, in any case, condition (14) implies that this term is bounded by log
n.
The conclusions of Theorem 3 do not include the estimate in probability similar to item (iii) in Theorem
2. We can obtain such an estimate in the case of a bounded noise, where we assume that h ∈ L∞ (X) and
11

η is a bounded random variable, or equivalently, assuming that β is a bounded random variable, that is we
use the noise model (26) with
|β| ≤ D, a.s.
(30)
For this bounded noise model we have the following result.
Theorem 4 For any r > 0, if m and n are such that condition (14) is satisfied, then the following hold
for the the noise model (26) under (30): if u ∈ L∞ (X, dρ), then the nontruncated weighted least-squares
estimator satisfies
√
√
ku − uW k ≤ (1 + 2)em (u)∞ + 2D,
(31)
with probability larger than 1 − 2n−r .
Proof: Similar to the proof of (iii) in Theorem 2, we place ourselves in the event where kG − Ik2 ≤
use the norm equivalence (23). We then write that for any v ∈ Vm ,

1
2

and

n
n
ku − uW k ≤ ku − vk + kv − Pm
uk + kPm
βk.

The first two terms already appeared in the noiseless case and can be treated in the same way. The new
n
term Pm
β corresponds to the weighted least-squares approximation from the noise vector, and satisfies
√
√
√
n
n
βkn ≤ 2kβkn ≤ 2D.
kPm
βk ≤ 2kPm
This leads to (31).

5

✷

Random sampling from µm

The analysis in the previous sections prescribes the use of the optimal sampling measure dµm defined in
(20) for drawing the samples x1 , . . . , xn in the weighted least-squares method. In this section we discuss
numerical methods for generating independent random samples according to this measure, in a specific
relevant multivariate setting.
Here, we make the assumption that X = ×di=1 Xi is a Cartesian product of univariate real domains Xi ,
and that dρ is a product measure, that is,
d
O
dρi ,
dρ =
i=1

where each dρi is a measure defined on Xi . We assume that each dρi is of the form
dρi (t) = ρi (t)dt,

for some nonnegative continuous function ρi , and therefore
dρ(x) = ρ(x) dx,

ρ(x) =

d
Y

ρi (xi ),

i=1

x = (x1 , . . . , xd ) ∈ X.

In particular dρ is absolutely continuous with respect to the Lebesgue measure.
We consider the following general setting: for each i = 1, . . . , d, we choose a univariate basis (φij )j≥0
orthonormal in L2 (Xi , dρi ). We then define the tensorized basis
Lν (x) :=

d
Y

φiνi (xi ),

i=1

12

ν ∈ Nd0 ,

which is orthonormal in L2 (X, dρ). We consider general subspaces of the form
Vm := span{Lν : ν ∈ Λ},
for some multi-index set Λ ⊂ Nd0 such that #(Λ) = m. Thus we may rename the (Lν )ν∈Λ as (Lj )j=1,...,m
after a proper ordering has been chosen, for example in the lexicographical sense. For the given set Λ of
interest, we introduce
λj := max νj and λΛ := max λj .
ν∈Λ

j=1,...,d

The measure dµm is thus given by dµm (x) = µm (x)dx, where
µm (x) :=

m
1 X
1 X
|Li (x)|2 ρ(x) =
|Lν (x)|2 ρ(x),
m i=1
#(Λ)
ν∈Λ

x ∈ X.

(32)

We now discuss our sampling method for generating n independent random samples x1 , . . . , xn identically
distributed according to the multivariate density (32). Note that this density does not have a product
structure, despite ρ is a product density. There exist many methods for sampling from multivariate densities.
In contrast to Markov Chain Monte Carlo methods mentioned in the introduction, the method that we next
propose exploits the particular structure of the multivariate density (32), in order to generate independent
samples in a straightforward manner, and sampling only from univariate densities.
Given the vector x = (x1 , . . . , xd ) of all the coordinates, for any A ⊆ {1, . . . , d}, we introduce the notation
xA := (xi )i∈A ,

Ā := {1, . . . , d} \ A,

xĀ := (xi )i∈Ā ,

and
dxA :=

O

dxi ,

dρA :=

i∈A

O

dρi ,

ρA (xA ) :=

i∈A

Y

ρi (xi ),

i∈A

In the following, we mainly use the particular sets

XA := × Xi .
i∈A

Aq := {1, . . . , q} and Āq := {q + 1, . . . , d},
so that any x ∈ X may be written as x = (xAq , xĀq ).
Using such a notation, for any q = 1, . . . , d, we associate to the joint density µm its marginal density ψq
of the first q variables, namely
Z
ψq (xAq ) :=

XĀq

µm (xAq , xĀq ) dxĀq .

(33)

Since (φij )j≥0 is an orthonormal basis of L2 (Xi , dρi ), for any q = 1, . . . , d and any ν ∈ Nd0 , we obtain that
Z

XĀq

2

|Lν (xAq , xĀq )| ρ(xAq , xĀq )dxĀq = ρAq (xAq )

q
Y

i=1

|φiνi (xi )|2 ,

xAq ∈ XAq .

Therefore, the marginal density (33) can be written in simple form as
q

ψq (xAq ) =

XY
1
|φiνi (xi )|2 .
ρAq (xAq )
#(Λ)
i=1
ν∈Λ

13

(34)

Sequential conditional sampling. Based on the previous notation and remarks, we propose an algorithm
which generates n samples xk = (xk1 , . . . , xkd ) ∈ X with k = 1, . . . , n, that are independent and identically
distributed realizations from the density µm in (32).
In the multivariate case the coordinates can be arbitrarily reordered. Start with the first coordinate x1
and sample n points x11 , . . . , xn1 ∈ X1 from the univariate density
ϕ1 : X1 → R : t 7→ ϕ1 (t) := ψ1 (t) =

ρ1 (t) X 1
|φν1 (t)|2 ,
#(Λ)

(35)

ν∈Λ

which coincides with the marginal ψ1 of x1 calculated in (34). In the univariate case d = 1 the algorithm
terminates. In the multivariate case d ≥ 2, by iterating q from 2 to d, consider the qth coordinate xq ,
and sample n points x1q , . . . , xnq ∈ Xq in the following way: for any k = 1, . . . , n, given the values xkAq−1 =
(xk1 , . . . , xkq−1 ) ∈ XAq−1 that have been calculated at the previous q − 1 steps, sample the point xkq ∈ Xq from
the univariate density
Qq−1 j k 2
P
2
q
j=1 |φνj (xj )|
ν∈Λ |φνq (t)|
k
ϕq : Xq → R : t 7→ ϕq (t|xAq−1 ) := ρq (t)
.
(36)
P
Qq−1 j k 2
ν∈Λ
j=1 |φνj (xj )|
The expression on the right-hand side of (36) is continuous at any t ∈ Xq and at any xkAq−1 ∈ XAq−1 .
Assumption 1 ensures that the denominator of (36) is strictly positive for any possible choice of xkAq−1 =
(xk1 , . . . , xkq−1 ) ∈ XAq−1 , and also ensures that the marginal ψq−1 is strictly positive at any point xkAq−1 ∈
XAq−1 such that ρAq−1 (xkAq−1 ) 6= 0. For any t ∈ Xq and any xkAq−1 ∈ XAq−1 such that ρAq−1 (xkAq−1 ) 6= 0, the
density ϕq satisfies
ϕq (t|xkAq−1 ) =

ψq (xkAq−1 , t)
,
ψq−1 (xkAq−1 )

(37)

where the densities ψq and ψq−1 are the marginals defined in (33) and evaluated at the points (xkAq−1 , t) ∈
X q and xkAq−1 ∈ XAq−1 , respectively. From (37), using (34) and simplifying the term ρAq−1 (xkAq−1 ) =
QAq−1
k
j=1 ρj (xj ) 6= 0, one obtains the right-hand side of (36). The right-hand side of equation (37) is well
defined for any t ∈ Xq and any xkAq−1 ∈ XAq−1 such that ρAq−1 (xkAq−1 ) 6= 0, and it is not defined at the
points xkAq−1 ∈ XAq−1 such that ρAq−1 (xkAq−1 ) = 0 where ψq−1 (xkAq−1 ) vanishes. Nonetheless, (37) has finite
limits at any point (xkAq−1 , t) ∈ XAq , and these limits equal expression (36).
According to technical terminology, the right-hand side of equation (37) is the conditional density of
xq given x1 , . . . , xq−1 with respect to the density ψq , and ϕq is the continuous extension to XAq of this
conditional density.
The densities ϕ1 , . . . , ϕd defined in (35)–(36) can be concisely rewritten for any q = 1, . . . , d as
X
(38)
αν (xkAq−1 )|φqνq (t)|2 ,
ϕq (t|xkAq−1 ) = ρq (t)
ν∈Λ

where the nonnegative weights (αν )ν∈Λ are defined as

1


if q = 1,

 #(Λ) ,
Q
q−1
αν = αν (zAq−1 ) :=
|φjνj (zj )|2

 P j=1
, if 2 ≤ q ≤ d,
Q

q−1
j

2
ν∈Λ
j=1 |φνj (zj )|
P
for any zAq−1 = (z1 , . . . , zq−1 ) ∈ XAq−1 . Since ν∈Λ αν = 1, each density ϕq in (38) is a convex combination
of the densities ρq |φq1 |2 , . . . , ρq |φqλq |2 . Note that if the orthonormal basis (φqj )j≥0 have explicit expressions
14

and can be evaluated at any point in Xq , then the same holds for the univariate densities (38). In particular,
in the polynomial case, for standards univariate densities ρi such as uniform, Chebyshev or Gaussian, the
orthonormal polynomials (φij )j≥1 have expressions which are explicitely computable, for example by recursion
formulas.
In Algorithm 1 we summarize our sampling method, that sequentially samples the univariate densities
(38) to generate independent samples from the multivariate density (32). In the univariate case d = 1 the
algorithm does not run the innermost loop, and only samples from ϕ1 . In the multivariate case d ≥ 2
the algorithm runs also the innermost loop, and conditionally samples also from ϕ2 , . . . , ϕd . Our algorithm
therefore relies on accurate sampling methods for the relevant univariate densities (38).

Algorithm 1 Sequential conditional sampling for µm .
INPUT: n, d, Λ, ρi , (φij )j≥0 for i = 1, . . . , d.
i.i.d.

OUTPUT: x1 , . . . , xn ∼ µm .
for k = 1 to n do
αν ← (#(Λ))−1 , for any ν ∈ Λ.

Sample xk1 from t 7→ ϕ1 (t) = ρ1 (t)
for q = 2 to d do
q−1
Q j k 2
|φνj (xj )|
αν ←

j=1

P q−1
Q

ν∈Λ j=1

|φjνj (xkj )|2

P

ν∈Λ

αν |φ1ν1 (t)|2 .

, for any ν ∈ Λ.

Sample xkq from t 7→ ϕq (t) = ρq (t)
end for
xk ← (xk1 , . . . , xkd ).
end for

P

ν∈Λ

αν |φqνq (t)|2 .

We close this section by discussing two possible methods for sampling from such densities: rejection
sampling and inversion transform sampling. Both methods equally apply to any univariate density ϕq , and
therefore we present them for any q arbitrarily chosen from 1 to d.
Rejection sampling (RS). For applying this method, one needs to find a suitable univariate density Θq ,
whose support contains the support of ϕq , and a suitable real Mq > 1 such that
ϕq (t) ≤ Mq Θq (t),

t ∈ supp(ϕq ).

The density Θq should be easier to sample than ϕq , i.e. efficient pseudorandom number generators for
sampling from Θq are available. The value of Mq should be the smallest possible. For sampling one point
from ϕq using RS: sample a point z from Θq , and sample u from the standard uniform U(0, 1). Then check if
u < ϕq (z)/Mq Θq (z): if this is the case then accept z as a realization from ϕq , otherwise reject z and restart
sampling z and u from beginning. On average, acceptance occurs once every Mq trials. Therefore, for a
given q, sampling one point from ϕq by RS requires on average Mq evaluations of the function
t 7→

ρq (t) X
ϕq (t)
=
αν |φqνq (t)|2 .
Mq Θq (t)
Mq Θq (t)
ν∈Λ

15

This amounts in evaluating Mq times the terms φq0 , φqλq and a subset of the terms φq1 , . . . , φqλq −1 , depending
on Λ. The coefficients αν depend on the terms φj0 , . . . , φjλj for j = 1, . . . , q − 1, which have been already
evaluated when sampling the previous coordinates 1, . . . , q −1. Thus, if we use RS for sampling the univariate
densities, the overall computational cost of Algorithm 1 for sampling n points x1 , . . . , xn ∈ X is on average
Pd
proportional to n q=1 Mq (λq + 1).
When the basis functions (φqj )j≥0 form a bounded orthonormal system, an immediate and simple choice
of the parameters in the algorithm is
Mq = max kφqνq k2L∞ ,

and

ν∈Λ

Θq (t) = ρq (t).

(39)

With such a choice, we can quantify more precisely the average computational cost of sampling n points
√
in dimension d. When (φqj )j≥0 are the Chebyshev polynomials, whose L∞ norms satisfy kφqj kL∞ ≤ 2, we
Pd
obtain the bound 2n q=1 (λq + 1) ≤ 2nd(λΛ + 1) ≤ 2ndm. When (φqj )j≥0 are the Legendre polynomials,
Pd
√
whose L∞ norms satisfy kφqj kL∞ ≤ 2j + 1, we have the crude estimate 2n q=1 (λq + 1)2 ≤ 2nd(λΛ + 1)2 ≤
2ndm2 . In general, when (φqj )j≥0 are Jacobi polynomials, similar upper bounds can be derived, and the
dependence of these bounds on n and d is linear.
Inversion transform sampling (ITS). Let Φq : Xq → [0, 1] be the cumulative distribution function
associated to the univariate density ϕq . In the following, only when using the ITS method, we make
the further assumption that ρq vanishes at most a finite number of times in Xq . Such an assumption is
fulfilled in many relevant situations, e.g. when ρq is the density associated to Jacobi or Hermite polynomials
orthonormal in L2 (Xq , dρq ). Together with Assumption 1, this ensures that the function t 7→ Φq (t) is
continuous and strictly increasing on Xq . Hence Φq is a bijection between Xq and [0, 1], and it has a
unique inverse Φ−1
: [0, 1] → Xq which is continuous and strictly increasing on [0, 1]. Sampling from ϕq
q
using ITS can therefore be performed as follows: sample n independent realizations u1 , . . . , un identically
distributed according to the standard uniform U(0, 1), and obtain the n independent samples from ϕq as
1
−1 n
(Φ−1
q (u ), . . . , Φq (u )).
For any u ∈ [0, 1], computing z = Φ−1
q (u) ∈ Xq is equivalent to find the unique solution z ∈ Xq to
Φq (z) = u. This can be executed by elementary root-finding numerical methods, e.g. the bisection method
or Newton’s method. In alternative to using root-finding methods, one can build an interpolant operator Iq
−1
of Φ−1
q , and then approximate Φq (u) ≈ Iq (u) for any u ∈ [0, 1]. Such an interpolant Iq can be constructed
for example by piecewise linear interpolation, from the data (Φq (tq1 ), tq1 ), . . . , (Φq (tqsq ), tqsq ) at sq suitable
points tq1 < . . . < tqsq in Xq .
Both root-finding methods and the interpolation method require evaluating the function Φq pointwise in
Xq . In general these evaluations can be computed using standard univariate quadrature formulas. When
(φqj )j≥0 are orthogonal polynomials, the explicit expression of the primitive of ϕq can be used for directly
evaluating the function Φq .
Finally we discuss the overall computational cost of Algorithm 1 for sampling n points x1 , . . . , xn ∈ X
when using ITS for sampling the univariate densities. With the bisection method, this overall cost amounts
Pd
to n q=1 γq Wq , where γq is the maximum number of iterations for locating the zero in Xq up to some desired
tolerance, and Wq is the computational cost of each iteration. With the interpolation of Φ−1
q , the overall
cost amounts to n evaluations of each interpolant Iq , in addition to the cost for building the interpolants
which does not depend on n.

16

6

Examples and numerical illustrations

This section presents the numerical performances of the weighted least-squares method compared to the
standard least-squares method, in three relevant situations where dρ can be either the uniform measure,
the Chebyshev measure, or the Gaussian measure. In each one of these three cases, we choose w and dµ
in the weighted least-squares method from (19) and (20), as prescribed by our analysis in Corollary 1. For
standard least squares we choose w and dµ as in (8). Our tests focus on the condition number of the Gramian
matrix, that quantifies the stability of the linear system (5) and the stability of the weighted and standard
least-squares estimators. A meaningful quantity is therefore the probability
Pr{cond(G) ≤ 3},

(40)

where, through (7), the value three of the threshold is related to the parameter δ = 1/2 in the previous
analysis. For any n and m, from (7) the probability (40) is larger than Pr{kG − Ik2 ≤ 21 }. From Corollary 1,
under condition (21) between n, m and r, the Gramian matrix of weighted least squares satisfies (15) and
therefore the probability (40) is larger than 1 − 2n−r . For standard least squares, from Theorem 1 the
Gramian matrix satisfies (40) with probability larger than 1 − 2n−r , but under condition (10).
In the numerical tests the probability (40) is approximated by empirical probability, obtained by counting
how many times the event cond(G) ≤ 3 occurs when repeating the random sampling one hundred times.
All the examples presented in this section confine to multivariate approximation spaces of polynomial
type. One natural assumption in this case is to require that the set Λ is downward closed, that is, satisfies
ν∈Λ

and ν̃ ≤ ν =⇒ ν̃ ∈ Λ,

where ν̃ ≤ ν means that ν̃j ≤ νj for all i = 1, . . . , d. Then Vm is the polynomial space spanned by the
monomials
d
Y
ν
zj j ,
z 7→ z ν :=
j=1

and the orthonormal basis Lν is provided by taking each (φij )j≥0 to be a sequence of univariate orthonormal
polynomials of L2 (Xi , dρi ).
In both the univariate and multivariate forthcoming examples, the random samples from the measure
dµm are generated using Algorithm 1. The univariate densities ϕ1 , . . . , ϕd are sampled using the inversion
transform sampling method. The inverse of the cumulative distribution function is approximated using the
interpolation technique.

6.1

Univariate examples

In the univariate case d = 1, let the index set be Λ = {0, . . . , m − 1} and Vm = PΛ = span{z k : k =
0, . . . , m − 1}. We report in Fig. 1 the probability (40), approximated by empirical probability, when G is
the Gramian matrix of the weighted least-squares method. Different combinations of values for m and n are
tested, with three choices of the measure dρ: uniform, Gaussian and Chebyshev. The results do not show
perceivable differences among the performances of weighted least squares with the three different measures.
In any of the three cases, n/ ln(n) ≥ 4m is enough to obtain an empirical probability equal to one that
cond(G) ≤ 3. This confirms that condition (21) with any choice of r > 0 ensures (40), since it demands for
a larger number of samples.
Fig. 2 shows the probability (40) when G is the Gramian matrix of standard least squares. With the
uniform measure, the condition n/ ln(n) ≥ m2 is enough to have (40) with empirical probability larger
17

dρ uniform measure

dρ Gaussian measure

dρ Chebyshev measure

Figure 1: Weighted least squares, P r{cond(G) ≤ 3}, d = 1. Left: dρ uniform measure. Center: dρ Gaussian
measure. Right: dρ Chebyshev measure.
dρ uniform measure

dρ Gaussian measure

dρ Chebyshev measure

Figure 2: Standard least squares, P r{cond(G) ≤ 3}, d = 1. Left: dρ uniform measure. Center: dρ Gaussian
measure. Right: dρ Chebyshev measure.
than 0.95. When dρ is the Gaussian measure, stability requires a very large number of evaluations, roughly
n/ ln(n) linearly proportional to exp(m/3). For the univariate Chebyshev measure, it is proven that standard
least squares are stable under the same minimal condition (21) as for weighted least squares. In accordance
with the theory, the numerical results obtained in this case with weighted and standard least squares are
indistinguishable, see Fig. 1-right and Fig. 2-right.

6.2

Multivariate examples

Afterwards we present some numerical tests in the multivariate setting. These tests are again based, as
in the previous section, on approximating the probability (40) by empirical probability. In dimension d
larger than one there are many possible ways to enrich the polynomial space PΛ . The number of different
downward closed sets whose cardinality equals m gets very large already for moderate values of m and d.
Therefore, we present the numerical results for a chosen sequence of polynomial spaces PΛ1 , . . . , PΛm such
that Λ1 ⊂ · · · ⊂ Λm , where each Λj ⊂ Nd0 is downward closed, #(Λj ) = dim(PΛj ) = j and the starting set
Λ1 contains only the null multi-index. All the tests in Fig. 3 and Fig. 4 have been obtained using the same
sequence of increasingly embedded polynomial spaces PΛ1 ⊂ . . . ⊂ PΛm , for both weighted and standard least
18

squares and for the three choices of the measures dρ. Such a choice allows us to establish a fair comparison
between the two methods and among different measures, without the additional variability arising from
modifications to the polynomial space.
dρ uniform measure

dρ Gaussian measure

dρ Chebyshev measure

Figure 3: Weighted least squares, P r{cond(G) ≤ 3}, d = 10. Left: dρ uniform measure. Center: dρ Gaussian
measure. Right: dρ Chebyshev measure.
dρ uniform measure

dρ Gaussian measure

dρ Chebyshev measure

Figure 4: Standard least squares, P r{cond(G) ≤ 3}, d = 10. Left: dρ uniform measure. Center: dρ Gaussian
measure. Right: dρ Chebyshev measure.
We report the results obtained for the tests in dimension d = 10. The results in Fig. 3 confirm that
weighted least squares always yield an empirical probability equal to one that cond(G) ≤ 3, provided that
n/ log(n) ≥ 2m. This condition ensures that (21) with any choice of r > 0 implies (40), thus verifying
Corollary 1. Again, the results do not show significant differences among the three choices of the measure
dρ: a straight line, with the same slope for all the three cases uniform, Chebyshev and Gaussian, separates
the two regimes corresponding to empirical probabilities equal to zero and one. Compared to the univariate
case in Fig. 1, the results in Fig. 3 exhibit a sharper transition between the two extreme regimes, and an
overall lower variability in the transition regime.
The results for standard least squares with d = 10 are shown in Fig. 4. In the case of the uniform
measure, in Fig. 4-right, stability is ensured if n/ ln(n) ≥ 3.5m, which is more demanding than the condition
n/ ln(n) ≥ 2m needed for the stability of weighted least squares in Fig. 3-right, but much less strict than
the condition required with standard least squares in the univariate case, where n/ ln(n) scales like m2 .
19

These phenomena have already been observed and described in [7]. Similar results as those with the uniform
measure are obtained with the Chebyshev measure in Fig. 4-left, where again standard least squares achieve
stability using more evaluations than weighted least squares in Fig. 3-left. The case of the Gaussian measure
drastically differs from the uniform and Chebyshev cases: the results in Fig. 4-center clearly indicate that a
very large number of evaluations n compared to m is required to achieve stability of standard least squares.
Let us mention that analogous results as those presented in Figs. 1 and 3 for weighted least squares have
been obtained also in other dimensions, and with many other sequences of increasingly embedded polynomial
spaces. In the next tables we report some of these results for selected values of d = 1, 2, 5, 10, 50, 100. We
choose n = 26599 and m = 200 that satisfy condition (21) with r = 1, and report in Table 1 the empirical
probabilities that approximate (40), again calculated over one hundred repetitions. This table provides
multiple comparisons: weighted least squares versus standard least squares, for the three choices of the
measure dρ (uniform, Gaussian and Chebyshev) and with d varying between 1 and 100.
method

dρ

d=1

d=2

d=5

d = 10

d = 50

d = 100

weighted LS
weighted LS
weighted LS

uniform
Gaussian
Chebyshev

1
1
1

1
1
1

1
1
1

1
1
1

1
1
1

1
1
1

standard LS
standard LS
standard LS

uniform
Gaussian
Chebyshev

0
0
1

0
0
1

0.54
0
1

1
0
1

1
0
1

1
0
1

Table 1: P r{cond(G) ≤ 3}, with n = 26559 and m = 200: weighted least squares versus standard least
squares, dρ uniform versus dρ Gaussian versus dρ Chebyshev, d = 1, 2, 5, 10, 50, 100.
method

dρ

d=1

d=2

d=5

d = 10

d = 50

d = 100

weighted LS
weighted LS
weighted LS

uniform
Gaussian
Chebyshev

1.5593
1.5994
1.5364

1.4989
1.5698
1.4894

1.4407
1.4743
1.4694

1.4320
1.4643
1.4105

1.4535
1.4676
1.4143

1.4179
1.4237
1.4216

standard LS
standard LS
standard LS

uniform
Gaussian
Chebyshev

19.9584
∼ 1019
1.5574

29.8920
∼ 1019
1.5367

3.0847
∼ 1019
1.5357

1.9555
∼ 1016
1.4752

1.7228
∼ 109
1.4499

1.5862
∼ 103
1.4625

Table 2: Average of cond(G), with n = 26559 and m = 200: weighted least squares versus standard least
squares, dρ uniform versus dρ Gaussian versus dρ Chebyshev, d = 1, 2, 5, 10, 50, 100.
In Table 1, all the empirical probabilities related to results for weighted least squares are equal to one,
and confirm the theory since, for the chosen values of n, m and r, the probability (40) is larger than
1 − 5.67 × 10−7 . This value is computed using estimate (22) from the proof of Theorem 2. In contrast
to weighted least squares, whose empirical probability equal one independently of dρ and d, the empirical
probability of standard least squares does depend on the chosen measure, and to some extent on the dimension
d as well. With the uniform measure, the empirical probability that approximates (40) equals zero when
d = 1 or d = 2, equals 0.54 when d = 5, and equals one when d = 10, d = 50 or d = 100. In the Gaussian
case, standard least squares always feature null empirical probabilities. With the Chebyshev measure, the
condition number of G for standard least squares is always lower than three for any tested value of d.
In addition to the results in Table 1, further information are needed for assessing how severe is the lack
20

of stability when obtaining null empirical probabilities. To this aim, in Table 2 we also report the average
value of cond(G), obtained when averaging the condition number of G over the same repetitions used to
estimate the empirical probabilities in Table 1. The information in Table 2 are complementary to those in
Table 1. On the one hand they point out the stability and robustness of weighted least squares, showing a
tamed condition number with any measure dρ and any dimension d. On the other hand they provide further
insights on stability issues of standard least squares and their dependence on dρ and d. For standard least
squares with the uniform measure, the average condition number reduces as the dimension d increases, in
agreement with the conclusion drawn from Table 1. The Gramian matrix of standard least squares with
the Gaussian measure is very ill-conditioned for all tested values of d. For standard least squares with the
Chebyshev measure, the averaged condition number of G is only slightly larger than the one for weighted
least squares.
It is worth remarking that, the results for standard least squares in Fig. 4, Table 1 and Table 2 are
sensitive to the chosen sequence of polynomial spaces. Testing different sequences might produce different
results, that however necessarily obey to the estimates proven in Theorem 1 with uniform and Chebyshev
measures, when n, m and r satisfy condition (10). Many other examples with standard least squares have
been extensively discussed in previous works e.g. [7, 2], also in situations where n, m and r do not satisfy
condition (10) and therefore Theorem 1 does not apply. In general, when n, m and r do not satisfy (10), there
exist multivariate polynomial spaces of dimension m such that the Gramian matrix of standard least squares
with the uniform and Chebyshev measures does not satisfy (11). Examples of such spaces are discussed in
[7, 2]. Using these spaces would yield null empirical probabilities in Table 1 for standard least squares with
the uniform and Chebyshev measures.
For weighted least squares, when n, m and r satisfy condition (21), any sequence of polynomial spaces
yields empirical probabilities close to one, according to Corollary 1. Indeed such a robustness with respect
to the choices of dρ, of the polynomial space and of the dimension d represents one of the main advantages
of the weighted approach.

References
[1] G. Chardon, A. Cohen, and L. Daudet, Sampling and reconstruction of solutions to the Helmholtz
equation, Sampl. Theory Signal Image Process., 13:67–89, 2014.
[2] A. Chkifa, A. Cohen, G. Migliorati, F. Nobile, and R. Tempone, Discrete least squares polynomial
approximation with random evaluations - application to parametric and stochastic elliptic PDEs, M2AN,
49(3):815–837, 2015.
[3] A. Cohen , M.A. Davenport, and D. Leviatan, On the stability and accuracy of least squares approximations, Found. Comput. Math., 13:819–834, 2013.
[4] A. Doostan and J. Hampton, Coherence motivated sampling and convergence analysis of least squares
polynomial Chaos regression, Comput. Methods Appl. Mech. Engrg., 290:73–97,2015.
[5] J.D. Jakeman, A. Narayan, and T. Zhou, A Christoffel function weighted least squares algorithm for
collocation approximations, preprint.
[6] G. Migliorati, Multivariate Markov-type and Nikolskii-type inequalities for polynomials associated with
downward closed multi-index sets, J. Approx. Theory, 189:137–159, 2015.

21

[7] G. Migliorati, F. Nobile, E. von Schwerin, and R. Tempone, Analysis of discrete L2 projection on
polynomial spaces with random evaluations, Found. Comput. Math., 14:419–456, 2014.
[8] G. Migliorati, F. Nobile, and R. Tempone, Convergence estimates in probability and in expectation for
discrete least squares with noisy evaluations at random points, J. Multivar. Analysis, 142:167–182, 2015.
[9] E.B. Saff, and V. Totik, Logarithmic Potentials with External Fields, Springer, 1997.
[10] P. Nevai, Géza Freud, orthogonal polynomials and Christoffel Functions. A case study, J. Approx. theory,
48:3–167, 1986.
[11] A. Máté, P. Nevai, and V. Totik, Szegö’s extremum problem on the unit circle, Annals of Mathematics,
134:433–453, 1991.
[12] J. Tropp, User friendly tail bounds for sums of random matrices, Found. Comput. Math., 12:389–434,
2012.

22

