Normal equations

1 Least-squares setup and objective

1.1 Linear model and residual definition

In a standard least-squares problem, one seeks an unknown parameter vector \(x\) that explains observed data \(b\) through a linear model \[ b \approx A x, \] where \(A\) is a known matrix (often called the design matrix). The residual is the mismatch between the observations and the model prediction: \[ r(x) = b - A x. \]

1.2 Squared-error minimization

Least-squares estimation chooses \(x\) to minimize the squared Euclidean norm of the residual: \[

\min_x \|b - A x\|_2^2.

\] This objective penalizes deviations in proportion to their magnitude and produces a smooth function suitable for analytic optimization.

1.3 Geometry of residuals and orthogonality

The residual has a geometric meaning: for any candidate \(x\), the vector \(A x\) lies in the column space (range) of \(A\). The least-squares choice selects the point in this subspace that is closest to \(b\). As a consequence, the residual at the optimum is orthogonal to the column space of \(A\). This geometric orthogonality is the foundation of both the normal equations and their numerical implications.

2 Derivation of normal equations (stationarity)

2.1 Objective gradient and stationarity condition

Let \[

f(x)=\|b-Ax\|_2^2 = (b-Ax)^T(b-Ax)

\] for real vectors (the analogous complex form uses conjugate transposes). Expanding and differentiating yields a gradient that vanishes at the minimizer. Setting \(\nabla f(x)=0\) produces the stationarity condition that defines the normal equations.

2.2 From residual norm to matrix equation

2.2.1 Expanding the squared norm in matrix form

Compute: \[ f(x) = (b-Ax)^T(b-Ax)= b^T b - 2x^T A^T b + x^T A^T A x. \] Only the last two terms depend on \(x\).

2.2.2 Identifying the Gram matrix structure

Differentiation gives \[ \nabla f(x)= -2A^T b + 2A^T A x. \] The stationarity condition \(\nabla f(x)=0\) becomes \[ A^T A x = A^T b. \] The matrix \(A^T A\) is a Gram matrix: it encodes inner products between columns of \(A\) and is symmetric positive semidefinite.

2.3 Interpretation as first-order optimality

The normal equations represent the first-order optimality condition for the least-squares objective. In convex quadratic problems, any point satisfying stationarity corresponds to a global minimizer, though uniqueness depends on rank.

3 Projection interpretation

3.1 Orthogonal projection onto the column space

Let \(\mathcal{C}(A)\) denote the column space of \(A\). The least-squares fit corresponds to projecting \(b\) orthogonally onto \(\mathcal{C}(A)\). If \(\hat{b}\) is the projected vector, then \(\hat{b}\in \mathcal{C}(A)\) and \[ \hat{b} = A\hat{x} \] for some \(\hat{x}\) (not necessarily unique in rank-deficient settings).

3.2 Fitted values as projected vectors

The fitted values \(A\hat{x}\) are exactly this closest vector in the column space. The residual is \[ \hat{r} = b - \hat{b}. \] Thus the optimization is equivalent to selecting \(\hat{b}\) as the orthogonal projection of \(b\).

3.3 Residual orthogonality to the column space

3.3.1 Pythagorean relationship in least squares

Orthogonality implies \[

\|b\|_2^2 = \|\hat{b}\|_2^2 + \|\hat{r}\|_2^2,

\] a direct form of the Pythagorean theorem in Euclidean space. More importantly for computation, orthogonality means \(\hat{r}\) has zero inner product with every vector in \(\mathcal{C}(A)\), which yields \[ A^T \hat{r} = A^T(b-A\hat{x})=0 \] and therefore the normal equations.

4 Solving the normal equations

4.1 Full-rank case and uniqueness

If \(A\) has full column rank, then \(A^T A\) is positive definite and invertible. The normal equations have a unique solution: \[ \hat{x} = (A^T A)^{-1}A^T b. \] In this case, both existence and uniqueness are guaranteed.

4.2 Rank-deficient case and non-uniqueness

When \(A\) is rank-deficient, \(A^T A\) becomes singular. The set of solutions to \[ A^T A x = A^T b \] may contain infinitely many vectors \(x\), all producing the same fitted values \(A x\), because they differ only within the null space of \(A\).

4.3 Conditions for existence of solutions

4.3.1 Consistency and solvability via subspace inclusion

The normal equations are consistent precisely when \(A^T b\) lies in the column space of \(A^T A\). Because \(A^T A\) shares its range with \(A^T\), least-squares always admits solutions in the original minimization problem. At the level of equations, this translates to solvability conditions tied to subspace inclusion relationships among the ranges of \(A\) and \(A^T\).

4.4 Computational approaches and stability concerns

4.4.1 Forming the Gram matrix vs. direct methods

A common approach forms \(A^T A\) explicitly and then solves the resulting linear system. This can be numerically problematic because it squares the condition number: small singular values of \(A\) become even smaller in \(A^T A\). For improved numerical stability, direct methods such as QR factorization or singular value decomposition (SVD) often avoid this amplification, while still producing the same least-squares fit.

5 Relationship to pseudoinverses and minimal-norm solutions

5.1 Moore–Penrose pseudoinverse overview

The Moore–Penrose pseudoinverse \(A^+\) generalizes the matrix inverse to singular or rectangular matrices. It yields a canonical way to map data \(b\) to parameter estimates even when \(A^T A\) is not invertible.

