Least-Squares Estimation for Polynomial Terms

Why Least-Squares Estimation Needs Equations

A polynomial regression model represents the mean response as a polynomial in a predictor. Once the degree and the predictor values are chosen, estimating the coefficients means choosing a curve that fits the observed responses as closely as possible under a specified measure of discrepancy. Least squares uses the sum of squared residuals: it squares each observed-minus-fitted difference and adds the results.

The polynomial model has the same linear-in-the-coefficients structure as a linear regression model, even when it curves as a function of the predictor. This structure lets us express the estimation problem with a design matrix and derive equations that characterize its solution. The equations are useful both for understanding what a least-squares fit does and for checking a computed fit.

Model and assumptions: For fixed observed predictor values \(x_1,\ldots,x_n\), a degree-\(d\) polynomial regression model is \(Y_i=\beta_0+\beta_1x_i+\cdots+\beta_dx_i^d+\varepsilon_i\). For the usual statistical interpretation, assume the errors have conditional mean zero, are independent, and have a common variance. Normal errors are not needed to define or find the least-squares estimates. In fact, minimizing the sum of squared residuals is an algebraic procedure and can be carried out without these distributional assumptions.

For observed responses \(y_1,\ldots,y_n\), the fitted value at \(x_i\) is \(\hat y_i=\hat\beta_0+\hat\beta_1x_i+\cdots+\hat\beta_dx_i^d\). The residual is \(e_i=y_i-\hat y_i\). The estimates are chosen to minimize \(S(\hat\beta_0,\ldots,\hat\beta_d)=\sum_{i=1}^n e_i^2\). The hats distinguish estimates calculated from the sample from the model coefficients \(\beta_0,\ldots,\beta_d\).

The Polynomial Design Matrix and Normal Equations

The design matrix records the polynomial terms for each observed predictor value. Its first column contains ones for the intercept, its next column contains \(x_i\), the next contains \(x_i^2\), and so on through \(x_i^d\). With \(\mathbf y\) denoting the vector of observed responses, the fitted values are \(\mathbf X\hat{\boldsymbol\beta}\), and the residual vector is \(\mathbf e=\mathbf y-\mathbf X\hat{\boldsymbol\beta}\).

