When do orthogonal transformations outperform Gaussian elimination?

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?..

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).

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.

• Upvoted, but QR should actually be on par with LU if you consider the unnecessary $n^3$ operations needed to form $A^H A$ (QR only requires $2/3 n^3$ more flops than LU). The SVD approach should still be slower though (one can think of its cost as roughly $6n^3$). – Jack Poulson Jan 29 '12 at 17:41
• In addition to the stability guaranteed by the use of orthogonal transformations, the great advantage of SVD is that the decomposition provides its own condition check, as the ratio of the largest to the smallest singular value is precisely the (2-norm) condition number. For the other decompositions, the use of a condition estimator (e.g. Hager-Higham) is, though not as expensive as the decomposition proper, somewhat "tacked on". – J. M. Jan 31 '12 at 0:25
• @JackPoulson Just out of curiosity, do you have a reference for your flop count for SVD? From what I can tell from a quick look in Golub & Van Loan (p. 254 3rd edition), the constant would seem higher for using the SVD in solving least-squares problems, but I could be mistaken. Thanks in advance. – OscarB Feb 8 '12 at 13:27
• @OscarB: It was a very rough number off the top of my head which is lower than forming the full SVD (because we can avoid backtransformation costs). $8/3 n^3$ work is needed for the reduction to bidiagonal form (say, $A=FBG^H$), some amount of work, say $C$, is needed for the bidiagonal SVD ($B=U\Sigma V^H$), and then $x := (G (V (\mathrm{inv}(\Sigma) (U^H (F^H b)))))$, which should require $O(n^2)$ work. Thus, it is all a matter of how big $C$ is...if MRRR ever works here it will be $O(n^2)$, but until then it is cubic and problem dependent. – Jack Poulson Feb 8 '12 at 15:36
• @J.M. Note, though, that the condition number of the least-squares problem is not the "classical" condition number $\frac{\sigma_1}{\sigma_n}$ of a matrix; it is a more complicated quantity. – Federico Poloni Mar 19 '19 at 18:57

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!).

• Any of the merits you've suggested are welcome. – faleichik Jan 27 '12 at 13:21

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.

• If the linear system is initially ill-conditioned then no matter which method you use: both LU and QR decomposition will give inaccurate results. QR can win only in cases when the process of Gaussian elimination "spoils" a good matrix. The main issue is that practical cases of such behavior are not known. – faleichik Jan 27 '12 at 14:48
• For most scientific applications, we generally obtain matrices that are sparse, symmetric, positive definite, and/or diagonally dominant. With very few exceptions, there is structure in the matrix that allows us to exploit certain techniques over traditional gaussian elimination. – Paul Jan 27 '12 at 15:46
• @Paul: On the other hand, dense Gaussian elimination is where most of the time is spent in the multifrontal method for sparse nonsymmetric matrices. – Jack Poulson Jan 27 '12 at 16:59
• @Paul It is just not true that "most applications produce SPD/diagonally dominant matrices". Yes, there is usually exploitable structure of some sort, but nonsymmetric and indefinite problems are extremely common. – Jed Brown Jan 27 '12 at 23:39
• "In fifty years of computing, no matrix problems that excite an explosive instability are known to have arisen under natural circumstances." - L.N. Trefethen and D. Bau They give an interesting probabilistic analysis in their book. – J. M. Jan 28 '12 at 2:44