Fast and accurate noise removal by curve fitting using orthogonal polynomials

Local polynomial smoothing is a widespread technique in data analysis, and Savitzky-Golay (SG) filters are one of its most well-known realizations. In real settings, the effectiveness of SG filtering depends critically on proper tuning of its parameters, constrained in turn by repeated polynomial fitting over large data windows and for varying polynomial degrees. Standard implementations based on monomial bases and Vandermonde matrix formulations are known to suffer from ill-conditioning and unfavorable scaling as the problem size increases. In this work, we present a fast and numerically stable method for computing polynomial fitting and differentiation matrices by reformulating the problem in terms of discrete orthogonal (Chebyshev) polynomials. Exploiting their recursive structure and the intrinsic symmetry properties of the resulting matrices, we derive two algorithms designed to reduce computational overhead. Both methods significantly reduce memory usage and improve scalability with respect to the polynomial degree and window length. A discussion of the performance demonstrates that the proposed algorithms achieve orders-of-magnitude improvements in numerical accuracy compared to standard matrix multiplication, while also providing potential gains in execution time for large-scale problems. These features make the approach particularly well suited for applications requiring repeated local polynomial fits, such as the optimization of SG filters in high-resolution spectral analyses, including axion dark matter searches and the ALPHA haloscope.