  Weighted least squares polynomial approximation uses random samples to
determine projections of functions onto spaces of polynomials. It has been
shown that, using an optimal distribution of sample locations, the number of
samples required to achieve quasi-optimal approximation in a given polynomial
subspace scales, up to a logarithmic factor, linearly in the dimension of this
space. However, in many applications, the computation of samples includes a
numerical discretization error. Thus, obtaining polynomial approximations with
a single level method can become prohibitively expensive, as it requires a
sufficiently large number of samples, each computed with a sufficiently small
discretization error. As a solution to this problem, we propose a multilevel
method that utilizes samples computed with different accuracies and is able to
match the accuracy of single-level approximations with reduced computational
cost. We derive complexity bounds under certain assumptions about polynomial
approximability and sample work. Furthermore, we propose an adaptive algorithm
for situations where such assumptions cannot be verified a priori. Finally, we
provide an efficient algorithm for the sampling from optimal distributions and
an analysis of computationally favorable alternative distributions. Numerical
experiments underscore the practical applicability of our method.
