How do least-squares, Chebyshev, and minimax approximation differ?

They differ in the norm used to measure the error. Least-squares approximation minimizes the squared-error norm ∫(f-p)², is robust to outliers, and only requires solving the normal equations. Minimax approximation minimizes the maximum absolute error max|f-p| (the ∞ norm), spreading the error as evenly as possible across the interval, but the computation is iterative. Chebyshev approximation sits in between: its coefficients can be computed directly from the Chebyshev polynomials, and because its error is nearly equioscillating it serves as a 'near-minimax' and a practical substitute for true minimax.

What does the coefficient of determination R² represent?

R² is a measure between 0 and 1 of how much of the variance of the response variable the regression model explains. It is defined as R² = 1 - (residual sum of squares SS_res / total sum of squares SS_tot); the closer to 1, the better the fit. sangi's linearRegression returns R² = Sxy² / (Sxx·Syy) for simple regression, and weightedLinearRegression returns 1 - SS_res/SS_tot from the weighted residual sum of squares.

Why is an approximation using Chebyshev polynomials nearly minimax?

The Chebyshev polynomials T_n(cosθ)=cos nθ oscillate equally between ±1 on the interval [-1,1]. When a Chebyshev series is truncated to finitely many terms, the leading error term is roughly proportional to this T_n, so the error oscillates with nearly equal amplitude across the whole interval. Equioscillation is the signature of the best uniform (minimax) approximation, so Chebyshev approximation achieves an error very close to true minimax (typically within a few percent).

Why is high-degree polynomial fitting dangerous?

When you least-squares fit a high-degree polynomial in the power basis 1, x, x², ... at equally spaced sample points, the condition number of the normal-equation coefficient matrix (the Gram matrix derived from a Vandermonde matrix) deteriorates exponentially with degree, so the coefficients lose significant digits and become meaningless. sangi's polynomialFit subtracts the sample mean to improve conditioning, but the fundamental remedy is to use an orthogonal basis such as the Chebyshev basis.

Approximation: Least Squares, Chebyshev & Minimax

Overview

"Approximation" means replacing a given function $f$ or a sequence of data points with a simpler, more tractable function $p$ (usually a polynomial). The character of the problem changes with the criterion by which we measure how close $f$ and $p$ are. sangi's approximation.hpp offers methods that address the following three criteria.

Squared norm (least squares): minimizes the sum of squared errors $\sum (f-p)^2$ or the integral $\int (f-p)^2\,dx$. The computation only solves a linear system (the normal equations) once, and because it averages out the influence of outliers and noise it is well suited to data fitting.
$\infty$ norm (minimax): minimizes the maximum absolute error $\max_x |f(x)-p(x)|$. Because it spreads the error as evenly as possible across the interval, it is well suited to numerical-library implementations that need to evaluate a function to a uniform accuracy, but the computation is iterative.
Chebyshev (near-minimax): expands the function in a series of Chebyshev polynomials of the first kind. The coefficients can be computed directly from orthogonality, and because the truncation error is nearly equioscillating it serves as a good substitute for minimax (near-minimax).

Linear Regression and Polynomial Fitting (Least Squares)

Normal equations for simple regression

We fit a line $y = a x + b$ to a sequence of data points $(x_i, y_i)$. Partially differentiating the residual sum of squares $\sum_i (y_i - a x_i - b)^2$ with respect to $a, b$ and setting it to $0$ yields a $2 \times 2$ system of normal equations. sangi's linearRegression obtains the answer in closed form from the means $\bar{x}, \bar{y}$ and the centered cross-products $S_{xx}, S_{xy}, S_{yy}$:

$$a = \frac{S_{xy}}{S_{xx}}, \qquad b = \bar{y} - a\,\bar{x}, \qquad S_{xy} = \sum_i (x_i - \bar{x})(y_i - \bar{y})$$

Normal equations for polynomial fitting

To least-squares fit a polynomial $p(x) = \sum_{j=0}^{d} a_j x^j$ of degree $d$, set the design matrix $X$ to $X_{ij} = x_i^{\,j}$. The coefficient vector $\beta = (a_0,\dots,a_d)^\top$ that minimizes the residual sum of squares satisfies the normal equations:

$$X^\top X\,\beta = X^\top y$$

sangi's polynomialFit builds this Gram matrix $X^\top X$ and its right-hand side $X^\top y$ from the power sums $\sum_k x_k^{\,p}$, and solves them by LU decomposition (algorithms::solve).

Coefficient of determination $R^2$

The quality of fit is measured by the coefficient of determination. From the total sum of squares $\mathrm{SS}_{\text{tot}}$ and the residual sum of squares $\mathrm{SS}_{\text{res}}$:

