What Is the Growth Factor for Gaussian Elimination?

Gaussian elimination is the process of reducing an n \times n matrix to upper triangular form by elementary row operations. It consists of n stages, in the kth of which multiples of row k are added to later rows to eliminate elements below the diagonal in the kth column. The result of Gaussian elimination (assuming it succeeds) is a factorization A = LU, where L is unit lower triangular (lower triangular with ones on the diagonal) and U is upper triangular. This factorization reduces the solution of a linear system Ax = b to the solution of two triangular systems.

James Wilkinson showed in the early 1960s that the numerical stability of Gaussian elimination depends on the growth factor

\rho_n = \displaystyle\frac{\max_{i,j,k} |a_{ij}^{(k)}|} {\max_{i,j}|a_{ij}|} \ge 1,

where the a_{ij}^{(k)} are the elements at the start of the kth stage of Gaussian elimination, with A = \bigl(a_{ij}^{(1)}\bigr). Specifically, he obtained a bound for the backward error of the computed solution that is proportional to \rho_nu, where u is the unit roundoff of the floating-point arithmetic. Wilkinson’s analysis focused attention on the size of the growth factor.

If no row interchanges are done during the elimination then \rho_n can be arbitrarily large, as is easily seen for n = 2. For some specific types of matrix more can be said.

  • If A is symmetric positive definite then \rho_n = 1. In this case one would normally use Cholesky factorization instead of LU factorization.
  • If A is nonsingular and totally nonnegative (every submatrix has nonnegative determinant) then \rho_n = 1.
  • If A is diagonally dominant by rows or columns then \rho_n \le 2.
  • If A is complex and its real and imaginary parts are both symmetric positive definite then \rho_n < 3.

The entry in the (k,k) position at the start of stage k of Gaussian elimination is called the pivot. By swapping rows and columns, any element with row and column indices greater than or equal to k can be brought into the pivot position and used as the pivot. A pivoting strategy is a strategy for interchanging rows and columns at the start of each stage in order to choose a good sequence of pivots. For dense matrices the main criterion is that the pivoting strategy keeps the growth factor small, while for sparse matrices minimizing fill-in (a zero entry becoming nonzero) is also important.

Partial Pivoting

Partial pivoting interchanges rows to ensure that the pivot element is the largest in magnitude in its column. Wilkinson showed that partial pivoting ensures \rho_n \le 2^{n-1} and that equality is attained for matrices of the form illustrated for n = 4 by

\left[\begin{array}{@{\mskip3mu}rrrr} 1 & 0 & 0 & 1 \\ -1 & 1 & 0 & 1 \\ -1 & -1 & 1 & 1 \\ -1 & -1 & -1 & 1 \\ \end{array} \right].

This matrix is a special case of a larger class of matrices for which equality is attained (Higham and Higham, 1989). Exponential growth can occur for matrices arising in practical applications, as shown by Wright for the multiple shooting method for two-point boundary value problems and Foster for a quadrature method for solving Volterra integral equations. Remarkably, though, \rho_n is usually of modest size in practice. Explaining this phenomenon is one of the outstanding problems in numerical analysis.

Rook Pivoting

Rook pivoting interchanges rows and columns to ensure that the pivot element is the largest in magnitude in both its row and its column. Foster has shown that

\rho_n \le 1.5 n^{\frac{3}{4} \log n}.

This bound grows much more slowly than the bound for partial pivoting, but it can still be large. In practice, the bound is pessimistic.

Complete Pivoting

Complete pivoting chooses as pivot the largest element in magnitude in the remaining part of the matrix. Wilkinson (1961) showed that

\rho_n \le n^{1/2} (2\cdot 3^{1/2} \cdots n^{1/(n-1)})^{1/2} \sim c \, n^{1/2} n^{\frac{1}{4}\log n}.

This bound grows much more slowly than that for rook pivoting, but again it is pessimistic. Wilkinson noted in 1965 that no matrix had been discovered for which \rho_n > n. Much interest has focused on the question of whether growth of n can be exceeded and what is the largest possible growth. For a Hadamard matrix, \rho_n \ge n, but such matrices do not exist for all n. Evidence suggests that the growth factor for complete pivoting on a Hadamard matrix is exactly n, and this has been proved for n = 12 and n = 16.

That \rho_n > n is possible was shown in the early 1990s by Gould (in floating-point arithmetic) and Edelman (in exact point arithmetic), through a 13\times 13 matrix with growth 13.02.

Growth Factor Bounds for Any Pivoting Strategy

While much attention has focused on investigating the growth factor for particular pivoting strategies, some results are known that provide growth factor lower bounds for any pivoting strategy. Higham and Higham (1989) obtained a result that implies that for the symmetric orthogonal matrix

S = \sqrt{ \displaystyle\frac{2}{n+1} } \biggl( \sin\Bigl( \displaystyle\frac{ij\pi}{n+1} \Bigr) \biggr)_{i,j=1}^n

the lower bound

\rho_n(S)\ge \displaystyle\frac{n+1}{2}

holds for any pivoting strategy. Recently, Higham, Higham, and Pranesh (2020) have shown that

\rho_n(Q) \gtrsim \displaystyle\frac{n}{4\log n}

with high probability for large n when Q is a random orthogonal matrix from the Haar distribution. They show that the same bound holds for matrices that have only one singular value different from 1 and have singular vector matrices that are random orthogonal matrices from the Haar distribution. These “randsvd matrices” can be generated as gallery('randsvd',n,kappa,2) in MATLAB.

Experimentation

If you want to experiment with growth factors in MATLAB, you can use the function gep from the Matrix Computation Toolbox. Examples:

>> rng(0); A = randn(1000);
>> [L,U,P,Q,rho] = gep(A,'p'); rho  % Partial pivoting.
rho =
   1.9437e+01
>> [L,U,P,Q,rho] = gep(A,'r'); rho  % Rook pivoting.
rho =
   9.9493e+00
>> [L,U,P,Q,rho] = gep(A,'c'); rho  % Complete pivoting.
rho =
   6.8363e+00
>> A = gallery('randsvd',1000,1e1,2);
>> [L,U,P,Q,rho] = gep(A,'c'); rho  % Complete pivoting.
rho =
   5.4720e+01

References

This is a minimal set of references, which contain further useful references within.

Related Blog Posts

This article is part of the “What Is” series, available from https://nhigham.com/category/what-is and in PDF form from the GitHub repository https://github.com/higham/what-is.

3 thoughts on “What Is the Growth Factor for Gaussian Elimination?”

  1. Hi Sirs.
    There must be some bug in the “gep” function when rook pivoting is chosen.
    In the example
    rng(0); A = randn(1000);
    [L,U,P,Q,rho] = gep(A,’r’); rho % Rook pivoting.
    a warning “too many FOR loop iterations” is issued.

    Reply

Leave a comment