7.2 Discrete least-squares approximation
Suppose the data \((x_i,y_i)\), \(i = 1,\dots ,m\), come from measurements and contain errors. An interpolating polynomial would reproduce the errors. Instead we choose a function from a simple family, for instance a straight line \(y = a_1x+a_0\), that minimizes the sum of the squared deviations \[E(a_0,a_1) = \sum _{i=1}^m\bigl (y_i-(a_1x_i+a_0)\bigr )^2 .\] Squares are used because they make \(E\) a smooth function of the parameters, they penalize large deviations heavily, and (as the statistics courses show) they give the best estimates when the errors are independent and normally distributed.
Theorem 7.2.1. The line \(y = a_1x+a_0\) minimizing \(E\) is given by the solution of the normal equations \[a_0m+a_1\sum x_i = \sum y_i,\qquad a_0\sum x_i+a_1\sum x_i^2 = \sum x_iy_i ,\] that is, \[a_1 = \frac {m\sum x_iy_i-\sum x_i\sum y_i}{m\sum x_i^2-\left (\sum x_i\right )^2},\qquad a_0 = \frac {\sum y_i-a_1\sum x_i}{m} .\] More generally, the polynomial \(P_n(x) = \sum _{k=0}^na_kx^k\) (\(n<m-1\)) minimizing \(E = \sum _i(y_i-P_n(x_i))^2\) satisfies the \(n+1\) normal equations \[\sum _{k=0}^na_k\sum _{i=1}^mx_i^{j+k} = \sum _{i=1}^my_ix_i^j,\qquad j = 0,1,\dots ,n,\] which have a unique solution when the \(x_i\) are distinct.
Proof. \(E\) is a quadratic function of \(a_0,\dots ,a_n\), bounded below by \(0\), so at a minimum \(\frac {\partial E}{\partial a_j} = 0\) for each \(j\): \[0 = \frac {\partial E}{\partial a_j} = -2\sum _i\bigl (y_i-P_n(x_i)\bigr )x_i^j ,\] which rearranges to the normal equations. In matrix form, with \(X\) the \(m\times (n+1)\) matrix of entries \(x_i^k\), the equations are \(X^TX\mathbf a = X^T\mathbf y\). If the \(x_i\) are distinct and \(m>n\), then \(X\mathbf v = \mathbf 0\) means that the polynomial with coefficients \(\mathbf v\) vanishes at \(m>n\) points, so \(\mathbf v = \mathbf 0\); hence \(X^TX\) is positive definite (\(\mathbf v^TX^TX\mathbf v = \|X\mathbf v\|_2^2>0\)), the solution is unique, and since \(E\) is then a strictly convex quadratic, it is the minimum. □
Example 7.2.2. Find the least-squares line for the data
| \(x_i\) | \(1\) | \(2\) | \(3\) | \(4\) | \(5\) | \(6\) | \(7\) | \(8\) | \(9\) | \(10\) |
| \(y_i\) | \(1.3\) | \(3.5\) | \(4.2\) | \(5.0\) | \(7.0\) | \(8.8\) | \(10.1\) | \(12.5\) | \(13.0\) | \(15.6\) |
Solution. The sums are \(\sum x_i = 55\), \(\sum x_i^2 = 385\), \(\sum y_i = 81.0\), \(\sum x_iy_i = 572.4\). So \begin {gather*} a_1 = \frac {10(572.4)-55(81.0)}{10(385)-55^2} = \frac {1269}{825} = 1.538,\\ a_0 = \frac {81.0-1.538(55)}{10} = -0.360, \end {gather*} and the line is \(y = 1.538x-0.360\), with \(E = 2.34\). The interpolating polynomial of degree \(9\) through these ten points would pass through every one of them, but would oscillate between them; the line captures the trend.
Example 7.2.3. Fit a quadratic \(P_2(x) = a_0+a_1x+a_2x^2\) by least squares to
| \(x_i\) | \(0\) | \(0.25\) | \(0.50\) | \(0.75\) | \(1.00\) |
| \(y_i\) | \(1.0000\) | \(1.2840\) | \(1.6487\) | \(2.1170\) | \(2.7183\) |
Solution. The power sums are \(\sum x_i^k = 5, 2.5, 1.875, 1.5625, 1.3828\) for \(k = 0,\dots ,4\), and \(\sum y_ix_i^k = 8.7680, 5.4514, 4.4015\) for \(k = 0,1,2\). The normal equations \begin {align*} 5a_0+2.5a_1+1.875a_2 &= 8.7680,\\ 2.5a_0+1.875a_1+1.5625a_2 &= 5.4514,\\ 1.875a_0+1.5625a_1+1.3828a_2 &= 4.4015 \end {align*}
have the solution \(a_0 = 1.00514\), \(a_1 = 0.86418\), \(a_2 = 0.84366\), so \(P_2(x) = 1.00514+0.86418x+0.84366x^2\), with \(E = 2.74\times 10^{-4}\). (The data are values of \(e^x\), and the quadratic follows them to within about \(0.01\).)
Fitting other forms
Data that grow exponentially are often fitted by \(y = be^{ax}\). Minimizing \(\sum (y_i-be^{ax_i})^2\) leads to non-linear equations, so in practice one takes logarithms, \(\ln y = \ln b+ax\), and fits a straight line to \((x_i,\ln y_i)\). This minimizes the squared errors of the logarithms rather than of the data, which is a slightly different criterion, but it is simple and usually adequate.
Example 7.2.4. Fit \(y = be^{ax}\) to the data
| \(x_i\) | \(1.00\) | \(1.25\) | \(1.50\) | \(1.75\) | \(2.00\) |
| \(y_i\) | \(5.10\) | \(5.79\) | \(6.53\) | \(7.45\) | \(8.46\) |
Solution. With \(Y_i = \ln y_i\): \(\sum x_i = 7.5\), \(\sum x_i^2 = 11.875\), \(\sum Y_i = 9.40534\), \(\sum x_iY_i = 14.42409\). The normal equations give \[a = \frac {5(14.42409)-7.5(9.40534)}{5(11.875)-7.5^2} = 0.5057,\qquad \ln b = \frac {9.40534-0.5057(7.5)}{5} = 1.1225,\] so \(b = 3.072\) and \(y = 3.072e^{0.5057x}\). The sum of squared deviations of the data from this curve is \(1.2\times 10^{-3}\).
Note. The normal-equation matrix for polynomial fitting on \([0,1]\) with many equally spaced points is close to a multiple of the Hilbert matrix of Example 5.5.4, and becomes very ill-conditioned as the degree grows. For degrees beyond about \(5\), least-squares problems are solved instead with orthogonal polynomials or with a QR factorization of \(X\), which avoid forming \(X^TX\).
Questions on this section
Stuck on something here? Ask below and it stays attached to this topic.