\mathbf X= \begin{bmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^d\\ 1 & x_2 & x_2^2 & \cdots & x_2^d\\ \vdots & \vdots & \vdots & & \vdots\\ 1 & x_n & x_n^2 & \cdots & x_n^d \end{bmatrix}, \qquad \mathbf y= \begin{bmatrix}y_1\\y_2\\\vdots\\y_n\end{bmatrix}, \qquad \hat{\boldsymbol\beta}= \begin{bmatrix}\hat\beta_0\\\hat\beta_1\\\vdots\\\hat\beta_d\end{bmatrix}.

In this notation, the sum of squared residuals is \(S(\hat{\boldsymbol\beta})=(\mathbf y-\mathbf X\hat{\boldsymbol\beta})^\mathsf{T}(\mathbf y-\mathbf X\hat{\boldsymbol\beta})\). A least-squares estimate satisfies the normal equations:

\mathbf X^\mathsf{T}\mathbf X\hat{\boldsymbol\beta} = \mathbf X^\mathsf{T}\mathbf y.

For polynomial terms, the same equations can be written one equation per coefficient. For \(j=0,\ldots,d\), multiply each residual by \(x_i^j\), sum over observations, and set the result to zero:

\sum_{i=1}^n x_i^j \left(y_i-\sum_{k=0}^d\hat\beta_kx_i^k\right)=0, \qquad j=0,\ldots,d.

Equivalently, the \(j\)th equation is \(\sum_i x_i^j y_i=\sum_{k=0}^d\hat\beta_k\sum_i x_i^{j+k}\). These are called normal equations because they require the residual vector to be orthogonal to every column of the design matrix. In ordinary language, the residuals have zero summed product with each included polynomial term.

For a quadratic model, the three equations are especially explicit. The residuals sum to zero, their products with \(x_i\) sum to zero, and their products with \(x_i^2\) sum to zero:

\begin{aligned} \sum_i y_i&=n\hat\beta_0+\hat\beta_1\sum_i x_i+\hat\beta_2\sum_i x_i^2,\\ \sum_i x_i y_i&=\hat\beta_0\sum_i x_i+\hat\beta_1\sum_i x_i^2+\hat\beta_2\sum_i x_i^3,\\ \sum_i x_i^2 y_i&=\hat\beta_0\sum_i x_i^2+\hat\beta_1\sum_i x_i^3+\hat\beta_2\sum_i x_i^4. \end{aligned}

Worked Example: Forming the Quadratic Normal Equations

Suppose three hypothetical observations have predictor values \(-1,0,1\) and responses \(1,2,5\). Fit a quadratic polynomial. The predictor values are distinct, so the three columns \(1,x,x^2\) are linearly independent for these observations; the coefficients are uniquely determined.

The design matrix and response vector are

\mathbf X= \begin{bmatrix} 1&-1&1\\ 1&0&0\\ 1&1&1 \end{bmatrix}, \qquad \mathbf y= \begin{bmatrix}1\\2\\5\end{bmatrix}.

The normal equations use \(\mathbf X^\mathsf{T}\mathbf X\) and \(\mathbf X^\mathsf{T}\mathbf y\). Multiplication gives

\begin{bmatrix} 3&0&2\\ 0&2&0\\ 2&0&2 \end{bmatrix} \begin{bmatrix}\hat\beta_0\\\hat\beta_1\\\hat\beta_2\end{bmatrix} = \begin{bmatrix}8\\4\\6\end{bmatrix}.

The second equation gives \(2\hat\beta_1=4\), so \(\hat\beta_1=2\). The first and third equations are \(3\hat\beta_0+2\hat\beta_2=8\) and \(2\hat\beta_0+2\hat\beta_2=6\). Subtracting the third from the first gives \(\hat\beta_0=2\); then \(4+2\hat\beta_2=6\), so \(\hat\beta_2=1\). The fitted polynomial is \(\hat y=2+2x+x^2\). It passes through all three observed points, so its residuals and sum of squared residuals are all zero.

Why the Normal Equations Give a Minimum

The normal equations are not merely a convenient recipe. They characterize a minimum of the sum of squared residuals. Here is a direct justification. Suppose \(\hat{\boldsymbol\beta}\) satisfies the normal equations, and consider any other coefficient vector \(\hat{\boldsymbol\beta}+\boldsymbol\delta\). Its residual vector is \(\mathbf e-\mathbf X\boldsymbol\delta\). Expanding the squared length gives

\begin{aligned} S(\hat{\boldsymbol\beta}+\boldsymbol\delta) &=(\mathbf e-\mathbf X\boldsymbol\delta)^\mathsf{T} (\mathbf e-\mathbf X\boldsymbol\delta)\\ &=\mathbf e^\mathsf{T}\mathbf e -2\boldsymbol\delta^\mathsf{T}\mathbf X^\mathsf{T}\mathbf e +\boldsymbol\delta^\mathsf{T}\mathbf X^\mathsf{T}\mathbf X\boldsymbol\delta. \end{aligned}

The normal equations imply \(\mathbf X^\mathsf{T}\mathbf e=\mathbf X^\mathsf{T}(\mathbf y-\mathbf X\hat{\boldsymbol\beta})=\mathbf 0\), so the middle term vanishes. The last term is a squared length:

\boldsymbol\delta^\mathsf{T}\mathbf X^\mathsf{T}\mathbf X\boldsymbol\delta =(\mathbf X\boldsymbol\delta)^\mathsf{T}(\mathbf X\boldsymbol\delta) \geq 0.

Thus \(S(\hat{\boldsymbol\beta}+\boldsymbol\delta)\geq S(\hat{\boldsymbol\beta})\) for every alternative coefficient change \(\boldsymbol\delta\). A solution of the normal equations is therefore a least-squares minimizer. If the columns of \(\mathbf X\) are linearly independent, then \(\mathbf X\boldsymbol\delta\neq\mathbf 0\) whenever \(\boldsymbol\delta\neq\mathbf 0\), making the inequality strict for every nonzero change. The minimizing coefficient vector is then unique.

For a degree-\(d\) polynomial, having at least \(d+1\) distinct predictor values ensures that the columns \(1,x,\ldots,x^d\) are linearly independent. If this condition fails, or if the number and pattern of observations do not provide enough independent information, coefficients may not be unique. A least-squares fit still exists, but different coefficient vectors can sometimes give the same fitted values.

Worked Example: Solving an Overdetermined Quadratic Fit

Consider five hypothetical observations at \(x=-2,-1,0,1,2\), with responses \(y=4,1,3,5,12\), respectively. There are more observations than coefficients, and there are five distinct predictor values, so the quadratic design matrix has full column rank. The least-squares equations are

\begin{aligned} 5\hat\beta_0+0\hat\beta_1+10\hat\beta_2&=25,\\ 0\hat\beta_0+10\hat\beta_1+0\hat\beta_2&=20,\\ 10\hat\beta_0+0\hat\beta_1+34\hat\beta_2&=70. \end{aligned}

For example, the right sides are \(\sum_i y_i=25\), \(\sum_i x_i y_i=20\), and \(\sum_i x_i^2y_i=70\). The second equation gives \(\hat\beta_1=2\). From the first, \(\hat\beta_0=5-2\hat\beta_2\). Substitution in the third gives \(10(5-2\hat\beta_2)+34\hat\beta_2=70\), so \(14\hat\beta_2=20\), \(\hat\beta_2=10/7\), and \(\hat\beta_0=15/7\). The fitted polynomial is

\hat y=\frac{15}{7}+2x+\frac{10}{7}x^2.

At the five observed predictor values, the residuals are \(1/7,-4/7,6/7,-4/7,1/7\). They are not all zero, but they satisfy the required orthogonality checks: \(\sum_i e_i=0\), \(\sum_i x_i e_i=0\), and \(\sum_i x_i^2e_i=0\). The sum of squared residuals is \[ SSE=\frac{1^2+(-4)^2+6^2+(-4)^2+1^2}{49} =\frac{70}{49}=\frac{10}{7}. \] The nonzero residuals reflect that a single quadratic does not pass through every point; the normal equations show that no other quadratic has a smaller sum of squared residuals.

Worked Example: A Linear Polynomial and Its Residual Checks

A linear polynomial is the degree-one case. For hypothetical data \((x,y)=(0,1),(1,2),(2,2)\), the sums are \(n=3\), \(\sum x_i=3\), \(\sum x_i^2=5\), \(\sum y_i=5\), and \(\sum x_i y_i=6\). The degree-one normal equations are

3\hat\beta_0+3\hat\beta_1=5, \qquad 3\hat\beta_0+5\hat\beta_1=6.

Subtract the first equation from the second: \(2\hat\beta_1=1\), giving \(\hat\beta_1=0.5\). Then \(3\hat\beta_0+1.5=5\), so \(\hat\beta_0=7/6\approx1.167\). The fitted line is \(\hat y=7/6+0.5x\). The residuals are \(-1/6,1/3,-1/6\), so their sum is \(-1/6+1/3-1/6=0\). Their predictor-weighted sum is \(0(-1/6)+1(1/3)+2(-1/6)=0\). These are precisely the two normal-equation checks, one for the intercept column and one for the \(x\) column. The sum of squared residuals is \(1/36+1/9+1/36=1/6\).

Reading and Using a Polynomial Least-Squares Fit

The normal equations explain two important features of a polynomial least-squares fit. First, each included column imposes an orthogonality condition on the residuals. With an intercept, residuals sum to zero; with a linear term, their products with \(x_i\) sum to zero; with a quadratic term, their products with \(x_i^2\) sum to zero. These statements concern sums over the observed data, not each residual individually.

Second, least squares minimizes the squared discrepancies in the response units squared. It does not guarantee that every observation is close to the fitted curve, nor does it show that the polynomial relationship is appropriate for all predictor values. The estimates describe the polynomial that minimizes the chosen objective for the observed data. Interpreting the model beyond that estimation task requires attention to the data, model assumptions, and context.

In practical computation, software typically solves the least-squares problem using numerical methods rather than explicitly calculating \((\mathbf X^\mathsf{T}\mathbf X)^{-1}\mathbf X^\mathsf{T}\mathbf y\). The inverse expression is valid when the matrix has full column rank, but forming it directly can be numerically less stable. The normal equations remain valuable as the mathematical characterization of the solution and as a way to check residual behavior.

Common Mistake

Common mistake: “Least squares makes every residual zero, so the fitted polynomial must pass through every observation.” This is true only in special cases, such as the three-point quadratic example above. In an overdetermined fit, residuals are generally nonzero. The correct condition is that residuals are orthogonal to each included design-matrix column, and that the total squared residual is no larger than for any other coefficient vector.

Key takeaway: Polynomial least squares builds columns \(1,x,x^2,\ldots,x^d\), then chooses coefficients to minimize the sum of squared residuals. The normal equations \(\mathbf X^\mathsf{T}\mathbf X\hat{\boldsymbol\beta}=\mathbf X^\mathsf{T}\mathbf y\) express the resulting residual orthogonality, and their solution is a global minimum because any coefficient change adds a nonnegative squared length.

Check Your Understanding

Answer each question using the normal-equation and minimization ideas developed above.

For a cubic polynomial, what four columns appear in the design matrix, and how many normal equations are there?
What does the intercept-column normal equation say about the residuals when the model includes an intercept?
Why does \(\mathbf X^\mathsf{T}\mathbf X\hat{\boldsymbol\beta}=\mathbf X^\mathsf{T}\mathbf y\) imply that no alternative coefficient vector has a smaller sum of squared residuals?
For a quadratic fit, what additional residual check is imposed by the \(x^2\) column?
What condition on the predictor values guarantees that the coefficients of a degree-\(d\) polynomial are uniquely determined?