Normal equation systems ATAx=ATb have special characteristics that distinguish them from ordinary linear systems. Imagine finding the best point on a line to represent scattered data points, this system provides a mathematical way to find that optimal solution.
For a matrix with , the normal equation system always has a solution. More specifically, this system has a unique solution exactly when matrix has full rank, that is when . Under this condition, the solution can be expressed as .
A∈Rm×n
m≥n
ATAx=ATb
A
Rank(A)=n
x^=(ATA)−1ATb
When matrix A does not have full rank, the solution set of the normal equation system takes the form x^+Kernel(A), where x^ is any particular solution of the system.
The fundamental reason why normal equation systems always have a solution lies in the concept of orthogonal projection. The orthogonal projection of vector b onto the column space {Ax:x∈Rn} always exists and is the solution to the linear least squares problem, which automatically is also a solution to the normal equation system.
To understand why other solutions take the form x^+Kernel(A), suppose x~ is another solution of the system ATAx~=ATb. Then x~ is a solution to the normal equation system if and only if ATA(x~−x^)=0, which is equivalent to (x~−x^)TATA(x~−x^)=0, which is further equivalent to A(x~−x^)=0, or in other words x~−x^∈Kernel(A).
For a matrix A∈Rm×n with m≥n and Rank(A)=n, we can define the Moore-Penrose pseudoinverse as
A†=(ATA)−1AT
The Moore-Penrose pseudoinverse functions like the "best inverse" of a non-square matrix. It provides an optimal way to "cancel" linear transformations in the least squares context.
The Moore-Penrose pseudoinverse satisfies four Penrose axioms that uniquely determine its characteristics
AA†A=A
A†AA†=A†
(AA†)T=AA†
(A†A)T=A†A
These four properties are unique, meaning if a matrix B satisfies all four axioms, then automatically B=A†. The Moore-Penrose pseudoinverse thus functions as the unique solution operator for linear least squares problems.
A more numerically stable approach to solving normal equation systems uses QR decomposition. For a matrix A∈Rm×n with full rank and m≥n, we can use the (thin) QR decomposition A=Q1R1.
With this decomposition, the normal equation system ATAx=ATb can be solved through
ATAx=R1TQ1TQ1R1x=R1TR1x=R1TQ1Tb=ATb
Since R1 is an invertible upper triangular matrix, this equation is equivalent to
R1x=Q1Tb
This upper triangular system can be solved using back substitution, providing the solution x directly and efficiently.
Let's apply this method to a concrete example. Suppose we have experimental data that we want to fit with a quadratic polynomial. We will use the following data
A=9410149−3−2−101231111111
b=−2.2−4.2−4.2−1.81.88.215.8
Each row in matrix A has the format [ti2,ti,1] to find the coefficients of polynomial y=at2+bt+c, while vector b contains the corresponding observation values.
To solve normal equation systems, there are two main approaches that can be compared in terms of computation and numerical stability.
The Cholesky approach involves explicitly forming the matrix ATA first, then applying Cholesky decomposition since this matrix is positive definite. This method requires approximately n2⋅m+61n3+O(n2)+O(m⋅n) arithmetic operations. However, multiplication and decomposition can become sources of large error propagation, especially when m=n where cond(ATA)≈cond(A)2.
The QR approach, conversely, can solve this problem with better numerical stability and comparable computational complexity. The main complexity is determined by n2⋅m operations for QR decomposition, making it comparable to the Cholesky approach. However, the significant advantage of QR lies in the fact that orthogonal transformations do not worsen the problem condition, unlike forming ATA in the Cholesky method.
The choice of the appropriate method depends on the data characteristics and the level of accuracy required in the specific application.
