Polynomial Fits as Least-Squares Problems
Polynomial regression uses a polynomial to summarize measured data. Given input values \(x_1,\dots,x_m\) and responses \(y_1,\dots,y_m\), a degree-\(d\) model has the form \(p(x)=c_0+c_1x+\cdots+c_dx^d\). The coefficients enter linearly, even though the model uses powers of the input. This makes fitting the polynomial a least-squares problem in a design matrix, and reduced QR gives a way to solve it without forming \(\mathbf{A}^{\top}\mathbf{A}\).
Throughout, the data vector is \(\mathbf{b}=(y_1,\dots,y_m)^{\top}\in\mathbb{R}^m\). The fitted values are the entries of \(\mathbf{A}\mathbf{c}\), and the residual vector is \(\mathbf{r}=\mathbf{b}-\mathbf{A}\mathbf{c}\). The residual sum of squares is \(\|\mathbf{r}\|^2\), the sum of the squares of its entries.
The order of the columns fixes the order of the coefficients: the first column corresponds to the constant term, the next to \(x\), and so on. A reduced QR factorization of a full-column-rank matrix \(\mathbf{A}\in\mathbb{R}^{m\times n}\) writes \(\mathbf{A}=\mathbf{Q}\mathbf{R}\), where \(\mathbf{Q}\in\mathbb{R}^{m\times n}\) has orthonormal columns and \(\mathbf{R}\in\mathbb{R}^{n\times n}\) is upper triangular with nonzero diagonal entries. In the polynomial setting, \(n=d+1\).
When the Polynomial Coefficients Are Unique
For a degree-\(d\) model, if there are at least \(d+1\) distinct input values, then the design matrix has independent columns. Indeed, if \(\mathbf{A}\mathbf{c}=\mathbf{0}\), the polynomial \(p(x)=c_0+c_1x+\cdots+c_dx^d\) vanishes at every input value. A nonzero polynomial of degree at most \(d\) has at most \(d\) distinct roots, so \(p\) must be the zero polynomial and all its coefficients are zero. Thus the columns are independent. In particular, the full-column-rank condition required for the reduced QR solve holds for these data.
For a full-column-rank design matrix, the QR method obtains the coefficient vector by solving \(\mathbf{R}\mathbf{c}=\mathbf{Q}^{\top}\mathbf{b}\). Because \(\mathbf{R}\) is upper triangular, this can be done by back substitution. The following result explains why this gives the least-squares fit, rather than merely a solution of a convenient equation.
Proof. Since \(\mathbf{Q}^{\top}\mathbf{Q}=\mathbf{I}_n\), for any \(\mathbf{z}\in\mathbb{R}^n\) we can write \(\mathbf{b}-\mathbf{A}\mathbf{z}=(\mathbf{b}-\mathbf{Q}\mathbf{Q}^{\top}\mathbf{b})+\mathbf{Q}(\mathbf{Q}^{\top}\mathbf{b}-\mathbf{R}\mathbf{z})\). The first summand is orthogonal to every column of \(\mathbf{Q}\), because \(\mathbf{Q}^{\top}(\mathbf{b}-\mathbf{Q}\mathbf{Q}^{\top}\mathbf{b})=\mathbf{Q}^{\top}\mathbf{b}-\mathbf{Q}^{\top}\mathbf{Q}\mathbf{Q}^{\top}\mathbf{b}=\mathbf{0}\). The second summand is a linear combination of those columns. By the Pythagorean identity, the squared residual therefore equals \(\|\mathbf{b}-\mathbf{Q}\mathbf{Q}^{\top}\mathbf{b}\|^2+\|\mathbf{Q}(\mathbf{Q}^{\top}\mathbf{b}-\mathbf{R}\mathbf{z})\|^2\). Orthonormal columns preserve the norm of coefficient vectors, so the second term is \(\|\mathbf{Q}^{\top}\mathbf{b}-\mathbf{R}\mathbf{z}\|^2\). It is nonnegative and is zero exactly when \(\mathbf{R}\mathbf{z}=\mathbf{Q}^{\top}\mathbf{b}\). Since \(\mathbf{R}\) is invertible, this equation has exactly one solution, which is the unique minimizer. Substituting that solution gives \(\mathbf{A}\mathbf{c}=\mathbf{Q}\mathbf{R}\mathbf{c}=\mathbf{Q}\mathbf{Q}^{\top}\mathbf{b}\); the first term in the decomposition is its residual, giving the stated residual sum of squares. \(\square\)
A Quadratic Fit by Reduced QR
Worked Example: A Quadratic Model with a Nonzero Residual
Fit a degree-two polynomial to the data with inputs \(-2,-1,1,2\) and responses \(2,-2,6,8\). The design matrix uses columns for \(1,x,x^2\), in that order:
Take the following reduced QR factors:
The columns of \(\mathbf{Q}\) have norm \(1\): their squared norms are respectively \(4(1/4)=1\), \((4+1+1+4)/10=1\), and \(4(1/4)=1\). Their pairwise dot products are \(0\): for columns 1 and 2 the numerator is \(-2-1+1+2=0\); for columns 1 and 3 the sum is \(1/2-1/2-1/2+1/2=0\); for columns 2 and 3 the numerator is \(-2+1-1+2=0\). To check the factorization, its columns multiply as follows:
These are exactly the columns of \(\mathbf{A}\), so \(\mathbf{A}=\mathbf{Q}\mathbf{R}\). Next compute the right-hand side of the triangular system:
Back substitution in \(\mathbf{R}\mathbf{c}=\mathbf{Q}^{\top}\mathbf{b}\) gives \(3c_2=3\), hence \(c_2=1\); then \(\sqrt{10}c_1=2\sqrt{10}\), hence \(c_1=2\); and \(2c_0+5c_2=7\), hence \(2c_0+5=7\) and \(c_0=1\). Thus the fitted polynomial is \(p(x)=1+2x+x^2\). Multiplication back checks the coefficients and fitted values:
The residual is \(\mathbf{r}=\mathbf{b}-\mathbf{A}\mathbf{c}=(1,-2,2,-1)^{\top}\). Therefore the residual sum of squares is \(1^2+(-2)^2+2^2+(-1)^2=10\). As a further check on the fit, the residual is orthogonal to each design column: its dot products with \((1,1,1,1)^{\top}\), \((-2,-1,1,2)^{\top}\), and \((4,1,1,4)^{\top}\) are respectively \(1-2+2-1=0\), \(-2+2+2-2=0\), and \(4-2+2-4=0\).
Two More Polynomial Fits
Worked Example: A Linear Polynomial Fit
Fit a degree-one polynomial to the inputs \(-1,0,1\) and responses \(0,1,3\). Here the design matrix has columns for \(1,x\):
A reduced QR factorization is \(\mathbf{Q}=\begin{bmatrix}1/\sqrt{3}&-1/\sqrt{2}\\1/\sqrt{3}&0\\1/\sqrt{3}&1/\sqrt{2}\end{bmatrix}\) and \(\mathbf{R}=\begin{bmatrix}\sqrt{3}&0\\0&\sqrt{2}\end{bmatrix}\). The columns of \(\mathbf{Q}\) have squared norms \(3(1/3)=1\) and \(1/2+1/2=1\), and their dot product is \((-1+0+1)/\sqrt{6}=0\). Multiplying \(\mathbf{Q}\mathbf{R}\) gives columns \((1,1,1)^{\top}\) and \((-1,0,1)^{\top}\), so the factorization reproduces \(\mathbf{A}\).
Now \(\mathbf{Q}^{\top}\mathbf{b}=(4/\sqrt{3},3/\sqrt{2})^{\top}\). The triangular equations give \(\sqrt{2}c_1=3/\sqrt{2}\), so \(c_1=3/2\), and \(\sqrt{3}c_0=4/\sqrt{3}\), so \(c_0=4/3\). The fitted polynomial is \(p(x)=4/3+(3/2)x\), and its fitted values are \(\mathbf{A}\mathbf{c}=(-1/6,4/3,17/6)^{\top}\). Subtracting from the data gives \(\mathbf{r}=(1/6,-1/3,1/6)^{\top}\), so the residual sum of squares is \(1/36+1/9+1/36=1/6\). The residual is orthogonal to both design columns: \(1/6-1/3+1/6=0\), and \((-1)(1/6)+1/6=0\).
Worked Example: A Different Response on the Same Quadratic Design
Keep the inputs \(-2,-1,1,2\) and the quadratic design matrix from the first example, but use \(\mathbf{b}=(1,0,0,0)^{\top}\). Since the design matrix is unchanged, use the same verified \(\mathbf{Q}\) and \(\mathbf{R}\). The transformed response is
Solving \(\mathbf{R}\mathbf{c}=\mathbf{Q}^{\top}\mathbf{b}\), the last equation gives \(3c_2=1/2\), so \(c_2=1/6\). The middle equation gives \(\sqrt{10}c_1=-2/\sqrt{10}\), so \(c_1=-1/5\). The first gives \(2c_0+5/6=1/2\), so \(c_0=-1/6\). Thus \(p(x)=-1/6-(1/5)x+(1/6)x^2\). Evaluating at the four inputs gives \(\mathbf{A}\mathbf{c}=(9/10,1/5,-1/5,1/10)^{\top}\). The residual is consequently \((1/10,-1/5,1/5,-1/10)^{\top}\), and its squared norm is \(1/100+1/25+1/25+1/100=1/10\). Its dot products with the three columns of the design matrix are \(0\), \(0\), and \(0\), respectively, confirming residual orthogonality.
Why the Triangular Solve Matters
The QR computation has three distinct stages: build the polynomial design matrix, transform the data to \(\mathbf{Q}^{\top}\mathbf{b}\), then solve the triangular system \(\mathbf{R}\mathbf{c}=\mathbf{Q}^{\top}\mathbf{b}\). The transformed data are not themselves the polynomial coefficients; the triangular solve converts them to coefficients in the original design columns. Once \(\mathbf{c}\) is known, \(\mathbf{A}\mathbf{c}\) gives fitted values, \(\mathbf{b}-\mathbf{A}\mathbf{c}\) gives residuals, and summing the squared residual entries gives the residual sum of squares.
One advantage of this route is that it does not require forming \(\mathbf{A}^{\top}\mathbf{A}\). Instead, the method uses the orthonormal columns of \(\mathbf{Q}\) and a triangular system in \(\mathbf{R}\). Full column rank matters: if the polynomial design matrix does not have independent columns, the invertible triangular solve in the theorem is not available as stated, and the least-squares coefficients need not be unique.
Common Mistake: Treating \(Q^{\top}b\) as the Coefficients
In the first quadratic example, a tempting but incorrect step would be to set the coefficient vector equal to \(\mathbf{Q}^{\top}\mathbf{b}=(7,2\sqrt{10},3)^{\top}\). These are the correct transformed-response entries, but they are not the coefficients of \(1,x,x^2\). Substitution into the required triangular equation exposes the error:
The correction is to solve \(\mathbf{R}\mathbf{c}=\mathbf{Q}^{\top}\mathbf{b}\), which gives \(\mathbf{c}=(1,2,1)^{\top}\) in that example. A second setup error is to put the data values \(y_i\) into the design matrix: they belong in \(\mathbf{b}\), while the design matrix entries are formed from the inputs and their powers.
Check Your Understanding
For each question, keep the coefficient order aligned with the design-matrix columns and check numerical answers by substitution.
- Write the design matrix for a degree-three polynomial using input values \(-1,0,2,3\).
- For a degree-two polynomial and four distinct inputs, explain why the design matrix has independent columns.
- In a reduced QR factorization \(\mathbf{A}=\mathbf{Q}\mathbf{R}\), what triangular system must be solved to obtain the least-squares coefficient vector?
- For the linear fit above, verify the residual sum of squares directly from the three residual entries.
- Explain why \(\mathbf{Q}^{\top}\mathbf{b}\) is generally not the coefficient vector in the original polynomial basis.
