Linear Least Squares and the Normal Equations
~9 min read · Tier toggle changes depth (T1 plain · T2 standard · T3 rigour)
You should be comfortable with matrix algebra, specifically the transpose, matrix multiplication, and the concept of a matrix inverse. You must also understand the definition of the chi-squared (\(\chi^2\)) statistic as a measure of goodness-of-fit.
Why this matters
This is the mathematical engine that turns raw data into physical constants.
In experimental physics, we rarely measure a single number; we measure a trend. Whether you are determining the slope of a velocity-time graph to find acceleration or fitting a decay curve to find a half-life, you are solving for unknown parameters that best describe a set of noisy observations. The method of least squares provides the unique, mathematically rigorous solution to this problem. By deriving the "normal equations," we move from a geometric intuition of "minimising the distance to the line" to a solvable linear system that can be computed by any modern computer, forming the backbone of data analysis in every subsequent unit.
The Geometric Intuition
Imagine you have plotted a set of experimental points on a graph. You suspect they follow a straight line, but because of measurement noise, they scatter around the true trend. Your goal is to draw the "best" line through them. What does "best" mean? In the context of least squares, it means the line that minimises the sum of the squared vertical distances (residuals) between the data points and the line. We square the distances to ensure that positive and negative deviations do not cancel each other out and to penalise large errors more heavily than small ones. This is not just a convention; as established in the previous lecture, if the noise is Gaussian, minimising the sum of squared residuals is mathematically equivalent to finding the maximum likelihood estimate. The power of this approach lies in its linearity: although the data is noisy, the relationship between the parameters and the data is linear, allowing us to solve for the parameters using standard linear algebra rather than complex iterative numerical methods.
Formulating the Problem
To apply this to general physics problems, we must express the model in matrix notation. Suppose we have \(N\) measurements of a variable \(y\), and our model depends on \(M\) unknown parameters \(\mathbf{p}\). We assume the model is linear in these parameters. The relationship can be written as:
Here, \(\mathbf{y}\) is an \(N \times 1\) column vector of observed values. \(A\) is the \(N \times M\) design matrix, where each row corresponds to a specific measurement condition (e.g., the time \(t\) in a linear fit \(y = mt + c\)). \(\mathbf{p}\) is the \(M \times 1\) vector of parameters we wish to find (e.g., slope and intercept). \(\boldsymbol{\epsilon}\) represents the noise. We define the residual vector \(\mathbf{r} = \mathbf{y} - A\mathbf{p}\). The quantity we wish to minimise is the weighted chi-squared:
The matrix \(W\) is the weight matrix, typically the inverse of the covariance matrix of the measurements. If all measurements have equal uncertainty \(\sigma^2\), then \(W = \frac{1}{\sigma^2} I\), where \(I\) is the identity matrix. Minimising \(\chi^2\) with respect to \(\mathbf{p}\) yields the optimal parameters.
Deriving the Normal Equations
To find the minimum, we take the derivative of \(\chi^2\) with respect to the parameter vector \(\mathbf{p}\) and set it to zero. Using matrix calculus rules, specifically \(\frac{\partial}{\partial \mathbf{p}} (\mathbf{a}^T \mathbf{p}) = \mathbf{a}\) and \(\frac{\partial}{\partial \mathbf{p}} (\mathbf{p}^T B \mathbf{p}) = 2B\mathbf{p}\) for symmetric \(B\), we expand the expression. The derivative of the cross terms yields \(A^T W \mathbf{y}\), and the derivative of the quadratic term in \(\mathbf{p}\) yields \(A^T W A \mathbf{p}\). Setting the gradient to zero gives the normal equations:
This is a system of \(M\) linear equations with \(M\) unknowns. Provided that the matrix \(A^T W A\) is non-singular (which is true if the columns of \(A\) are linearly independent and \(W\) is positive definite), we can solve for \(\mathbf{p}\) by multiplying both sides by the inverse:
This result, \(\hat{\mathbf{p}}\), is the least-squares estimate of the parameters. It is unique; for independent errors of known variance it is the best linear unbiased estimator (Gauss–Markov theorem), and for Gaussian noise it is also the maximum-likelihood estimate.
Uncertainty in the Parameters
Finding the best-fit values is only half the story; a physicist must also report the uncertainty in those values. The covariance matrix of the estimated parameters, \(\mathrm{Cov}(\hat{\mathbf{p}})\), is directly related to the inverse of the matrix that appears in the normal equations. Specifically:
If the weights are based on known variances \(\sigma_i^2\), this formula gives the actual variances of the parameters. The diagonal elements of this matrix are the variances of the individual parameters, and the square roots of these diagonal elements are the standard uncertainties. The off-diagonal elements represent the covariances, indicating how the uncertainty in one parameter is correlated with the uncertainty in another. For example, in a linear fit \(y = mx + c\), the slope and intercept are often negatively correlated; if the slope is slightly higher than the true value, the intercept must be slightly lower to keep the line passing through the data.
The derivation
Linear Least Squares and the Normal Equations
For a model that is linear in its parameters, y = Ap + noise, the weighted least-squares estimate p̂ that minimises χ2 = (y − Ap)TW(y − Ap) satisfies the normal equations (ATWA)p̂ = ATWy, where A is the N×M design matrix, y is the vector of N observations, and W is a symmetric positive-definite weight matrix.
A common pitfall is assuming that the diagonal elements of the covariance matrix are independent of the choice of origin. While the variance of the slope in a linear fit is invariant to shifting the x-axis, the variance of the intercept is not. If you fit \(y = mx + c\) versus \(y = m(x - x_0) + c'\), the parameter \(c'\) represents the value at \(x_0\). The uncertainty in \(c'\) will be smaller if \(x_0\) is chosen near the centroid of the data, because the lever arm for the slope uncertainty is reduced. Always report parameters relative to a physically meaningful origin to minimise artificial inflation of uncertainties due to correlation. C
Worked example
Problem. A student measures the position \(x\) (in metres) of a particle at three time points \(t\) (in seconds). The data is: \(t_1 = 1.0\), \(x_1 = 2.1\) \(t_2 = 2.0\), \(x_2 = 4.0\) \(t_3 = 3.0\), \(x_3 = 6.2\) Assume the motion is linear, \(x = vt + x_0\), and that all measurements have an equal standard uncertainty \(\sigma = 0.5\ \mathrm{m}\). Calculate the best-fit velocity \(v\) and initial position \(x_0\), and the uncertainty in the velocity.
- Construct the design matrix \(A\). Since the model is \(x = v t + 1 \cdot x_0\), the columns of \(A\) correspond to \(t\) and \(1\). \(A = \begin{pmatrix} 1.0 & 1 \\ 2.0 & 1 \\ 3.0 & 1 \end{pmatrix}\)
- Construct the observation vector \(\mathbf{y}\). \(\mathbf{y} = \begin{pmatrix} 2.1 \\ 4.0 \\ 6.2 \end{pmatrix}\)
- Determine the weight matrix \(W\). Since \(\sigma = 0.5\), \(W = \frac{1}{0.5^2} I = 4 I\). \(W = \begin{pmatrix} 4 & 0 & 0 \\ 0 & 4 & 0 \\ 0 & 0 & 4 \end{pmatrix}\)
- Compute \(A^T W A\). \(A^T W = \begin{pmatrix} 1.0 & 2.0 & 3.0 \\ 1 & 1 & 1 \end{pmatrix} \begin{pmatrix} 4 & 0 & 0 \\ 0 & 4 & 0 \\ 0 & 0 & 4 \end{pmatrix} = \begin{pmatrix} 4 & 8 & 12 \\ 4 & 4 & 4 \end{pmatrix}\) \(A^T W A = \begin{pmatrix} 4 & 8 & 12 \\ 4 & 4 & 4 \end{pmatrix} \begin{pmatrix} 1.0 & 1 \\ 2.0 & 1 \\ 3.0 & 1 \end{pmatrix} = \begin{pmatrix} 4(1)+8(2)+12(3) & 4(1)+8(1)+12(1) \\ 4(1)+4(2)+4(3) & 4(1)+4(1)+4(1) \end{pmatrix} = \begin{pmatrix} 56 & 24 \\ 24 & 12 \end{pmatrix}\)
- Compute \(A^T W \mathbf{y}\). \(A^T W \mathbf{y} = \begin{pmatrix} 4 & 8 & 12 \\ 4 & 4 & 4 \end{pmatrix} \begin{pmatrix} 2.1 \\ 4.0 \\ 6.2 \end{pmatrix} = \begin{pmatrix} 4(2.1)+8(4.0)+12(6.2) \\ 4(2.1)+4(4.0)+4(6.2) \end{pmatrix} = \begin{pmatrix} 8.4+32+74.4 \\ 8.4+16+24.8 \end{pmatrix} = \begin{pmatrix} 114.8 \\ 49.2 \end{pmatrix}\)
- Solve the system \(\begin{pmatrix} 56 & 24 \\ 24 & 12 \end{pmatrix} \begin{pmatrix} v \\ x_0 \end{pmatrix} = \begin{pmatrix} 114.8 \\ 49.2 \end{pmatrix}\). The determinant is \(56(12) - 24(24) = 672 - 576 = 96\). The inverse is \(\frac{1}{96} \begin{pmatrix} 12 & -24 \\ -24 & 56 \end{pmatrix}\). \(\hat{\mathbf{p}} = \frac{1}{96} \begin{pmatrix} 12 & -24 \\ -24 & 56 \end{pmatrix} \begin{pmatrix} 114.8 \\ 49.2 \end{pmatrix} = \frac{1}{96} \begin{pmatrix} 12(114.8) - 24(49.2) \\ -24(114.8) + 56(49.2) \end{pmatrix}\) \(v = \frac{1377.6 - 1180.8}{96} = \frac{196.8}{96} = 2.05\ \mathrm{m/s}\) \(x_0 = \frac{-2755.2 + 2755.2}{96} = 0.00\ \mathrm{m}\)
- Calculate the uncertainty in \(v\). \(\mathrm{Cov}(\hat{\mathbf{p}}) = (A^T W A)^{-1} = \frac{1}{96} \begin{pmatrix} 12 & -24 \\ -24 & 56 \end{pmatrix}\) The variance of \(v\) is the first diagonal element: \(\sigma_v^2 = \frac{12}{96} = 0.125\ \mathrm{m^2/s^2}\). \(\sigma_v = \sqrt{0.125} \approx 0.354\ \mathrm{m/s}\). Likewise \(\sigma_{x_0} = \sqrt{56/96} \approx 0.76\ \mathrm{m}\), and the off-diagonal \(-24/96 = -0.25\ \mathrm{m^2/s}\) shows that slope and intercept are anticorrelated.
Units check. The design matrix \(A\) has units of seconds in the first column and dimensionless in the second. \(W\) has units of \(\mathrm{m}^{-2}\). \(A^T W A\) has units of \(\mathrm{s^2 m^{-2}}\) and \(\mathrm{m^{-2}}\). The inverse has units \(\mathrm{m^2 s^{-2}}\) (first diagonal element) and \(\mathrm{m^2}\) (second). The variance of \(v\) (first diagonal) has units \(\mathrm{m^2/s^2}\), so the standard deviation has units \(\mathrm{m/s}\). This is consistent.
Three quick tests
- What is the dimension of the matrix \(A^T W A\) if we have \(N\) data points and \(M\) parameters?
Answer
It is an \(M \times M\) matrix. - If the weight matrix \(W\) is the identity matrix, what does this imply about the measurement uncertainties?
Answer
It implies that all measurements have equal variance (and that variance is normalised to 1, or absorbed into the scale of \(\chi^2\)). - Why do we require the matrix \(A^T W A\) to be invertible?
Answer
To ensure a unique solution for the parameters \(\hat{\mathbf{p}}\). If it is singular, the parameters are not uniquely determined by the data (e.g., perfect collinearity).
Problem set for PU-106 →
Each derivation in this unit ends with five graded problems with full solutions.
Lecture notes drafted with a local model and checked line by line by the editor; every number in the worked example was recomputed.
Derivations cited: D-121 · Threads: chance