5.2 Minimum-norm solution characterization

In rank-deficient settings, among all vectors \(x\) that satisfy the least-squares objective, there is typically a preferred one: the solution with smallest Euclidean norm \(\|x\|_2\). This minimal-norm solution is given by

\[ \hat{x} = A^+ b. \] It selects a specific member from the affine set of minimizers.

5.3 Connection to least-squares via projections

The pseudoinverse characterization ties back to the geometric view: \(A A^+ b\) equals the orthogonal projection of \(b\) onto \(\mathcal{C}(A)\). Therefore, using \(A^+\) produces the same projected fitted vector \(\hat{b}\), while the minimal-norm property resolves ambiguity in \(x\).

6 Alternative derivations and equivalent formulations

6.1 Using inner products and orthogonality conditions

The orthogonality condition can be derived directly without differentiating. If \(\hat{b}=A\hat{x}\) is the projection of \(b\) onto \(\mathcal{C}(A)\), then for every column direction \(A y\) one has \[ (b-\hat{b})^T(Ay)=0. \] Rewriting gives \(A^T(b-A\hat{x})=0\), which again becomes the normal equations.

6.2 Lagrange-multiplier perspective (constraint form)

One can introduce auxiliary variables, for example letting \(y=Ax\) and constraining \(y\) to lie in \(\mathcal{C}(A)\). Minimizing \(\|b-y\|_2^2\) subject to \(y=Ax\) leads, via Lagrange multipliers, to the same optimality condition \(A^T(b-Ax)=0\).

6.3 Equivalent normal equations from residual minimization

All these routes—gradient stationarity, orthogonality, and constrained minimization—converge to the same linear system. The normal equations therefore serve as an equivalence hub: they are not merely an algebraic convenience but express the optimality geometry in coordinate form.

7 Numerical linear algebra considerations

7.1 Conditioning and sensitivity to scaling

Because \(A^T A\) depends quadratically on the entries of \(A\), the resulting system can be sensitive to scaling and to the magnitude range of columns. Poor scaling can worsen numerical errors and increase the risk of inaccurate solutions, especially when \(A\) is ill-conditioned.

7.2 Cholesky factorization for symmetric positive semidefinite systems

When \(A^T A\) is well-conditioned and symmetric positive definite (for full column rank), it can be factored using Cholesky decomposition: \[ A^T A = R^T R. \] For the semidefinite case, Cholesky may require modifications or pivoting strategies. Factorization-based solvers are often efficient, but they still inherit conditioning effects from forming \(A^T A\).

7.3 Iterative methods for least-squares stationarity

7.3.1 Conjugate gradients on the normal equations

Iterative solvers such as conjugate gradients can be applied to symmetric positive definite systems. Applied to \(A^T A x = A^T b\), conjugate gradients exploits the symmetry and can be effective for large problems. However, the same conditioning concerns apply: ill-conditioning in \(A^T A\) can slow convergence, so preconditioning is often important in practice.

8 Extensions and variants

8.1 Weighted least squares and modified normal equations

In weighted least squares, the objective becomes \[

\min_x \|W^{1/2}(b-Ax)\|_2^2

\] for a positive semidefinite weight matrix \(W\). The stationarity condition becomes a modified normal equation: \[ A^T W A x = A^T W b. \] This reweights residual directions according to the chosen metric.

8.2 Regularized least squares (ridge/Tikhonov) connections

Regularization modifies the objective by adding a penalty on the parameter size: \[

\min_x \|b-Ax\|_2^2 + \lambda \|x\|_2^2,

\quad \lambda>0. \] The stationarity system becomes \[ (A^T A + \lambda I)x = A^T b, \] which is symmetric positive definite for \(\lambda>0\), improving numerical behavior and mitigating overfitting in statistical interpretations.

8.2.1 How regularization alters the stationarity system

By adding \(\lambda I\), regularization shifts eigenvalues of \(A^T A\) upward, reducing sensitivity to noise and stabilizing inversion-like steps. The fit is no longer the pure projection associated with the unregularized least-squares objective, but it remains closely related to projection in an altered geometry.

8.3 Total least squares distinction (overview-level contrast)

Total least squares (TLS) differs in that it models errors in both \(A\) and \(b\), whereas standard least squares assumes \(A\) is exact and only residuals in \(b\) are minimized. Consequently, TLS does not lead to the same normal equations derived from \(\|b-Ax\|_2^2\); it instead relies on singular vector structure of the augmented data.

9 Summary and key takeaways

9.1 Core equivalences: stationarity, normal equations, projection

Normal equations are the coordinate form of least-squares optimality. They arise from setting the gradient of the squared residual norm to zero, and they encode the geometric fact that the residual at the optimum is orthogonal to the column space of the design matrix. This yields a projection interpretation of the fitted response.

9.2 Practical guidance: when and how to use normal equations

Normal equations can be used to compute solutions efficiently when the system is well-conditioned and \(A^T A\) can be handled reliably. For ill-conditioned or poorly scaled problems, methods that avoid forming \(A^T A\) (such as QR or SVD) often produce more accurate results. When rank deficiency or ambiguity is present, pseudoinverse-based solutions provide a principled minimal-norm estimate.