Polynomial approximation of noisy functions

Approximating a univariate function on the interval $$[-1,1]$$ with a polynomial is among the most classical problems in numerical analysis. When the function evaluations come with noise, a least-squares fit is known to reduce the effect of noise as more samples are taken. The generic algorithm for the least-squares problem requires $$O(Nn^2)$$ operations, where $$N+1$$ is the number of sample points and n is the degree of the polynomial approximant. This algorithm is unstable when n is large, for example $$n\gg \sqrt{N}$$ for equispaced sample points. In this study, we blend numerical analysis and statistics to introduce a stable and fast $$O(N\log N)$$ algorithm called NoisyChebtrunc based on Chebyshev interpolation. It has the same error reduction effect as least-squares and the convergence is spectral until the error reaches $$O(\sigma \sqrt{{n}/{N}})$$ , where $$\sigma$$ is the noise level, after which the error continues to decrease at the Monte-Carlo $$O(1/\sqrt{N})$$ rate. To determine the polynomial degree, NoisyChebtrunc employs a statistical criterion, namely Mallows’ $$C_p$$ . We analyze NoisyChebtrunc in terms of the variance and concentration in the infinity norm to the underlying noiseless function. These results show that with high probability the infinity-norm error is bounded by a small constant times $$\sigma \sqrt{{n}/{N}}$$ , when the noise is independent and follows a subgaussian or subexponential distribution. We illustrate the performance of NoisyChebtrunc with numerical experiments.