Numerically stable solution methods

Factorization-based techniques make it possible to obtain the least-squares solution directly from the matrix $\mathbf{A}$, avoiding the explicit formation of $\mathbf{A}^{\top}\mathbf{A}$.

One possibility is to use QR factorization, an algorithm known to be numerically stable. For a matrix $\mathbf{A}$ of dimensions $m\times n$, with $m\geq n$, we can write

$$\mathbf{A}=\mathbf{Q}\mathbf{R}$$ (1.18)

where $\mathbf{Q}$ has orthonormal columns and $\mathbf{R}$ is an upper-triangular matrix. The least-squares problem therefore becomes

$$\mathbf{Q}\mathbf{R}\mathbf{x}\simeq\mathbf{y}.$$ (1.19)

Using the orthogonality of $\mathbf{Q}$, the solution can be obtained by solving

$$\mathbf{R}\mathbf{x} = \mathbf{Q}^{\top}\mathbf{y}.$$ (1.20)

This approach avoids forming the matrix $\mathbf{A}^{\top}\mathbf{A}$ and is therefore numerically more stable than directly solving the normal equations.

The SVD is an even more general method and also naturally handles rank-deficient matrices. Indeed, from decomposition (1.9) it is possible to obtain the pseudoinverse directly using (1.12) and hence the solution

$$\hat{\mathbf{x}} = \mathbf{V}\mathbf{\Sigma}^{+}\mathbf{U}^{\top}\mathbf{y}.$$ (1.21)

The SVD also makes it possible to identify directly the problem directions associated with very small singular values. These directions correspond to solution components that are weakly determined by the data and are particularly important in the analysis of conditioning and homogeneous linear systems.

From a numerical-computation perspective, whenever possible, the condition number of the matrix can also be reduced through suitable normalization or scaling of the columns of $\mathbf{A}$.