$$R^2 = 1 - \frac{\mathrm{SS}_{\text{res}}}{\mathrm{SS}_{\text{tot}}} = 1 - \frac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i (y_i - \bar{y})^2}$$

$R^2 = 1$ means a perfect fit, and $0$ means it is no better than predicting the mean. For simple regression this is equivalent to $R^2 = S_{xy}^2 / (S_{xx} S_{yy})$, and it is stored in the r_squared field of LinearRegressionResult.

Weighted regression

When the reliability of each sample differs, we minimize $\sum_i w_i (y_i - a x_i - b)^2$, attaching a weight $w_i$ to each residual. weightedLinearRegression assembles the following weighted normal equations and solves the $2\times 2$ system:

$$\begin{bmatrix} \sum w_i x_i^2 & \sum w_i x_i \\ \sum w_i x_i & \sum w_i \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} \sum w_i x_i y_i \\ \sum w_i y_i \end{bmatrix}$$

Taking the inverse of the variance $w_i = 1/\sigma_i^2$ as the weight gives the best linear unbiased estimate (generalized least squares) in the Gauss-Markov sense.

Beware Vandermonde ill-conditioning

The design matrix $X$ is essentially a Vandermonde matrix, and as the degree rises the condition number of the Gram matrix $X^\top X$ deteriorates rapidly. Using the power basis directly with equally spaced samples can make the high-degree coefficients lose significance and become meaningless.

sangi's polynomialFit solves in a space centered by subtracting the sample means $\mu_x, \mu_y$, then maps back to the original coefficients via the binomial expansion, which somewhat improves the conditioning. Even so, there are limits at high degree (roughly $d \gtrsim 6$), so for high degrees the right approach is to switch to the Chebyshev basis of the next section.

Chebyshev Approximation

Chebyshev polynomials of the first kind

The Chebyshev polynomials of the first kind $T_n$ are defined by the substitution $x = \cos\theta$ as:

$$T_n(\cos\theta) = \cos(n\theta), \qquad n = 0, 1, 2, \dots$$

From the trigonometric addition formulas, a three-term recurrence is obtained:

$$T_0(x) = 1, \quad T_1(x) = x, \quad T_{n+1}(x) = 2x\,T_n(x) - T_{n-1}(x)$$

On the interval $[-1, 1]$, $T_n$ oscillates $n$ times between values of $\pm 1$, and is orthogonal with respect to the weight $1/\sqrt{1-x^2}$. This orthogonality is the key to computing the coefficients.

Chebyshev coefficients (computation via discrete orthogonality)

When a function $f$ is expressed as a Chebyshev series $f(x) \approx \tfrac{1}{2} c_0 + \sum_{j\ge 1} c_j T_j(x)$, the coefficients $c_j$ can be written as integrals from orthogonality. In the actual computation we use a discrete orthogonality sum, sampling $f$ at $N$ Chebyshev nodes (Chebyshev-Gauss points of the first kind):

$$x_k = \cos\!\left(\frac{\pi (k + \tfrac{1}{2})}{N}\right), \qquad c_j = \frac{2}{N} \sum_{k=0}^{N-1} f(x_k)\,\cos\!\left(\frac{\pi j (k + \tfrac{1}{2})}{N}\right)$$

This is exactly a discrete cosine transform (DCT-II); note that the half-integer shift $k+\tfrac{1}{2}$ is used so that the nodes do not land on the zeros of $f$. ChebyshevApprox computes and stores these $c_j$ in its constructor.

Linear transformation to the interval $[a, b]$

Since the domain of the Chebyshev polynomials is $[-1, 1]$, a general interval $[a, b]$ is mapped onto it by a linear transformation:

$$y = \frac{2x - a - b}{b - a} \;\in [-1, 1], \qquad x = \frac{(b-a)\,y + (a+b)}{2}$$

ChebyshevApprox(f, a, b, n) performs this transformation internally, so the user can pass $x$ in the original interval $[a, b]$ directly.

Evaluation (Clenshaw recurrence) and truncation

When computing the approximation $\sum_j c_j T_j(y)$, expanding $T_j$ explicitly is unstable. ChebyshevApprox evaluates stably using the Clenshaw recurrence. Starting from $d_{m} = d_{m+1} = 0$, for $j = m-1, \dots, 1$:

$$d_j = 2y\,d_{j+1} - d_{j+2} + c_j, \qquad p(x) = y\,d_1 - d_2 + \tfrac{1}{2} c_0$$

The coefficients $c_j$ shrink rapidly with degree (exponentially if $f$ is smooth), so even truncating (economization) at a suitable number of terms $m$ barely loses accuracy. evaluate(x, m) can evaluate using only the leading $m$ terms, trading accuracy against cost as the application requires.

If you want to convert the coefficients back to an ordinary polynomial $\sum_j a_j x^j$, toPolynomialCoefficients() expands each $T_k$ into the power basis using the three-term recurrence, composes the interval transformation $y = \alpha x + \beta$ via the binomial expansion, and returns the coefficient sequence.

