“Of all the principles which can be proposed for that purpose, I think there is none more general, more exact, and more easy of application, than that of which we made use in the preceding researches, and which consists of rendering the sum of squares of the errors a minimum. By this means, there is established among the errors a sort of equilibrium which, preventing the extremes from exerting an undue influence, is very well fitted to reveal that state of the system which most nearly approaches the truth.”
— Adrien-Marie Legendre, Nouvelles méthodes pour la détermination des orbites des comètes, 1805, translated by Henry A. Ruger and Helen M. Walker
A system \(\bm A\bm x = \bm b\) with more equations than unknowns, \(m > n\), is overdetermined. Measurements produce such systems: every measured point adds an equation while the model keeps its few unknowns, and the errors in the measurements make the equations contradict each other. In the column picture of Chapter 3.12, \(\bm A\bm x\) is a combination of the \(n\) columns of \(\bm A\) and runs over their column space as \(\bm x\) varies. A right-hand side \(\bm b\) outside the column space is reached by no \(\bm x\), and the system has no solution.
We then ask for the \(\bm x\) that comes closest. The residual
\[ \bm r = \bm b - \bm A\bm x \tag{3.15.1}\]
holds one miss per equation, and we choose \(\hat{\bm x}\) so that the sum of the squared misses, the squared length of \(\bm r\), is least,
\[ \boxed{\lvert\bm b - \bm A\hat{\bm x}\rvert^2 = \min_{\bm x}\,\lvert\bm b - \bm A\bm x\rvert^2} \tag{3.15.2}\]
The hat marks the least-squares solution. It satisfies none of the equations exactly; it is the best compromise between all of them.
The three equations
\[ x_1 = 1, \qquad x_2 = 1, \qquad x_1 + x_2 = 3 \]
in the two unknowns \(x_1\) and \(x_2\) show both pictures at a small size. In the row picture, Figure 3.15.1, each equation is a line in the \(x_1x_2\) plane. Any two of the lines meet, at \((1, 1)\), \((2, 1)\) or \((1, 2)\), and the third line misses that point. No point lies on all three.
On matrix form the system reads
\[ \bm A\bm x = \begin{bmatrix}1 & 0\\ 0 & 1\\ 1 & 1\end{bmatrix}\begin{bmatrix}x_1\\ x_2\end{bmatrix} = \begin{bmatrix}1\\ 1\\ 3\end{bmatrix} = \bm b \]
and in the column picture, Figure 3.15.2, it asks for the numbers \(x_1\) and \(x_2\) that make \(x_1\bm a_1 + x_2\bm a_2 = \bm b\). The columns \(\bm a_1 = (1, 0, 1)\) and \(\bm a_2 = (0, 1, 1)\) span a plane through the origin of space, the vectors whose third component is the sum of the first two. For \(\bm b = (1, 1, 3)\) the sum is \(2\), not \(3\), so \(\bm b\) lies above the plane.
The column picture turns the minimum of 3.15.2 into geometry. The length \(\lvert\bm b - \bm A\bm x\rvert\) is the distance from the tip of \(\bm b\) to the point \(\bm A\bm x\) in the column space, and the shortest distance from a point to a plane is measured along the perpendicular, as for the point and the plane in Chapter 3.8. The best residual \(\hat{\bm r} = \bm b - \bm A\hat{\bm x}\) is therefore perpendicular to the column space, which means perpendicular to every column,
\[ \bm a_k\cdot\hat{\bm r} = \bm a_k^\mathsf{T}\hat{\bm r} = 0, \qquad k = 1, \dots, n \]
The rows of \(\bm A^\mathsf{T}\) are the columns of \(\bm A\), so these \(n\) dot products are the \(n\) rows of the product \(\bm A^\mathsf{T}\hat{\bm r}\), and the conditions read \(\bm A^\mathsf{T}(\bm b - \bm A\hat{\bm x}) = \bm 0\). Multiplying out gives the normal equations
\[ \bm A^\mathsf{T}\bm A\,\hat{\bm x} = \bm A^\mathsf{T}\bm b \tag{3.15.3}\]
named after the residual, which is normal to the column space. The matrix \(\bm A^\mathsf{T}\bm A\) is \(n\times n\) and symmetric, Chapter 3.10, so the \(m\) contradictory equations have become \(n\) equations in \(n\) unknowns, solved as in Chapter 3.12. When the columns of \(\bm A\) are linearly independent, Chapter 3.6, the solution is unique: \(\bm A^\mathsf{T}\bm A\bm x = \bm 0\) gives \(\bm x^\mathsf{T}\bm A^\mathsf{T}\bm A\bm x = \lvert\bm A\bm x\rvert^2 = 0\), so \(\bm A\bm x = \bm 0\) and \(\bm x = \bm 0\), and \(\bm A^\mathsf{T}\bm A\) is invertible, Chapter 3.11.
The perpendicular residual is also the shortest one. Any other \(\bm x\) leaves the residual \(\bm b - \bm A\bm x = \hat{\bm r} + \bm A(\hat{\bm x} - \bm x)\), whose second part lies in the column space and is perpendicular to the first. Multiplying out the dot product of the sum with itself, the mixed terms vanish and
\[ \lvert\bm b - \bm A\bm x\rvert^2 = \lvert\hat{\bm r}\rvert^2 + \lvert\bm A(\hat{\bm x} - \bm x)\rvert^2 \ge \lvert\hat{\bm r}\rvert^2 \]
with equality only for \(\bm A\bm x = \bm A\hat{\bm x}\). In the video, \(\bm x\) wanders and \(\bm A\bm x\) runs over the plane of the columns of Figure 3.15.2. The length of the residual falls and rises again, and it is least where the residual meets the plane at a right angle.
With a single column \(\bm a\) the column space shrinks to a line. Then \(\bm A^\mathsf{T}\bm A = \bm a\cdot\bm a\) and \(\bm A^\mathsf{T}\bm b = \bm a\cdot\bm b\) are numbers, and 3.15.3 gives
\[ \hat x = \frac{\bm a\cdot\bm b}{\bm a\cdot\bm a} \tag{3.15.4}\]
the factor \(s\) of the projection in 3.6.7. Least squares is the projection of Chapter 3.6 with a plane, or a space of \(n\) columns, in place of the line, and \(\bm p = \bm A\hat{\bm x}\) is the projection of \(\bm b\) on the column space.
The derivative route of Chapter 3.8 arrives at the same equations. Multiplied out, the squared length of the residual is
\[ E(\bm x) = (\bm b - \bm A\bm x)^\mathsf{T}(\bm b - \bm A\bm x) = \bm b^\mathsf{T}\bm b - 2\,\bm x^\mathsf{T}\bm A^\mathsf{T}\bm b + \bm x^\mathsf{T}\bm A^\mathsf{T}\bm A\,\bm x \]
where the two mixed terms are equal because a \(1\times 1\) matrix is its own transpose. \(E\) is quadratic in the unknowns, so its partial derivatives are linear in them, and setting the \(n\) derivatives to zero gives 3.15.3 again. Example 1 checks this with SymPy.
Practise 3.15.4 below: the model \(y = c\,x\) has one column, the \(x\) of the points, and the least-squares slope is a projection.
Find the least-squares solution of the system of Figure 3.15.1,
\[ \begin{bmatrix}1 & 0\\ 0 & 1\\ 1 & 1\end{bmatrix}\begin{bmatrix}x_1\\ x_2\end{bmatrix} = \begin{bmatrix}1\\ 1\\ 3\end{bmatrix} \]
and check that the residual is perpendicular to both columns.
The two products of the normal equations are
\[ \bm A^\mathsf{T}\bm A = \begin{bmatrix}1 & 0 & 1\\ 0 & 1 & 1\end{bmatrix}\begin{bmatrix}1 & 0\\ 0 & 1\\ 1 & 1\end{bmatrix} = \begin{bmatrix}2 & 1\\ 1 & 2\end{bmatrix}, \qquad \bm A^\mathsf{T}\bm b = \begin{bmatrix}1 & 0 & 1\\ 0 & 1 & 1\end{bmatrix}\begin{bmatrix}1\\ 1\\ 3\end{bmatrix} = \begin{bmatrix}4\\ 4\end{bmatrix} \]
so 3.15.3 is the square system
\[ \begin{aligned} 2\hat x_1 + \hat x_2 &= 4\\ \hat x_1 + 2\hat x_2 &= 4 \end{aligned} \]
whose solution is \(\hat x_1 = \hat x_2 = 4/3\). The residual follows from 3.15.1, \(\hat{\bm r} = \bm b - \bm A\hat{\bm x}\), and the check is \(\bm A^\mathsf{T}\hat{\bm r} = \bm 0\). As a second check, the derivative route sets \(\partial E/\partial x_1 = \partial E/\partial x_2 = 0\) with \(E = \lvert\bm b - \bm A\bm x\rvert^2\).
x_1, x_2 = sp.symbols('x_1 x_2', real=True)
AA = sp.Matrix([[1, 0], [0, 1], [1, 1]])
bb = sp.Matrix([1, 1, 3])
AtA = AA.T*AA
Atb = AA.T*bb
xx_hat = AtA.solve(Atb) # the n x n normal equations, n = 2
rr = bb - AA*xx_hat # one miss per equation
check = AA.T*rr # perpendicular to both columns
# the derivative route: the same point
xx = sp.Matrix([x_1, x_2])
E = (bb - AA*xx).dot(bb - AA*xx)
sol = sp.solve([E.diff(x_1), E.diff(x_2)], [x_1, x_2])The normal equations give the solution, the residual and the check
\[ \begin{aligned}\bm A^\mathsf{T}\bm A &= \left[\begin{matrix}2 & 1\\1 & 2\end{matrix}\right],\quad \bm A^\mathsf{T}\bm b = \left[\begin{matrix}4\\4\end{matrix}\right]\\ \hat{\bm x} &= \left[\begin{matrix}\dfrac{4}{3}\\\dfrac{4}{3}\end{matrix}\right],\quad \hat{\bm r} = \left[\begin{matrix}- \dfrac{1}{3}\\- \dfrac{1}{3}\\\dfrac{1}{3}\end{matrix}\right],\quad \bm A^\mathsf{T}\hat{\bm r} = \left[\begin{matrix}0\\0\end{matrix}\right]\end{aligned} \]
while the derivative route and the least sum of squares are
\[ \begin{aligned}E &= 2 x_{1}^{2} + 2 x_{1} x_{2} - 8 x_{1} + 2 x_{2}^{2} - 8 x_{2} + 11\\ \hat{\bm x} &= \left[\begin{matrix}\dfrac{4}{3} & \dfrac{4}{3}\end{matrix}\right]^\mathsf{T},\quad \lvert\hat{\bm r}\rvert^2 = \dfrac{1}{3}\end{aligned} \]
Each equation is missed by one third, and the least-squares point \((4/3, 4/3)\) is the centroid of the triangle in Figure 3.15.1. The residual \((-1, -1, 1)/3\) points along \(\bm a_1\times\bm a_2 = (-1, -1, 1)\), the normal of the column plane.
Practise the normal equations of a system with three equations and two unknowns below.
A probe measured the height \(z\) of a plate at six points, Figure 3.15.3, the data of Example 11 in Chapter 3.8, all in millimetres.
| \(x\) | 0 | 200 | 400 | 0 | 200 | 400 |
|---|---|---|---|---|---|---|
| \(y\) | 0 | 0 | 0 | 300 | 300 | 300 |
| \(z\) | 0.10 | 0.22 | 0.31 | 0.43 | 0.48 | 0.62 |
Fit the plane \(z = c_0 + c_1 x + c_2 y\) with the normal equations and compare with the derivative route of that chapter.
Each measured point gives one equation \(c_0 + c_1 x_i + c_2 y_i = z_i\), six equations in three unknowns,
\[ \begin{bmatrix} 1 & 0 & 0\\ 1 & 200 & 0\\ 1 & 400 & 0\\ 1 & 0 & 300\\ 1 & 200 & 300\\ 1 & 400 & 300 \end{bmatrix} \begin{bmatrix}c_0\\ c_1\\ c_2\end{bmatrix} = \begin{bmatrix}0.10\\ 0.22\\ 0.31\\ 0.43\\ 0.48\\ 0.62\end{bmatrix} \]
The columns of \(\bm A\) are a column of ones, the \(x\) and the \(y\) of the points. The first row of \(\bm A^\mathsf{T}\hat{\bm r} = \bm 0\) is the dot product of the column of ones with the residual, \(\sum e_i = 0\), the property that Chapter 3.8 found from the derivative with respect to \(c_0\). The derivative route, 3.15.3 and the difference between the two are computed below.
xs = [0, 200, 400, 0, 200, 400]
ys = [0, 0, 0, 300, 300, 300]
zs = [sp.Rational(h) for h in ('0.10', '0.22', '0.31', '0.43', '0.48', '0.62')]
AA = sp.Matrix([[1, x_i, y_i] for x_i, y_i in zip(xs, ys)]) # one row per point
bb = sp.Matrix(zs)
cc = (AA.T*AA).solve(AA.T*bb)
ee = bb - AA*cc
# the derivative route of the lines-and-planes chapter
c_0, c_1, c_2 = sp.symbols('c_0 c_1 c_2', real=True)
E = sum((z_i - (c_0 + c_1*x_i + c_2*y_i))**2 for x_i, y_i, z_i in zip(xs, ys, zs))
sol = sp.solve([E.diff(c) for c in (c_0, c_1, c_2)], [c_0, c_1, c_2])
difference = cc - sp.Matrix([sol[c_0], sol[c_1], sol[c_2]])The normal equations and the coefficients, with the difference from the derivative route, are
\[ \begin{aligned}\bm A^\mathsf{T}\bm A &= \left[\begin{matrix}6 & 1200 & 900\\1200 & 400000 & 180000\\900 & 180000 & 270000\end{matrix}\right],\quad \bm A^\mathsf{T}\bm b = \left[\begin{matrix}\dfrac{54}{25}\\512\\459\end{matrix}\right]\\ \hat{\bm c} &= \left[\begin{matrix}\dfrac{11}{100}\\\dfrac{1}{2000}\\\dfrac{1}{1000}\end{matrix}\right]\approx\left[\begin{matrix}0.11\\0.0005\\0.001\end{matrix}\right],\quad \hat{\bm c} - \bm c_{\text{derivatives}} = \left[\begin{matrix}0\\0\\0\end{matrix}\right]\end{aligned} \]
and the residuals, in millimetres, add up to zero:
\[ \hat{\bm e} = \left[\begin{matrix}-0.01 & 0.01 & 0 & 0.02 & -0.03 & 0.01\end{matrix}\right]^\mathsf{T},\qquad \textstyle\sum e_i = 0 \]
The two routes give the same exact fractions, \(z = 0.11 + 0.0005\,x + 0.001\,y\) to three digits: the plane of Example 11 and its flatness of 0.05 mm follow from either. The normal equations ask for no derivatives, and a seventh measured point would only add a row to \(\bm A\) and \(\bm b\).
The residuals here are vertical: only the height \(z\) is taken to carry errors, while \(x\) and \(y\) are exact. When every coordinate of the points is measured with error, the distance to minimise is the perpendicular distance to the plane, and the normal of the best plane is an eigenvector of a \(3\times 3\) matrix built from the points, Chapter 3.14.
Nothing in Example 2 was special to planes. Least squares fits any model that is linear in its unknown coefficients,
\[ y = c_1 f_1(x) + c_2 f_2(x) + \dots + c_n f_n(x) \]
with known functions \(f_k\), which may be as curved as we like. Each data point \((x_i, y_i)\) gives one row, and the matrix of the system, the design matrix, holds \(f_k(x_i)\) in row \(i\) and column \(k\),
\[ \bm A = \begin{bmatrix} f_1(x_1) & f_2(x_1) & \cdots & f_n(x_1)\\ \vdots & \vdots & & \vdots\\ f_1(x_m) & f_2(x_m) & \cdots & f_n(x_m) \end{bmatrix}, \qquad \bm b = \begin{bmatrix}y_1\\ \vdots\\ y_m\end{bmatrix} \tag{3.15.5}\]
For a polynomial the functions are \(1, x, \dots, x^{n-1}\), and for a signal that repeats with the angular frequency \(\omega\) they are \(1\), \(\cos\omega t\) and \(\sin\omega t\). Only the design matrix changes from model to model; 3.15.3 stays.
Fit the straight line \(y = c_0 + c_1 x\) to the six points of Figure 3.15.4, \(x_i = 0, 1, \dots, 5\) and \(y_i = 1, 2, 4, 4, 6, 7\).
With \(f_1 = 1\) and \(f_2 = x\), the design matrix has a column of ones and the column of the \(x_i\). Its products are sums over the points,
\[ \bm A^\mathsf{T}\bm A = \begin{bmatrix}m & \sum x_i\\ \sum x_i & \sum x_i^2\end{bmatrix} = \begin{bmatrix}6 & 15\\ 15 & 55\end{bmatrix}, \qquad \bm A^\mathsf{T}\bm b = \begin{bmatrix}\sum y_i\\ \sum x_i y_i\end{bmatrix} = \begin{bmatrix}24\\ 81\end{bmatrix} \]
and the normal equations
\[ \begin{aligned} 6c_0 + 15c_1 &= 24\\ 15c_0 + 55c_1 &= 81 \end{aligned} \]
give \(c_1 = 6/5\) and \(c_0 = 1\).
xs = [0, 1, 2, 3, 4, 5]
ys = [1, 2, 4, 4, 6, 7]
AA = sp.Matrix([[1, x_i] for x_i in xs]) # columns f_1 = 1 and f_2 = x
bb = sp.Matrix(ys)
cc = (AA.T*AA).solve(AA.T*bb)
rr = bb - AA*ccSolving the normal equations gives the line and its residuals,
\[ \begin{aligned}\hat{\bm c} &= \left[\begin{matrix}1 & \dfrac{6}{5}\end{matrix}\right]^\mathsf{T},\quad y = \dfrac{6 x}{5} + 1\\ \hat{\bm r} &= \left[\begin{matrix}0 & - \dfrac{1}{5} & \dfrac{3}{5} & - \dfrac{3}{5} & \dfrac{1}{5} & 0\end{matrix}\right]^\mathsf{T},\quad \textstyle\sum r_i^2 = \dfrac{4}{5}\end{aligned} \]
The residuals add up to zero, as the column of ones demands, and the line passes exactly through the first and the last point. Figure 3.15.5 draws the fit over the data.
x = np.array(xs, float)
y = np.array(ys, float)
c = np.array(cc, float).ravel()
x_line = np.array([0, 5.5])
fig, ax = plt.subplots(figsize=(5.0, 3.2))
ax.plot(x_line, c[0] + c[1]*x_line, color='#4472c4', lw=2, label='$y = 1 + 1.2x$')
# each residual as the vertical gap between a point and the line
ax.vlines(x, c[0] + c[1]*x, y, color='#2e7d32', lw=2, label='residuals')
ax.plot(x, y, 'o', color='#c00000', ms=7, zorder=5, label='data')
ax.set_xlim(0, 5.5)
ax.set_ylim(0, 8)
ax.set_xticks([0, 1, 2, 3, 4, 5, 5.5], ['0', '1', '2', '3', '4', '5', '5.5'])
ax.set_yticks([0, 2, 4, 6, 8])
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.grid(True, alpha=0.35)
ax.legend(loc='upper left', frameon=False)
fig.tight_layout()
plt.show()
In the video the line moves over the points of Example 3. On every residual stands a square, so the shaded area is the sum that 3.15.2 makes least. The line passes a point where the sum is smaller than at the start and still not least, and it settles at \(c_0 = 1\), \(c_1 = 1.2\), where no move of the line shrinks the total area further.
Practise the straight-line fit below. Four points give a \(4\times 2\) design matrix, and the normal equations are built from the four sums \(m\), \(\sum x_i\), \(\sum x_i^2\) and \(\sum x_i y_i\) together with \(\sum y_i\).
Fit the parabola \(y = c_0 + c_1 x + c_2 x^2\) to the six points of Figure 3.15.6, \(x_i = 0, 1, \dots, 5\) and \(y_i = 6, 3, 2, 2, 4, 7\).
The design matrix gains the column of the \(x_i^2\), and the normal equations hold sums of powers up to \(x^4\),
\[ \bm A = \begin{bmatrix}1 & x_1 & x_1^2\\ \vdots & \vdots & \vdots\\ 1 & x_6 & x_6^2\end{bmatrix}, \qquad \bm A^\mathsf{T}\bm A = \begin{bmatrix}6 & \sum x_i & \sum x_i^2\\ \sum x_i & \sum x_i^2 & \sum x_i^3\\ \sum x_i^2 & \sum x_i^3 & \sum x_i^4\end{bmatrix}, \qquad \bm A^\mathsf{T}\bm b = \begin{bmatrix}\sum y_i\\ \sum x_i y_i\\ \sum x_i^2 y_i\end{bmatrix} \]
a \(3\times 3\) system that we leave to SymPy.
ys = [6, 3, 2, 2, 4, 7]
AA = sp.Matrix([[1, x_i, x_i**2] for x_i in xs]) # columns 1, x, x^2
bb = sp.Matrix(ys)
cc = (AA.T*AA).solve(AA.T*bb)
rr = bb - AA*ccThe normal equations and their solution are
\[ \begin{aligned}\bm A^\mathsf{T}\bm A &= \left[\begin{matrix}6 & 15 & 55\\15 & 55 & 225\\55 & 225 & 979\end{matrix}\right],\quad \bm A^\mathsf{T}\bm b = \left[\begin{matrix}24\\64\\268\end{matrix}\right]\\ \hat{\bm c} &= \left[\begin{matrix}\dfrac{83}{14}\\- \dfrac{493}{140}\\\dfrac{3}{4}\end{matrix}\right]\approx\left[\begin{matrix}5.93\\-3.52\\0.75\end{matrix}\right],\quad \textstyle\sum r_i^2 = \dfrac{3}{35}\end{aligned} \]
The parabola \(y \approx 5.93 - 3.52x + 0.75x^2\) misses no point by more than 0.16. Every column we add enlarges the column space, so the least sum of squares can only fall. In the video the points of Example 4 are fitted with one to six columns, the polynomials of degree 0 to 5. The sum drops sharply at the third column, where the model first has the shape of the data, and barely after it. With six columns \(\bm A\) is square and invertible, the curve passes through every point and the sum is zero, but the curve now follows the scatter of the points as well as their trend and bends between them. A good model has few columns compared with the number of points.
Practise the design matrix below: write \(\bm A\) for a model and a set of points, row by row, as a list such as [[1, 0, 0], [1, 1, 1], [1, 2, 4]].
A thermometer outdoors reads the temperatures of Figure 3.15.7 every four hours,
| \(t\) (h) | 0 | 4 | 8 | 12 | 16 | 20 |
|---|---|---|---|---|---|---|
| \(T\) (°C) | 8 | 7 | 11 | 16 | 17 | 13 |
Fit the daily cycle \(T(t) = c_0 + c_1\cos\omega t + c_2\sin\omega t\) with \(\omega = 2\pi/24\ \text{h}^{-1}\), and find the mean temperature, the amplitude of the cycle and the time of day of the maximum.
The design matrix has the rows \((1, \cos\omega t_i, \sin\omega t_i)\). At the six times, \(\omega t_i\) runs through \(0°, 60°, \dots, 300°\), so
\[ \bm A = \begin{bmatrix} 1 & 1 & 0\\ 1 & \tfrac12 & \tfrac{\sqrt3}{2}\\ 1 & -\tfrac12 & \tfrac{\sqrt3}{2}\\ 1 & -1 & 0\\ 1 & -\tfrac12 & -\tfrac{\sqrt3}{2}\\ 1 & \tfrac12 & -\tfrac{\sqrt3}{2} \end{bmatrix}, \qquad \bm A^\mathsf{T}\bm A = \begin{bmatrix}6 & 0 & 0\\ 0 & 3 & 0\\ 0 & 0 & 3\end{bmatrix} \]
The columns are orthogonal, as they are for any samples evenly spread over one period, so \(\bm A^\mathsf{T}\bm A\) is diagonal and each normal equation holds one unknown,
\[ \begin{aligned} 6c_0 &= \textstyle\sum T_i = 72\\ 3c_1 &= \textstyle\sum T_i\cos\omega t_i = 8 + \tfrac72 - \tfrac{11}{2} - 16 - \tfrac{17}{2} + \tfrac{13}{2} = -12\\ 3c_2 &= \textstyle\sum T_i\sin\omega t_i = \tfrac{\sqrt3}{2}(7 + 11 - 17 - 13) = -6\sqrt3 \end{aligned} \]
so \(c_0 = 12\), \(c_1 = -4\) and \(c_2 = -2\sqrt3\). The sum of a cosine and a sine of the same frequency is one cosine, shifted,
\[ c_1\cos\omega t + c_2\sin\omega t = R\cos\big(\omega(t - t_{\max})\big), \qquad R = \sqrt{c_1^2 + c_2^2}, \qquad \omega t_{\max} = \operatorname{atan2}(c_2, c_1) \]
with the amplitude \(R\) and the time \(t_{\max}\) of the maximum; atan2 picks the angle in the right quadrant.
t = sp.symbols('t', real=True)
w = 2*sp.pi/24
ts = [0, 4, 8, 12, 16, 20]
Ts = [8, 7, 11, 16, 17, 13]
AA = sp.Matrix([[1, sp.cos(w*t_i), sp.sin(w*t_i)] for t_i in ts])
bb = sp.Matrix(Ts)
cc = sp.simplify((AA.T*AA).solve(AA.T*bb))
R = sp.sqrt(cc[1]**2 + cc[2]**2)
phase = sp.atan2(cc[2], cc[1]) % (2*sp.pi) # in [0, 2 pi): a time of day
t_max = phase/wThe coefficients, the amplitude and the time of the maximum are
\[ \begin{aligned}\bm A^\mathsf{T}\bm A &= \left[\begin{matrix}6 & 0 & 0\\0 & 3 & 0\\0 & 0 & 3\end{matrix}\right],\quad \hat{\bm c} = \left[\begin{matrix}12\\-4\\- 2 \sqrt{3}\end{matrix}\right]\\ R &= 2 \sqrt{7}\approx5.29\ ^\circ\text{C},\quad t_{\max} = 14.73\ \text{h}\end{aligned} \]
The day averages 12 °C and swings 5.3 °C to either side, with the warmest moment at 14.7 h, about a quarter to three in the afternoon. The model \(T_0 + R\cos(\omega(t - t_{\max}))\) is not linear in \(t_{\max}\); written with a cosine and a sine it is linear in \(c_1\) and \(c_2\), and least squares applies. A model such as \(y = a\,e^{kx}\) is made linear the same way, by fitting \(\ln y = \ln a + kx\). Figure 3.15.8 draws the cycle through the readings.
T_fit = sp.lambdify(t, cc[0] + cc[1]*sp.cos(w*t) + cc[2]*sp.sin(w*t), 'numpy')
t_plot = np.linspace(0, 24, 200)
fig, ax = plt.subplots(figsize=(5.0, 3.2))
ax.plot(t_plot, T_fit(t_plot), color='#4472c4', lw=2, label='least-squares cycle')
ax.plot(ts, Ts, 'o', color='#c00000', ms=7, zorder=5, label='readings')
ax.axvline(float(t_max), color='0.35', lw=0.8, ls='--')
ax.set_xlim(0, 24)
ax.set_ylim(0, 20)
ax.set_xticks([0, 4, 8, 12, 16, 20, 24])
ax.set_yticks([0, 5, 10, 15, 20])
ax.set_xlabel('time of day $t$ (h)')
ax.set_ylabel('temperature $T$ (°C)')
ax.grid(True, alpha=0.35)
ax.legend(loc='lower right', frameon=False)
fig.tight_layout()
plt.show()
Example 6 of Chapter 3.8 found the point \(\bm p\) nearest to three sight lines, Figure 3.15.9, by setting derivatives to zero. The cameras are at \(C_1 = (6, 0, 0)\), \(C_2 = (0, 6, 0)\) and \(C_3 = (0, 0, 5)\) and look along \(\bm v_1 = (-4, 2, 1)\), \(\bm v_2 = (2, -4, 1)\) and \(\bm v_3 = (2, 2, -3)\). Write the problem on the form \(\bm A\bm p = \bm b\), solve it with the normal equations, and add a fourth camera at \(C_4 = (6, 6, 0)\) looking along \(\bm v_4 = (-3, -3, 1)\).
The distance vector from line \(i\) to \(\bm p\) is what is left of \(\bm a = \bm p - \bm r_{C_i}\) after its projection on the unit vector \(\bm e_i = \bm v_i/\lvert\bm v_i\rvert\) is removed, \(\bm d_i = \bm a - (\bm a\cdot\bm e_i)\,\bm e_i\). The projection is a matrix product: \(\bm e_i\bm e_i^\mathsf{T}\) is a \(3\times 1\) column times a \(1\times 3\) row, a \(3\times 3\) matrix, and \((\bm e_i\bm e_i^\mathsf{T})\bm a = \bm e_i(\bm e_i^\mathsf{T}\bm a) = (\bm a\cdot\bm e_i)\,\bm e_i\) is the cheap projection of 3.6.8. With
\[ \bm P_i = \bm I - \bm e_i\bm e_i^\mathsf{T} \]
the distance vector is \(\bm d_i = \bm P_i(\bm p - \bm r_{C_i})\). A point on every line would have \(\bm d_i = \bm 0\) for all \(i\), three equations \(\bm P_i\bm p = \bm P_i\bm r_{C_i}\) per line, and stacked they form an overdetermined system with nine rows,
\[ \begin{bmatrix}\bm P_1\\ \bm P_2\\ \bm P_3\end{bmatrix}\bm p = \begin{bmatrix}\bm P_1\bm r_{C_1}\\ \bm P_2\bm r_{C_2}\\ \bm P_3\bm r_{C_3}\end{bmatrix} \]
Its residual is the stack of the three distance vectors, so \(\lvert\bm r\rvert^2 = \sum\lvert\bm d_i\rvert^2\) is the sum that Chapter 3.8 minimised. The fourth camera adds three rows. \(\bm A^\mathsf{T}\bm A\) and \(\bm A^\mathsf{T}\bm b\) come out as sums over the lines, since \(\bm P_i\) is symmetric and \(\bm P_i\bm P_i = \bm P_i\) (removing the part along \(\bm e_i\) twice removes nothing more), so the normal equations read
\[ \Big(\sum_i \bm P_i\Big)\,\bm p = \sum_i \bm P_i\,\bm r_{C_i} \]
cameras = [sp.Matrix([6, 0, 0]), sp.Matrix([0, 6, 0]), sp.Matrix([0, 0, 5]), sp.Matrix([6, 6, 0])]
sight = [sp.Matrix([-4, 2, 1]), sp.Matrix([2, -4, 1]), sp.Matrix([2, 2, -3]), sp.Matrix([-3, -3, 1])]
def stacked(n):
"""A and b of the first n lines: the block P_i removes the part along line i."""
PP = [sp.eye(3) - vv*vv.T/vv.dot(vv) for vv in sight[:n]] # e e^T = v v^T / v.v
AA = sp.Matrix.vstack(*PP)
bb = sp.Matrix.vstack(*[P*rr_C for P, rr_C in zip(PP, cameras[:n])])
return AA, bb
AA, bb = stacked(3)
pp_3 = (AA.T*AA).solve(AA.T*bb)
same_as_sums = AA.T*AA - sum((sp.eye(3) - vv*vv.T/vv.dot(vv) for vv in sight[:3]), sp.zeros(3))
AA, bb = stacked(4)
pp_4 = (AA.T*AA).solve(AA.T*bb)
dd_4 = [(bb - AA*pp_4)[3*i:3*i + 3, 0].norm() for i in range(4)]With the three lines of Chapter 3.8, and with the fourth line added, the normal equations give
\[ \begin{aligned}\bm p_{3\text{ lines}} &= \left[\begin{matrix}\dfrac{244}{115}\\\dfrac{244}{115}\\\dfrac{53}{46}\end{matrix}\right]\approx\left[\begin{matrix}2.122\\2.122\\1.152\end{matrix}\right],\quad \bm A^\mathsf{T}\bm A - \textstyle\sum_i\bm P_i = \left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]\\ \bm p_{4\text{ lines}} &\approx \left[\begin{matrix}2.121\\2.121\\1.193\end{matrix}\right],\quad \lvert\bm d_i\rvert \approx \left[\begin{matrix}0.26 & 0.26 & 0.43 & 0.098\end{matrix}\right]\end{aligned} \]
The three-line point is the exact point \((244/115, 244/115, 53/46)\) of Chapter 3.8. The fourth sight line passes close to it, within 0.10, and moves the point by less than 0.05. In matrix form a further measurement is three more rows and one more term in each sum, and the system stays \(3\times 3\) however many lines are measured.
The normal equations are the way to solve by hand and to reason about least squares, and they serve for small, well-scaled problems. On a computer, forming \(\bm A^\mathsf{T}\bm A\) squares the condition number of the problem, Chapter 3.12: \(\kappa(\bm A^\mathsf{T}\bm A) = \kappa(\bm A)^2\).
x = np.arange(6.0)
y = np.array([6, 3, 2, 2, 4, 7.0])
conds = []
for n in (2, 3, 5):
A = np.vander(x, n, increasing=True) # columns 1, x, ..., x^(n-1)
conds.append([n, np.linalg.cond(A), np.linalg.cond(A.T @ A)])
conds = np.array(conds)
A = np.vander(x, 3, increasing=True)
c, res, rank, sv = np.linalg.lstsq(A, y, rcond=None)For the points of Example 4 fitted with \(n\) columns, the condition numbers are
\[ \begin{array}{c|cc} n & \kappa(\bm A) & \kappa(\bm A^\mathsf{T}\bm A) \\ \hline 2 & 5.78 & 33.4 \\ 3 & 37.8 & 1.43 \cdot 10^{3} \\ 5 & 3.59 \cdot 10^{3} & 1.29 \cdot 10^{7}\end{array} \]
With five columns \(\bm A^\mathsf{T}\bm A\) loses about seven of the sixteen digits of a double-precision number before the solution starts. np.linalg.lstsq(A, b, rcond=None) avoids \(\bm A^\mathsf{T}\bm A\) altogether: it factorises \(\bm A\) itself into orthogonal and triangular factors and solves 3.15.2 with the condition number of \(\bm A\). It returns the coefficients, the least sum of squares, the rank of \(\bm A\) and its singular values, and for the parabola of Example 4 it gives
\[ \hat{\bm c} \approx \begin{bmatrix}5.929 \\ -3.521 \\ 0.75\end{bmatrix},\qquad \lvert\hat{\bm r}\rvert^2 \approx 0.08571,\qquad \text{rank} = 3 \]
the coefficients and the sum \(3/35 \approx 0.0857\) of Example 4. In SymPy, exact arithmetic has no round-off and (AA.T*AA).solve(AA.T*bb) is safe, as in the examples above.