Sixty years ago James Wilkinson published his backward error analysis of Gaussian elimination for solving a linear system , where is a nonsingular matrix. He showed that in floating-point arithmetic the computed solution satisfies
where is the unit roundoff and is a low degree polynomial. The term is the growth factor, defined by
where the are the elements at the th stage of Gaussian elimination. The growth factor measures how much elements grow during the elimination. We would like the product to be of order 1, so that is a small relative perturbation of . We therefore need not to be too large.
With partial pivoting, in which row interchanges are used to ensure that at each stage the pivot element is the largest in its column, Wilkinson showed that and that equality is possible. Such exponential growth implies a large (unless we are lucky), meaning a severe loss of numerical stability. However, seventy years of digital computing experience have shown that is usually of modest size in practice. Explaining why this is the case is one of the outstanding problems in numerical analysis.
It is easy to experiment with growth factors in MATLAB. I will use the function
function g = gf(A) %GF Approximate growth factor. % g = GF(A) is an approximation to the % growth factor for LU factorization % with partial pivoting. [~,U] = lu(A); g = max(abs(U),,'all')/max(abs(A),,'all');
It computes a lower bound on the growth factor (since it only considers in the numerator in the definition), but it is entirely adequate for our purposes here. Let’s compute the growth factor for a random matrix of order 10,000 with elements from the standard normal distribution (mean 0, variance 1):
>> rng(1); n = 10000; gf(randn(n)) ans = 6.1335e+01
Growth of 61 is unremarkable for a matrix of this size. Now we try a matrix of the same size generated by the
>> A = gallery('randsvd',n,1e6,2,,,1); >> gf(A) ans = 9.7544e+02
This function generates an matrix with known singular value distribution and with singular vector matrices that are random orthogonal matrices from the Haar distribution. The parameter
1e6 specifies the 2-norm condition number, while the
2 (the mode parameter) specifies that there is only one small singular value, so the singular values are 1 repeated times and 1e-6. Growth of 975 is exceptional! These matrices have been in MATLAB since the 1990s, but this large growth property has apparently not been noticed before.
It turns out that mode 2 randsvd matrices generate with high probability growth factors of size at least for any condition number and for any pivoting strategy, not just partial pivoting. One way to check this is to randomly permute the columns of before doing the LU factorization with partial pivoting:
>> gf(A(:,randperm(n))) ans = 7.8395e+02
Here is a plot showing the maximum over 12 randsvd matrices for each of the growth factors for three different pivoting strategies, along with the maximum growth factors for partial pivoting for
randn matrices. The black curve is . This plot emphasizes the unusually large growth for mode 2 randsvd matrices.
What is the explanation for this large growth? It stems from three facts.
- Haar distributed orthogonal matrices have the property that that their elements are fairly small with high probability, as shown by Jiang in 2005.
- If the largest entries in magnitude of and are both small, in the sense that their product is , then will produce a growth factor of at least for any pivoting strategy. This was proved by Des Higham and I in the paper Large Growth Factors in Gaussian Elimination with Pivoting.
- If is an orthogonal matrix generating large growth then a rank-1 perturbation of 2-norm at most 1 tends to preserve the large growth.
For full details see the new EPrint Random Matrices Generating Large Growth in LU Factorization with Pivoting by Des Higham, Srikara Pranesh and me.
Is growth of order a problem in practice? It can be for two reasons.
- The largest dense linear systems solved today are of dimension . If we work in single precision then and so LU factorization can potentially be completely unstable if there is growth of order .
- For IEEE half precision arithmetic growth of order will cause overflow once exceeds . It was overflow in half precision LU factorization on randsvd matrices that alerted us to the large growth.