When do orthogonal transformations outperform Gaussian elimination?

Question

As we know, orthogonal transformations methods (Givens rotations and Housholder reflections) for systems of linear equations are more expensive than Gaussian elimination, but theoretically have nicer stability properties in the sense that they do not change the condition number of the system. Although I know just one academic example of a matrix which is spoiled by Gaussian elimination with partial pivoting. And there's common opinion that it is very unlikely to meet this kind of behavior in practice (see this lecture notes [pdf]).

So, where shall we look for the answer on the topic? Parallel implementations? Updating?..

Answer 1

Accuracy

Trefethen and Schreiber wrote an excellent paper, Average-case Stability of Gaussian Elimination, which discusses the accuracy side of your question. Here are a few of its conclusions:

1. "For QR factorization with or without column pivoting, the average maximal element of the residual matrix is $O(n^{1/2})$, whereas for Gaussian elimination it is $O(n)$. This comparison reveals that Gaussian elimination is mildly unstable, but the instability would only be detectable for very large matrix problems solved in low precision. For most practical problems, Gaussian elimination is highly stable on average." (Emphasis mine)

2. "After the first few steps of Gaussian elimination, the remaining matrix elements are approximately normally distributed, regardless of whether they started out that way."

There is much more to the paper that I can't capture here, including the discussion of the worst-case matrix you mentioned, so I strongly recommend that you read it.

Performance

For square real matrices, LU with partial pivoting requires roughly $2/3 n^3$ flops, whereas Householder-based QR requires roughly $4/3 n^3$ flops. Thus, for reasonably large square matrices, QR factorization will only be about twice as expensive as LU factorization.

For $m \times n$ matrices, where $m \ge n$, LU with partial pivoting requires $mn^2 - n^3/3$ flops, versus QR's $2mn^2 - 2n^3/3$ (which is still twice that of LU factorization). However, it is surprisingly common for applications to produce very tall skinny matrices ($m \gg n$), and Demmel et al. have a nice paper, Communication-avoiding parallel and sequential QR factorization, which (in section 4) discusses a clever algorithm which only requires $\log p$ messages to be sent when $p$ processors are used, versus the $n \log p$ messages of traditional approaches. The expense is that $O(n^3 \log p)$ extra flops are performed, but for very small $n$ this is often preferred to the latency cost of sending more messages (at least when only a single QR factorization needs to be performed).

Answer 2

I'm surprised no one mentioned linear least squares problems, which occur frequently in scientific computing. If you want to use Gaussian elimination, you have to form and solve the normal equations, which look like:

$$A^{T}Ax = A^{T}b,$$

where $A$ is a matrix of data points corresponding to observations of independent variables, $x$ is a vector of parameters to be found, and $b$ is a vector of data points corresponding to observations of a dependent variable.

As Jack Poulson frequently points out, the condition number of $A^{T}A$ is the square of the condition number of $A$, so the normal equations can be disastrously ill-conditioned. In such cases, although QR- and SVD-based approaches are slower, they yield much more accurate results.

Answer 3

How do you measure performance? Speed? Accuracy? Stability? A quick test in Matlab gives the following:

>> N = 100;
>> A = randn(N); b = randn(N,1);
>> tic, for k=1:10000, [L,U,p] = lu(A,'vector'); x = U\(L\b(p)); end; norm(A*x-b), toc
ans =
   1.4303e-13
Elapsed time is 2.232487 seconds.
>> tic, for k=1:10000, [Q,R] = qr(A); x = R\(Q'*b); end; norm(A*x-b), toc             
ans =
   5.0311e-14
Elapsed time is 7.563242 seconds.

So solving a single system with an LU-decomposition is about three times as fast as solving it with a QR-decomposition, at the cost of half a decimal digit of accuracy (this example!).

Answer 4

The article you cite defends Gaussian Elimination by saying that even though it is numerically unstable it tends to do well on random matrices and since most matrices one can think of are like random matrices, we should be ok. This same statement can be said of many numerically unstable methods.

Consider the space of all matrices. These methods work fine almost everywhere. That is 99.999...% of all matrices one could create will have no problems with unstable methods. There is only a very small fraction of matrices for which GE and others will have difficulty.

The problems that researchers care about tend to be in that small fraction.

We don't construct matrices randomly. We construct matrices with very special properties that correspond to very special, non-random systems. These matrices are often ill-conditioned.

Geometrically you can consider the linear space of all matrices. There is a zero volume/measure subspace of singular matrices cuts through this space. Many problems that we construct are clustered around this subspace. They are not distributed randomly.

As an example consider the heat equation or dispersion. These systems tend to remove information from the system (all initial states gravitate to a single final state) and as a result matrices that describe these equations are enormously singular. This process is very unlikely in a random situation yet ubiquitous in physical systems.