📂Matrix Algebra

Stability Comparison of LU, QR, and SVD in the Least Squares Method

Summary

For a matrix $A \in \mathbb{R}^{m \times n}$, suppose $m \ge n$. When performing the least squares method via LU decomposition, QR decomposition, and singular value decomposition, all three methods are backward stable. However, the conditions under which backward stability holds differ. Saying that a matrix is full rank means that the rank $r$ of $A$ equals $\min (m, n)$.

• LU decomposition¹: it is unstable even if $A$ is full rank.
• QR decomposition²: it is backward stable if $A$ is full rank.
• Singular value decomposition³: it is backward stable even if $A$ is not full rank.

Explanation

There are various ways to use matrix decomposition as a practical implementation of the least squares method, but three in particular are especially widely used: the method via LU decomposition (Cholesky decomposition), the method via QR decomposition, and the method via SVD decomposition. For reference, whether it is LU decomposition, Cholesky decomposition, or Gaussian elimination, they are essentially methods of solving a linear system through forward substitution and back substitution, so from the standpoint of implementing the least squares method they may all be regarded as the same thing. In fact, the flops of Cholesky decomposition is $mn^{2} + {\frac{ 1 }{ 3 }} n^{3}$, which is a bit better than LU decomposition, but in practice one ends up using PLU decomposition, which includes pivoting, for the sake of stability.

Looking at the implementations of Julia⁴ and MATLAB⁵, one can essentially regard it as using QR decomposition when $m \ne n$ in $A$, and LU decomposition when $m = n$.

LU Decomposition

If we consider the condition number of the matrix $\kappa = \kappa(A)$, the upper bound of the algorithm’s relative error depends on $\kappa$; and whether it is LU decomposition or Cholesky decomposition, at the moment the normal equation $A^{\top} A \mathbf{x} = A^{\top} \mathbf{b}$ is set up, one already has $\kappa \left( A^{\top} A \right) = \kappa^{2}$, so the condition number, which may already be poor, worsens further.

For $A = LU$, we call $\rho = {\frac{ \max \left| U_{ij} \right| }{ \max \left| A_{ij} \right| }}$ the growth factor of $A$; at best, in pivoting LU one has $\rho = O (1)$, and it has backward stability when the relative error of $A$ is $O \left( \rho \epsilon \right)$ with respect to $\rho$ and the machine epsilon $\epsilon$.

As you can probably sense at a glance, it is backward stable in the end, but there are so many long-winded conditions required that it is effectively unstable. This is because, when performing the elementary row operations of Gaussian elimination, the errors arising from multiplying and adding numbers can grow larger and larger.

QR Decomposition

The process of QR decomposition:

• Step 2-4. Compute the following. $$ \mathbf{q}_{j} = {{ \mathbf{v}_{j} } \over {r_{jj} }} $$

QR decomposition too has various ways of being implemented, so we cannot say that it loses stability strictly for this reason; but at least when using Gram-Schmidt orthogonalization as above, problems can arise in processes such as dividing a vector by $r_{jj}$. It literally depends on whether $A$ is ill-conditioned or not; in other words, on whether such singularity lurks within the problem.

However, it is nothing compared to LU, and since it is the most general-purpose one in use, there is no need to worry excessively about stability. For reference, in QR decomposition too, pivoting is usually included by default.

SVD

The only algorithm that is stable even in the presence of rank-deficiency is SVD. Mathematically, when computing the inverse matrix $\Sigma^{-1}$ of $\Sigma$ in $A = U \Sigma V^{\top}$, one can attempt a reduced SVD that entirely removes the parts of $A$’s singular values that are close to $0$. Unlike other decomposition methods that must rely on pivoting, this is tantamount to explicitly setting a sufficiently small threshold to fundamentally block numerical problems at their source. Because its basic cost is high, it is not commonly used for implementing the least squares method, but by that same measure it never falls behind in stability.

Trade-off with Speed

If we compare the speed of the three methods and compare the stability, we can very clearly confirm a trade-off as follows.

Method	Speed	Stability
LU	Fast	Unstable
QR	Moderate	Stable
SVD	Slow	Very stable