Why "near-minimax"?

When a smooth $f$ is truncated to a degree-$n$ Chebyshev series, the leading error term is roughly proportional to $c_{n+1} T_{n+1}(x)$. Since $T_{n+1}$ oscillates equally between $\pm 1$ across the whole interval, the approximation error $f - p$ also oscillates with nearly equal amplitude. As we will see, equioscillation is precisely the signature of the best uniform (minimax) approximation, so Chebyshev approximation achieves an error very close to true minimax (typically within a few percent in maximum error). To distinguish this from true minimax we call it "near-minimax".

Related article: Chebyshev approximation

Minimax Approximation

Minimizing the uniform norm

When approximating a continuous function $f$ in the space $\mathcal{P}_n$ of polynomials of degree $n$, the polynomial $p^*$ that minimizes the maximum absolute error (the $\infty$ norm, the uniform norm) is called the best uniform approximation (minimax approximation):

$$p^* = \arg\min_{p \in \mathcal{P}_n} \;\max_{x \in [a,b]} |f(x) - p(x)|$$

If $f$ is continuous, the best approximation $p^*$ exists and is moreover unique (Chebyshev's theorem).

Chebyshev's equioscillation theorem

What characterizes the best approximation is the equioscillation theorem. A necessary and sufficient condition for $p \in \mathcal{P}_n$ to be the best uniform approximation of $f$ is that the error function $e(x) = f(x) - p(x)$ attains the maximum absolute error $E = \max_x |e(x)|$ with alternating sign at at least $n+2$ points $x_0 < x_1 < \dots < x_{n+1}$ within the interval:

$$e(x_k) = (-1)^k\,\sigma\,E, \qquad \sigma = \pm 1, \qquad k = 0, 1, \dots, n+1$$

That is, $n+2$ points line up where the error swings to its extreme $E$ with alternating sign. As seen in the previous section, the error of Chebyshev approximation is close to this equioscillating form, which is why it is a good approximation to minimax.

Remez exchange algorithm (concept)

The classical iterative method for finding the true minimax polynomial is the Remez exchange algorithm. Conceptually it repeats the following:

Tentatively choose $n+2$ reference points (the reference set). Chebyshev nodes make good initial points.
At those points, solve the equioscillation condition $e(x_k) = (-1)^k E$ as a linear system to find the coefficients and the equal amplitude $E$.
Find the extremal points of the error $e(x)$ of the resulting $p$, and "exchange" the reference set toward those extremal points.
Repeat until the reference set stops moving (equioscillation is achieved).

The Remez method converges quadratically to the best approximation, but each iteration requires an extremum search and the implementation is heavy.

sangi's ChebyshevApprox provides a near-minimax approximation via Chebyshev series expansion, not the Remez exchange algorithm. Its coefficients can be computed in one shot by a DCT, and because its error is nearly equioscillating it is a practical substitute for true minimax in many applications. When strict equioscillation is required, it is standard practice to use the Chebyshev-approximation coefficients as the initial guess for the Remez iteration.

Comparison Table

Criterion	Error minimized	Method	Complexity	Main use
Least squares	Squared norm $\sum (f-p)^2$	Solve the normal equations $X^\top X\beta = X^\top y$ once	$O(n^2 m + n^3)$ ($m$ samples, degree $n$)	Fitting/regression of noisy data
Chebyshev (near-minimax)	Nearly the $\infty$ norm (close to equioscillation)	Compute coefficients directly by DCT of node samples	$O(N^2)$ (naive DCT)	High-accuracy approximation of smooth functions, coefficient analysis
Minimax	$\infty$ norm $\max\|f-p\|$	Remez exchange algorithm (iterative)	Iterations × (extremum search + linear system)	Function-evaluation routines that must guarantee a uniform accuracy

A rough guide: to fit noisy observed data, use least squares (linearRegression / polynomialFit). To approximate a known smooth function to high accuracy across the whole interval and also work with its coefficients, use Chebyshev (ChebyshevApprox). Only when you must strictly minimize the maximum error should you reach for minimax (Remez), seeding it with a Chebyshev approximation as the initial guess — that is the practical workflow.

References

Trefethen, L. N. (2013). Approximation Theory and Approximation Practice. SIAM.
Cheney, E. W. (1966). Introduction to Approximation Theory. McGraw-Hill.
Powell, M. J. D. (1981). Approximation Theory and Methods. Cambridge University Press.
Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. (2007). Numerical Recipes: The Art of Scientific Computing (3rd ed.). Cambridge University Press. (Chebyshev approximation and the Clenshaw recurrence)
Björck, Å. (1996). Numerical Methods for Least Squares Problems. SIAM.