Accuracy and Stability of Numerical Algorithms

The front matter includes the title page, copyright page, TOC, list of figures, list of tables, preface to the second edition, preface to the first edition, about the dedication, and dedication.

Numerical precision is the very soul of science.

— SIR D'ARCY WENTWORTH THOMPSON, On Growth and Form (1942)

There will always be a small but steady demand for error-analysts to … expose bad algorithms' big errors and, more important, supplant bad algorithms with provably good ones.

— WILLIAM M. KAHAN, Interval Arithmetic Options in the Proposed IEEE Floating Point Arithmetic Standard (1980)

Since none of the numbers which we take out from logarithmic and trigonometric tables admit of absolute precision, but are all to a certain extent approximate only, the results of all calculations performed by the aid of these numbers can only be approximately true … It may happen, that in special cases the effect of the errors of the tables is so augmented that we may be obliged to reject a method, otherwise the best, and substitute another in its place.

— CARL FRIEDRICH GAUSS, Theoria Motus (1809)

Backward error analysis is no panacea; it may explain errors but not excuse them.

— HEWLETT-PACKARD, HP-15C Advanced Functions Handbook (1982)

This book is concerned with the effects of finite precision arithmetic on numerical algorithms, particularly those in numerical linear algebra. Central to any understanding of high-level algorithms is an appreciation of the basic concepts of finite precision arithmetic. This opening chapter briskly imparts the necessary background material. Various examples are used for illustration, some of them familiar (such as the quadratic equation) but several less well known. Common misconceptions and myths exposed during the chapter are highlighted towards the end, in §1.19.

This chapter has few prerequisites and few assumptions are made about the nature of the finite precision arithmetic (for example, the base, number of digits, or mode of rounding, or even whether it is floating point arithmetic). The second chapter deals in detail with the specifics of floating point arithmetic.

A word of warning: some of the examples from §1.12 onward are special ones chosen to illustrate particular phenomena. You may never see in practice the extremes of behaviour shown here. Let the examples show you what can happen, but do not let them destroy your confidence in finite precision arithmetic!

1.1. Notation and Background

We describe the notation used in the book and briefly set up definitions needed for this chapter.

Generally, we use

 capital letters A, B, C, Δ, Λ for matrices,

 subscripted lower case letters a_ij, b_ij, c_ij, δ_ij, λ_ij for matrix elements,

 lower case letters x, y, z, c, g, h for vectors,

 lower case Greek letters α, β, γ, θ,π for scalars,

following the widely used convention originally introduced by Householder [644, 1964].

The vector space of all real m × n matrices is denoted by ℝ^{m × n} and the vector space of real n-vectors by ℝⁿ. Similarly, ℂ^{m × n} denotes the vector space of complex m × n matrices. A superscript “T” denotes transpose and a superscript “*” conjugate transpose.

Algorithms are expressed using a pseudocode based on the MATLAB language [576, 2000], [824]. Comments begin with the % symbol.

Submatrices are specified with the colon notation, as used in MATLAB and Fortran 90/95: A(p: q, r: s) denotes the submatrix of A formed by the intersection of rows p to q and columns r to s. As a special case, a lone colon as the row or column specifier means to take all entries in that row or column; thus A(:, j) is the jth column of A, and A(i, :) the ith row. The values taken by an integer variable are also described using the colon notation: “i = 1: n” means the same as “i = 1, 2, …, n”.

Evaluation of an expression in floating point arithmetic is denoted fl(·), and we assume that the basic arithmetic operations op = +, −, *, / satisfy fl (x op y)= (x op y) (1+δ) , |δ| ≤u. 1.1

From 1946–1948 a great deal of quite detailed coding was done. The subroutines for floating-point arithmetic were … produced by Alway and myself in 1947 … They were almost certainly the earliest floating-point subroutines.

— J. H. WILKINSON, Turing's Work at the National Physical Laboratory … (1980)

MATLAB's creator Cleve Moler used to advise foreign visitors not to miss the country's two most awesome spectacles: the Grand Canyon, and meetings of IEEE p754.

— MICHAEL L. OVERTON, Numerical Computing with IEEE Floating Point Arithmetic (2001)

Arithmetic on Cray computers is interesting because it is driven by a motivation for the highest possible floating-point performance … Addition on Cray computers does not have a guard digit, and multiplication is even less accurate than addition … At least Cray computers serve to keep numerical analysts on their toes!

— DAVID GOLDBERG, Computer Arithmetic (1996)

It is rather conventional to obtain a “realistic” estimate of the possible overall error due to k roundoffs, when k is fairly large, by replacing k by k in an expression for (or an estimate of) the maximum resultant error.

— F. B. HILDEBRAND, Introduction to Numerical Analysis (1974)

2.1. Floating Point Number System

A floating point number system F ⊂ ℝ is a subset of the real numbers whose elements have the form y=±m× βe−t . 2.1 The system F is characterized by four integer parameters:

• the base β (sometimes called the radix),

• the precision t, and

• the exponent range e_min ≤ e ≤ e_max.

The significand m is an integer satisfying 0 ≤ m ≤ β^t − 1. To ensure a unique representation for each nonzero y ∈ F it is assumed that m ≥ β^t−1 if y ≠ 0, so that the system is normalized. The number 0 is a special case in that it does not have a normalized representation. The range of the nonzero floating point numbers in F is given by β^{e_min − 1} ≤ |y| ≤ β^e_max(l − β^−t). Values of the parameters for some machines of historical interest are given in Table 2.1 (the unit roundoff u is defined on page 38).

Note that an alternative (and more common) way of expressing y is y=± βe ( d1 β + d2 β2 +⋯+ dt βt ) =± βe ×. d1 d2 … dt , 2.2 where each digit d_i satisfies 0 ≤ d_i ≤ β − 1, and d₁ ≠ 0 for normalized numbers. We prefer the more concise representation (2.1), which we usually find easier to work with. This “nonpositional” representation has pedagogical advantages, being entirely integer based and therefore simpler to grasp. In the representation (2.2), d₁ is called the most significant digit and d_t the least significant digit.

It is important to realize that the floating point numbers are not equally spaced. If β = 2, t = 3, e_min = −1, and e_max = 3 then the nonnegative floating point numbers are 0,0.25,0.3125,0.3750,0.4375,0.5,0.625,0.750,0.875,1.0,1.25,1.50,1.75,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0.

A method of inverting the problem of round-off error is proposed which we plan to employ in other contexts and which suggests that it may be unwise to separate the estimation of round-off error from that due to observation and truncation.

— WALLACE J. GIVENS, Numerical Computation of the Characteristic Values of a Real Symmetric Matrix (1954)

The enjoyment of one's tools is an essential ingredient of successful work.

— DONALD E. KNUTH, The Art of Computer Programming, Volume 2, Seminumerical Algorithms (1998)

The subject of propagation of rounding error, while of undisputed importance in numerical analysis, is notorious for the difficulties which it presents when it is to be taught in the classroom in such a manner that the student is neither insulted by lack of mathematical content nor bored by lack of transparence and clarity.

— PETER HENRICI, A Model for the Propagation of Rounding Error in Floating Arithmetic (1980)

The two main classes of rounding error analysis are not, as my audience might imagine, ‘backwards’ and ‘forwards’, but rather ‘one's own’ and ‘other people's’. One's own is, of course, a model of lucidity; that of others serves only to obscure the essential simplicity of the matter in hand.

— J. H. WILKINSON, The State of the Art in Error Analysis (1985)

Having defined a model for floating point arithmetic in the last chapter, we now apply the model to some basic matrix computations, beginning with inner products. This first application is simple enough to permit a short analysis, yet rich enough to illustrate the ideas of forward and backward error. It also raises the thorny question of what is the best notation to use in an error analysis. We introduce the “γ_n” notation, which we use widely, though not exclusively, in the book. The inner product analysis leads immediately to results for matrix-vector and matrix-matrix multiplication.

Also in this chapter we determine a model for rounding errors in complex arithmetic, derive some miscellaneous results of use in later chapters, and give a high-level overview of forward error analysis and backward error analysis.

3.1. Inner and Outer Products

Consider the inner product s_n = x^Ty, where x, y ∈ ℝⁿ. Since the order of evaluation of s_n = x₁y₁ + ⋯ + x_ny_n affects the analysis (but not the final error bounds), we will assume that the evaluation is from left to right. (The effect of particular orderings is discussed in detail in Chapter 4, which considers the special case of summation.) In the following analysis, and throughout the book, a hat denotes a computed quantity.

Let s_i = x₁y₁ + ⋯ + x_iy_i denote the ith partial sum. Using the standard model (2.4), we have s ̑ 1 =fl ( x1 y1 )= x1 y1 (1+ δ1 ) , s ̑ 2 =fl ( s ̑ 1 + x2 y2 )= ( s ̑ 1 + x2 y2 (1+ δ2 ) )(1+ δ3 ) = ( x1 y1 (1+ δ1 )+ x2 y2 (1+ δ2 ) )(1+ δ3 ) = x1 y1 (1+ δ1 ) (1+ δ3 )+ x2 y2 (1+ δ2 ) (1+ δ3 ) , 3.1 where |δ_i| ≤ u, i = 1: 3. For our purposes it is not necessary to distinguish between the different δ_i terms, so to simplify the expressions let us drop the subscripts on the δ_i and write 1 + δ_i ≡ 1 ± δ. Then s ̑ 3 =fl ( s ̑ 2 + x3 y3 )= ( s ̑ 2 + x3 y3 (1±δ) )(1±δ) = ( x1 y1 (1±δ)2 + x2 y2 (1±δ)2 + x3 y3 (1±δ) )(1±δ) = x1 y1 (1±δ)3 + x2 y2 (1±δ)3 + x3 y3 (1±δ)2 . The pattern is clear. Overall, we have s ̑ n = x1 y1 (1±δ)n + x2 y2 (1±δ)n + x3 y3 (1±δ)n−1 +⋯+ xn yn (1±δ)2 . 3.2 There are various ways to simplify this expression. A particularly elegant way is to use the following result.

I do hate sums. There is no greater mistake than to call arithmetic an exact science. There are … hidden laws of Number which it requires a mind like mine to perceive. For instance, if you add a sum from the bottom up, and then again from the top down, the result is always different.

— MRS. LA TOUCHE

Joseph Fourier introduced this delimited ∑-notation in 1820, and it soon took the mathematical world by storm.

— RONALD L. GRAHAM, DONALD E. KNUTH, and OREN PATASHNIK, Concrete Mathematics (1994)

One of the major difficulties in a practical [error] analysis is that of description. An ounce of analysis follows a pound of preparation.

— BERESFORD N. PARLETT, Matrix Eigenvalue Problems (1965)

Sums of floating point numbers are ubiquitous in scientific computing. They occur when evaluating inner products, means, variances, norms, and all kinds of nonlinear functions. Although at first sight summation might appear to offer little scope for algorithmic ingenuity, the usual “recursive summation” (with various orderings) is just one of a variety of possible techniques. We describe several summation methods and their error analyses in this chapter. No one method is uniformly more accurate than the others, but some guidelines can be given on the choice of method in particular cases.

4.1. Summation Methods

In most circumstances in scientific computing we would naturally translate a sum ∑ i=1 n xi into code of the form s=0for i=1:ns=s+ xi end This is known as recursive summation. Since the individual rounding errors depend on the operands being summed, the accuracy of the computed sum s̑ varies with the ordering of the x_i. (Hence Mrs. La Touche, quoted at the beginning of the chapter, was correct if we interpret her remarks as applying to floating point arithmetic.) Two interesting choices of ordering are the increasing order |x₁| ≤ |x₂| ≤ ⋯ ≤ |x_n|, and the decreasing order |x₁| ≥ |x₂| ≥ ⋯ ≥ |x_n|.

Another method is pairwise summation (also known as cascade summation, or fan-in summation), in which the x_i are summed in pairs according to yi = x2i−1 + x2i, i=1: ⌊ n 2 ⌋ ( y (n+1) /2 = xn if n is odd), and this pairwise summation process is repeated recursively on the y_i, i = 1: ⌊(n + 1)/2⌋. The sum is obtained in ⌈log₂ n⌉ stages. For n = 6, for example, pairwise summation forms S6 = ( ( x1 + x2 ) +( x3 + x4 ))+( x5 + x6 ). Pairwise summation is attractive for parallel computing, because each of the ⌈log₂ n⌉ stages can be done in parallel [629, 1988, §5.2.2].

A third summation method is the insertion method. First, the x_i are sorted by order of increasing magnitude (alternatively, some other ordering could be used). Then x₁ + x₂ is formed, and the sum is inserted into the list x₂, … , x_n, maintaining the increasing order. The process is repeated recursively until the final sum is obtained. In particular cases the insertion method reduces to one of the other two. For example, if x_i = 2^{i − 1}, the insertion method is equivalent to recursive summation, since the insertion is always to the bottom of the list: 1248→ 3 ̲ 48→ 7 ̲ 8→15.

The polynomial (z − 1)(z − 2) … (z − 20) is not a ‘difficult’ polynomial per se … The ‘difficulty’ with the polynomial Π(z − i) is that of evaluating the explicit polynomial accurately. If one already knows the roots, then the polynomial can be evaluated without any loss of accuracy.

— J. H. WILKINSON, The Perfidious Polynomial (1984)

I first used backward error analysis in connection with simple programs for computing zeros of polynomials soon after the PILOT ACE came into use.

— J. H. WILKINSON, The State of the Art in Error Analysis (1985)

The Fundamental Theorem of Algebra asserts that every polynomial equation over the complex field has a root. It is almost beneath the dignity of such a majestic theorem to mention that in fact it has precisely n roots.

— J. H. WILKINSON, The Perfidious Polynomial (1984)

It can happen … that a particular polynomial can be evaluated accurately by nested multiplication, whereas evaluating the same polynomial by an economical method may produce relatively inaccurate results.

— C. T. FIKE, Computer Evaluation of Mathematical Functions (1968)

Two common tasks associated with polynomials are evaluation and interpolation: given the polynomial find its values at certain arguments, and given the values at certain arguments find the polynomial. We consider Horner's rule for evaluation and the Newton divided difference polynomial for interpolation. A third task not considered here is finding the zeros of a polynomial. Much research was devoted to polynomial zero finding up until the late 1960s; indeed, Wilkinson devotes a quarter of Rounding Errors in Algebraic Processes [1232, 1963] to the topic. Since the development of the QR algorithm for finding matrix eigenvalues there has been less demand for polynomial zero finding, since the problem either arises as, or can be converted to (see §28.6 and [383, 1995], [1144, 1994]), the matrix eigenvalue problem.

5.1. Horner's Method

The standard method for evaluating a polynomial p (x) = a0 + a1 x+⋯+ an xn 5.1 is Horner's method (also known as Horner's rule and nested multiplication), which consists of the following recurrence: qn (x)= an for i=n−1:−1:0 qi (x)=x qi+1 (x)+ ai endp (x) = q0 (x) The cost is 2n flops, which is n less than the more obvious method of evaluation that explicitly forms powers of x (see Problem 5.2).

To analyse the rounding errors in Horner's method it is convenient to use the relative error counter notation <k> (see (3.10)). We have q ̑ n−1 = (x q ̑ n <1>+ an−1 )<1> = x an <2>+ an−1 <1> , q ̑ n−2 = (x q ̑ n−1 <1>+ an−2 )<1> = x2 an <4>+ xan−1 <3>+ an−2 <1>. It is easy to either guess or prove by induction that q ̑ 0 = a0 <1>+ a1 x<3>+⋯+ an−1 xn−1 <2n−1>+ an xn <2n>= (1+ θ1 ) a0 + (1+ θ3 ) a1 x+⋯+ (1+ θ2n−1 ) an−1 xn−1 + (1+ θ2n ) an xn , 5.2 where we have used Lemma 3.1, and where |θ_k| ≤ ku/(l − ku) ≕ γ_k. This result shows that Horner's method has a small backward error: the computed q̑₀ is the exact value at x of a polynomial obtained by making relative perturbations of size at most γ_2n to the coefficients of p(x).

While it is true that all norms are equivalent theoretically, only a homely one like the ∞-norm is truly useful numerically

— J. H. WILKINSON, Lecture at Stanford University (1984)

Matrix norms are defined in many different ways in the older literature, but the favorite was the Euclidean norm of the matrix considered as a vector in n²-space. Wedderburn (1934) calls this the absolute value of the matrix and traces the idea back to Peano in 1887.

— ALSTON S. HOUSEHOLDER,

The Theory of Matrices in Numerical Analysis (1964)

Norms are an indispensable tool in numerical linear algebra. Their ability to compress the mn numbers in an m × n matrix into a single scalar measure of size enables perturbation results and rounding error analyses to be expressed in a concise and easily interpreted form. In problems that are badly scaled, or contain a structure such as sparsity, it is often better to measure matrices and vectors componentwise. But norms remain a valuable instrument for the error analyst, and in this chapter we describe some of their most useful and interesting properties.

6.1. Vector Norms

A vector norm is a function ‖ · ‖ : ℂⁿ → ℝ satisfying the following conditions:

1. ‖x‖ ≥ 0 with equality iff x = 0.

2. ‖αx‖ = |α| ‖x‖ for all α ∈ ℂ, x ∈ ℂⁿ.

3. ‖x + y‖ ≤ ‖x‖ + ‖y‖ for all x, y ∈ ℂⁿ (the triangle inequality).

The three most useful norms in error analysis and in numerical computation are ‖x ‖1 = ∑ i=1 n | xi | ,"Manhattan" or "taxi cab" norm,‖x ‖2 = ( ∑ i=1 n | xi |2 )1/2 = (x*x)1/2 ,Euclidean length,‖x ‖∞ = max 1≤i≤n | xi | . These are all special cases of the Hölder p-norm: ‖x ‖p = ( ∑ i=1 n | xi |p )1/p , p≥1.

The 2-norm has two properties that make it particularly useful for theoretical purposes. First, it is invariant under unitary transformations, for if Q*Q = I, then ‖Qx ‖ 2 2 =x*Q*Qx=x*x= ‖x ‖ 2 2 . Second, the 2-norm is differentiable for all x, with gradient vector ∇‖x‖₂ = x/ ‖x‖₂.

A fundamental inequality for vectors is the Hölder inequality (see, for example, [547, 1967, App. 1]) |x*y|≤ ‖x ‖p ‖y ‖q , 1 p + 1 q =1. 6.1 This is an equality when p, q > 1 if the vectors (|x_i|^p) and (|y_i|^q) are linearly dependent and x_iy_i lies on the same ray in the complex plane for all i; equality is also possible when p = 1 and p = ∞, as is easily verified. The special case with p = q = 2 is called the Cauchy-Schwarz inequality: |x*y|≤ ‖x ‖2 ‖y ‖2 .

Our hero is the intrepid, yet sensitive matrix A. Our villain is E, who keeps perturbing A. When A is perturbed he puts on a crumpled hat: Ã = A + E.

— G. W. STEWART and JI-GUANG SUN, Matrix Perturbation Theory (1990)

The expression ‘ill-conditioned’ is sometimes used merely as a term of abuse applicable to matrices or equations … It is characteristic of ill-conditioned sets of equations that small percentage errors in the coefficients given may lead to large percentage errors in the solution.

— A. M. TURING, Rounding-Off Errors in Matrix Processes (1948)

In this chapter we are concerned with a linear system Ax = b, where A ∈ ℝ^{n × n}. In the context of uncertain data or inexact arithmetic there are three important questions:

(1) How much does x change if we perturb A and b; that is, how sensitive is the solution to perturbations in the data?

(2) How much do we have to perturb the data A and b for an approximate solution y to be the exact solution of the perturbed system—in other words, what is the backward error of y?

(3) What bound should we compute in practice for the forward error of a given approximate solution?

To answer these questions we need both normwise and componentwise perturbation theory.

7.1. Normwise Analysis

First, we present some classical normwise perturbation results. We denote by ‖ · ‖ any vector norm and the corresponding subordinate matrix norm. As usual, κ(A) = ‖A‖ ‖A⁻¹‖ is the matrix condition number. Throughout this chapter the matrix E and the vector ƒ are arbitrary and represent tolerances against which the perturbations are measured (their role becomes clear when we consider componentwise results).

Our first result makes precise the intuitive feeling that if the residual is small then we have a “good” approximate solution.

Theorem 7.1 (Rigal and Gaches). The normwise backward error ηE,ƒ (y)≔min {ϵ : (A+ΔA)y=b+Δb, ‖ΔA‖≤ϵ ‖E‖, ‖Δb‖ ≤ϵ ‖ƒ‖ } 7.1 is given by ηE,ƒ (y)= ‖r‖ ‖E‖ ‖y‖ + ‖ƒ‖ , 7.2 where r = b − Ay.

Proof. It is straightforward to show that the right-hand side of (7.2) is a lower bound for η_{E, ƒ} (y). This lower bound is attained for the perturbations Δ Amin = ‖E‖ ‖y‖ ‖E‖ ‖y‖ + ‖ƒ‖ rzT ,Δ bmin = ‖ƒ‖ ‖E‖ ‖y‖ + ‖ƒ‖ r, 7.3 where z is a vector dual to y (see §6.1).

For the particular choice E = A and ƒ = b, η_{E, ƒ}(y) is called the normwise relative backward error.

The next result measures the sensitivity of the system.

Theorem 7.2. Let Ax = b and (A + ΔA)y = b + Δb, where ‖ΔA‖ ≤ ϵ‖E‖ and ‖Δb‖ ≤ ϵ‖ ƒ ‖, and assume that ϵ‖ A⁻¹‖ ‖E‖ < 1.

In the end there is left the coefficient of one unknown and the constant term. An elimination between this equation and one from the previous set that contains two unknowns yields an equation with the coefficient of another unknown and another constant term, etc. The quotient of the constant term by the unknown yields the value of the unknown in each case.

— JOHN V. ATANASOFF, Computing Machine for the Solution of Large Systems of Linear Algebraic Equations (1940)

The solutions of triangular systems are usually computed to high accuracy. This fact … cannot be proved in general, for counter examples exist. However, it is true of many special kinds of triangular matrices and the phenomenon has been observed in many others. The practical consequences of this fact cannot be over-emphasized.

— G. W. STEWART, Introduction to Matrix Computations (1973)

In practice one almost invariably finds that if L is ill-conditioned, so that ‖L‖ ‖L⁻¹‖ ≫ 1, then the computed solution of Lx = b (or the computed inverse) is far more accurate than [standard norm bounds] would suggest.

— J. H. WILKINSON, Rounding Errors in Algebraic Processes (1963)

Triangular systems play a fundamental role in matrix computations. Many methods are built on the idea of reducing a problem to the solution of one or more triangular systems, including virtually all direct methods for solving linear systems. On serial computers triangular systems are universally solved by the standard back and forward substitution algorithms. For parallel computation there are several alternative methods, one of which we analyse in §8.4.

Backward error analysis for the substitution algorithms is straightforward and the conclusion is well known: the algorithms are extremely stable. The behaviour of the forward error, however, is intriguing, because the forward error is often surprisingly small—much smaller than we would predict from the normwise condition number κ, or, sometimes, even the componentwise condition number cond. The quotes from Stewart and Wilkinson at the start of this chapter emphasize the high accuracy that is frequently observed in practice. The analysis we give in this chapter provides a partial explanation for the observed accuracy of the substitution algorithms. In particular, it reveals three important but nonobvious properties:

• the accuracy of the computed solution from substitution depends strongly on the right-hand side;

• a triangular matrix may be much more or less ill conditioned than its transpose; and

• the use of pivoting in LU, QR, and Cholesky factorizations can greatly improve the conditioning of a resulting triangular system.

As well as deriving backward and forward error bounds, we show how to compute upper and lower bounds for the inverse of a triangular matrix.

8.1. Backward Error Analysis

Recall that for an upper triangular matrix U ∈ ℝ^{n × n} the system Ux = b can be solved using the formula xi = ( bi − ∑ j=i+1 n uij xj )/ uii , which yields the components of x in order from last to first.

Algorithm 8.1 (back substitution). Given a nonsingular upper triangular matrix U ∈ ℝ^{n × n} this algorithm solves the system Ux = b. xn = bn / unn for i=n−1:−1:1s= bi for j=i+1:ns=s− uij xj end xi =s/ uii end

Cost: n² flops.

It appears that Gauss and Doolittle applied the method only to symmetric equations. More recent authors, for example, Aitken, Banachiewicz, Dwyer, and Crout … have emphasized the use of the method, or variations of it, in connection with non-symmetric problems … Banachiewicz … saw the point … that the basic problem is really one of matrix factorization, or “decomposition” as he called it.

— PAUL S. DWYER, Linear Computations (1951)

Intolerable pivot-growth [with partial pivoting] is a phenomenon that happens only to numerical analysts who are looking for that phenomenon.

— WILLIAM M. KAHAN, Numerical Linear Algebra (1966)

By 1949 the major components of the Pilot ACE were complete and undergoing trials … During 1951 a programme for solving simultaneous linear algebraic equations was used for the first time. 26th June, 1951 was a landmark in the history of the machine, for on that day it first rivalled alternative computing methods by yielding by 3 p.m. the solution to a set of 17 equations submitted the same morning.

— MICHAEL WOODGER, The History and Present Use of Digital Computers at the National Physical Laboratory (1958).

The closer one looks, the more subtle and remarkable Gaussian elimination appears.

— LLOYD N. TREFETHEN, Three Mysteries of Gaussian Elimination (1985)

9.1. Gaussian Elimination and Pivoting Strategies

We begin by giving a traditional description of Gaussian elimination (GE) for solving a linear system Ax = b, where A ∈ ℝ^{n × n} is nonsingular.

The strategy of GE is to reduce a problem we can't solve (a full linear system) to one that we can (a triangular system), using elementary row operations. There are n − 1 stages, beginning with A⁽¹⁾ ≔ A, b⁽¹⁾ ≔ b, and finishing with the upper triangular system A⁽ⁿ⁾x = b⁽ⁿ⁾

At the kth stage we have converted the original system to A^(k)x = b^(k), where A (k) = [ A 11 (k) A 12 (k) 0 A 22 (k) ], with A 11 (k) ∈ ℝ (k−1)× (k−1) upper triangular. The purpose of the kth stage of the elimination is to zero the elements below the diagonal in the kth column of A^(k). This is accomplished by the operations a ij (k+1) = a ij (k) − mik a kj (k) ,i=k+1:n,j=k+1:n, b i (k+1) = b i (k) − mik b k (k) ,i=k+1:n, where the multipliers mik = a ik (k) / a kk (k) , i = k + 1: n. At the end of the (n — 1)st stage we have the upper triangular system A⁽ⁿ⁾x = b⁽ⁿ⁾, which is solved by back substitution. For an n × n matrix, GE requires 2n³/3 flops.

There are two problems with the method as described. First, there is a breakdown with division by zero if a kk (k) =0 . Second, if we are working in finite precision and some multiplier m_ik is large, then there is a possible loss of significance: in the subtraction a ik (k) − mik a kj (k) , low-order digits of a ik (k) could be lost. Losing these digits could correspond to making a relatively large change to the original matrix A. The simplest example of this phenomenon is for the matrix ϵ111 ; here, a 22 (2) =1−1/ϵ and fl ( a 22 (2) ) =−1/ϵ if ϵ<u if ϵ < u, which would be the exact answer if we changed a₂₂ from 1 to 0.

These observations motivate the strategy of partial pivoting. At the start of the kth stage, the kth and rth rows are interchanged, where | a rk (k) | ≔ max k≤i≤n | a ik (k) | . Partial pivoting ensures that the multipliers are nicely bounded: | mik | ≤1,i=k+1:n.

A more expensive pivoting strategy, which interchanges both rows and columns, is complete pivoting.

At the start of the kth stage, rows k and r and columns k and s are interchanged, where | a rs (k) | ≔ max k≤i,j≤n | a ij (k) | . Note that this requires O(n³) comparisons in total, compared with O(n²) for partial pivoting.

The matrix of that equation system is negative definite—which is a positive definite system that has been multiplied through by −1. For all practical geometries the common finite difference Laplacian operator gives rise to these, the best of all possible matrices. Just about any standard solution method will succeed, and many theorems are available for your pleasure.

— FORMAN S. ACTON, Numerical Methods That Work (1970)

Many years ago we made out of half a dozen transformers a simple and rather inaccurate machine for solving simultaneous equations—the solutions being represented as flux in the cores of the transformers. During the course of our experiments we set the machine to solve the equations— X+Y+Z=1X+Y+Z=2X+Y+Z=3 The machine reacted sharply—it blew the main fuse and put all the lights out.

— B. V. BOWDEN, The Organization of a Typical Machine (1953)

There does seem to be some misunderstanding about the purpose of an a priori backward error analysis. All too often, too much attention is paid to the precise error bound that has been established. The main purpose of such an analysis is either to establish the essential numerical stability of an algorithm or to show why it is unstable and in doing so to expose what sort of change is necessary to make it stable. The precise error bound is not of great importance.

— J. H. WILKINSON, Numerical Linear Algebra on Digital Computers (1974)

10.1. Symmetric Positive Definite Matrices

Symmetric positive definiteness is one of the highest accolades to which a matrix can aspire. Symmetry confers major advantages and simplifications in the eigen-problem and, as we will see in this chapter, positive definiteness permits economy and numerical stability in the solution of linear systems.

A symmetric matrix A ∈ ℝ^{n × n} is positive definite if x^TAx > 0 for all nonzero x ∈ ℝⁿ. Well-known equivalent conditions to A = A^T being positive definite are

• det(A_k) > 0, k = 1: n, where A_k = A(1: k, 1: k) is the leading principal submatrix of order k.

• λ_k(A) > 0, k = 1: n, where λ_k denotes the kth largest eigenvalue.

The first of these conditions implies that A has an LU factorization, A = LU (see Theorem 9.1). Another characterization of positive definiteness is that the pivots in LU factorization are positive, since u_kk = det(A_k)/ det(A_{k − 1}). By factoring out the diagonal of U and taking its square root, the LU factorization can be converted into a Cholesky factorization: A = R^TR, where R is upper triangular with positive diagonal elements. This factorization is so important that it merits a direct proof.

Theorem 10.1. If A ∈ ℝ^{n × n} is symmetric positive definite then there is a unique upper triangular R ∈ ℝ^{n × n} with positive diagonal elements such that A = R^TR.

Proof. The proof is by induction. The result is clearly true for n = 1. Assume it is true for n − 1. The leading principal submatrix A_{n − 1} = A(1: n − 1, 1: n − 1) is positive definite, so it has a unique Cholesky factorization An−1 = R n−1 T Rn−1 . We have a factorization A= [ An−1 c cT α]= [ R n−1 T 0 rT β] [ Rn−1 r0β]≔ RT R 10.1 if R n−1 T r=c, 10.2 rT r+ β2 =α. 10.3 Equation (10.2) has a unique solution since R_{n − 1} is nonsingular. Then (10.3) gives β² = α − r^Tr. It remains to check that there is a unique real, positive β satisfying this equation. From the equation 0<det (A) =det ( RT ) det (R)=det ( Rn−1 ) 2 β2 we see that β² > 0, hence there is a unique β > 0.

The proof of the theorem is constructive, and provides a way to compute the Cholesky factorization that builds R a column at a time.

Unfortunately, there is no stable scheme exactly analogous to Gaussian elimination with partial pivoting; one cannot construct an algorithm for which there is a bound on the element growth of the sequence A^(k) when at each stage only one column of A^(k) is examined.

— JAMES R. BUNCH and LINDA KAUFMAN, Some Stable Methods for Calculating Inertia and Solving Symmetric Linear Systems (1977)

The comparison count is much less … than for the complete pivoting scheme, and in practice this fact has had a much larger impact than originally anticipated, sometimes cutting the execution time by about 40%.

— JAMES R. BUNCH and LINDA KAUFMAN, Some Stable Methods for Calculating Inertia and Solving Symmetric Linear Systems (1977)

This fixed inertia property is why skew-symmetric matrices are easier to decompose than symmetric indefinite matrices.

— JAMES R. BUNCH, A Note on the Stable Decomposition of Skew-Symmetric Matrices (1982)

A symmetric matrix A ∈ ℝ^{n × n} is indefinite if (x^TAx)(y^TAy) < 0 for some x, y ∈ ℝⁿ, or, equivalently, if A has both positive and negative eigenvalues. Linear systems with symmetric indefinite coefficient matrices arise in many applications, including least squares problems, optimization, and discretized incompressible Navier-Stokes equations. How can we solve Ax = b efficiently?

Gaussian elimination with partial pivoting (GEPP) can be used to compute the factorization PA = LU, but it does not take advantage of the symmetry to reduce the cost and storage. We might try to construct a factorization A = LDL^T, where L is unit lower triangular and D is diagonal. But this factorization may not exist, even if symmetric pivoting is allowed, and if it does exist its computation may be unstable. For example, consider A≡ PAPT = [ϵ11ϵ]=[101/ϵ1] [ϵ00ϵ−1/ϵ] [11/ϵ01]. There is arbitrarily large element growth for 0 < ϵ ≪ 1, and the factorization does not exist for ϵ = 0.

For solving dense symmetric indefinite linear systems Ax = b, two types of factorization are used. The first is the block LDL^T factorization (or symmetric indefinite factorization) PAPT = LDLT , 11.1 where P is a permutation matrix, L is unit lower triangular, and D is block diagonal with diagonal blocks of dimension 1 or 2. This factorization is essentially a symmetric block form of GE, with symmetric pivoting. Note that by Sylvester's inertia theorem, A and D have the same inertia, which is easily determined from D (see Problem 11.2).

The second factorization is that produced by Aasen's method, PAPT = LTLT , where P is a permutation matrix, L is unit lower triangular with first column e₁, and T is tridiagonal. Aasen's factorization is much less widely used than block LDL^T factorization, but it is mathematically interesting. We describe both factorizations in this chapter.

A matrix A ∈ ℝ^{n × n} is skew-symmetric if A^T = −A. In the final section of this chapter we describe a form of block LDL^T factorization specialized to skew-symmetric matrices.

The ILLIAC's memory is sufficient to accommodate a system of 39 equations when used with Routine 51. The additional length of Routine 100 restricts to 37 the number of equations that it can handle. With 37 equations the operation time of Routine 100 is about 4 minutes per iteration.

— JAMES N. SNYDER, On the Improvement of the Solutions to a Set of Simultaneous Linear Equations Using the ILLIAC (1955)

In a short mantissa computing environment the presence of an iterative improvement routine can significantly widen the class of solvable Ax = b problems.

— GENE H. GOLUB and CHARLES F. VAN LOAN, Matrix Computations (1996)

Most problems involve inexact input data and obtaining a highly accurate solution to an imprecise problem may not be justified.

— J. J. DONGARRA, J. R. BUNCH, C. B. MOLER, and G. W. STEWART, LINPACK Users' Guide (1979)

It is shown that even a single iteration of iterative refinement in single precision is enough to make Gaussian elimination stable in a very strong sense.

— ROBERT D. SKEEL, Iterative Refinement Implies Numerical Stability for Gaussian Elimination (1980)

Iterative refinement is an established technique for improving a computed solution x̑ to a linear system Ax = b. The process consists of three steps:

1. Compute r = b − Ax̑.

2. Solve Ad = r.

3. Update y = x̑ + d.

(Repeat from step 1 if necessary, with x̑ replaced by y.)

If there were no rounding errors in the computation of r, d, and y, then y would be the exact solution to the system. The idea behind iterative refinement is that if r and d are computed accurately enough then some improvement in the accuracy of the solution will be obtained. The economics of iterative refinement are favourable for solvers based on a factorization of A, because the factorization used to compute x̑ can be reused in the second step of the refinement.

Traditionally, iterative refinement is used with Gaussian elimination (GE), and r is computed in extended precision before being rounded to working precision. Iterative refinement for GE was used in the 1940s on desk calculators, but the first thorough analysis of the method was given by Wilkinson in 1963 [1232, 1963]. The behaviour of iterative refinement for GE is usually summarized as follows: if double the working precision is used in the computation of r, and A is not too ill conditioned, then the iteration produces a solution correct to working precision and the rate of convergence depends on the condition number of A. In the next section we give a componentwise analysis of iterative refinement that confirms this summary and provides some further insight.

Iterative refinement is of interest not just for improving the accuracy or stability of an already stable linear system solver, but also for ameliorating instability in a less reliable solver. Examples of this usage are in Li and Demmel [786, 1998] and Dongarra, Eijkhout, and Luszczek [345, 2000], where sparse GE is performed without pivoting, for speed, and iterative refinement is used to regain stability.

12.1. Behaviour of the Forward Error

Let A ∈ ℝ^{n × n} be nonsingular and let x̑ be a computed solution to Ax = b. Define x₁ = x̑ and consider the following iterative refinement process: r_i = b − Ax_i (precision u‾), solve Ad_i = r_i (precision u), x_{i + 1} = x_i + d_i (precision u), i = 1,2, …. For traditional iterative refinement, u‾ = u². Note that in this chapter subscripts specify members of a vector sequence, not vector elements.

We henceforth define r_i, d_i, and x_i to be the computed quantities (to avoid a profusion of hats). The only assumption we will make on the solver is that the computed solution y̑ to a system Ay = c satisfies (A+ΔA) x ̑ =b, |ΔA| ≤ uW, 12.1 where W is a nonnegative matrix depending on A, n, and u (but not on b). Thus the solver need not be LU factorization or even a factorization method.

The page or so of analysis that follows is straightforward but tedious. The reader is invited to jump straight to (12.5), at least on first reading.

Block algorithms are advantageous for at least two important reasons. First, they work with blocks of data having b² elements, performing O(b³) operations. The O(b) ratio of work to storage means that processing elements with an O(b) ratio of computing speed to input/output bandwidth can be tolerated. Second, these algorithms are usually rich in matrix multiplication. This is an advantage because nearly every modern parallel machine is good at matrix multiplication.

— ROBERT S. SCHREIBER, Block Algorithms for Parallel Machines (1988)

It should be realized that, with partial pivoting, any matrix has a triangular factorization. DECOMP actually works faster when zero pivots occur because they mean that the corresponding column is already in triangular form.

— GEORGE E. FORSYTHE, MICHAEL A. MALCOLM, and CLEVE B. MOLER, Computer Methods for Mathematical Computations (1977)

It was quite usual when dealing with very large matrices to perform an iterative process as follows: the original matrix would be read from cards and the reduced matrix punched without more than a single row of the original matrix being kept in store at any one time; then the output hopper of the punch would be transferred to the card reader and the iteration repeated.

— MARTIN CAMPBELL-KELLY, Programming the Pilot ACE (1981)

13.1. Block Versus Partitioned LU Factorization

As we noted in Chapter 9 (Notes and References), Gaussian elimination (GE) comprises three nested loops that can be ordered in six ways, each yielding a different algorithmic variant of the method. These variants involve different computational kernels: inner product and saxpy operations (level-1 BLAS), or outer product and gaxpy operations (level-2 BLAS). To introduce matrix-matrix operations (level-3 BLAS), which are beneficial for high-performance computing, further manipulation beyond loop reordering is needed. We will use the following terminology, which emphasises an important distinction.

A partitioned algorithm is a scalar (or point) algorithm in which the operations have been grouped and reordered into matrix operations.

A block algorithm is a generalization of a scalar algorithm in which the basic scalar operations become matrix operations (α + β, αβ, and α/β become A + B, AB, and AB⁻¹), and a matrix property based on the nonzero structure becomes the corresponding property blockwise (in particular, the scalars 0 and 1 become the zero matrix and the identity matrix, respectively). A block factorization is defined in a similar way and is usually what a block algorithm computes.

A partitioned version of the outer product form of LU factorization may be developed as follows. For A ∈ ℝ^{n × n} and a given block size r, write [ A11 A12 A21 A22 ] =[ L11 0 L21 In−r ] [ Ir 00S] [ U11 U12 0 In−r ], 13.1 where A₁₁ is r × r.

Algorithm 13.1 (partitioned LU factorization). This algorithm computes an LU factorization A = LU ∈ ℝ^{n × n} using a partitioned outer product implementation, using block size r and the notation (13.1).

1. Factor A₁₁ = L₁₁U₁₁.

2. Solve L₁₁U₁₂ = A₁₂ for U₁₂.

3. Solve L₂₁U₁₁ = A₂₁ for L₂₁.

4. Form S = A₂₂ − L₂₁U₁₂.

5. Repeat steps 1–4 on S to obtain L₂₂ and U₂₂.

Note that in step 4, S= A22 − A21 A 11 −1 A12 is the Schur complement of A₁₁ in A. Steps 2 and 3 require the solution of the multiple right-hand side triangular systems, so steps 2–4 are all level-3 BLAS operations. This partitioned algorithm does precisely the same arithmetic operations as any other variant of GE, but it does the operations in an order that permits them to be expressed as matrix operations.

A genuine block algorithm computes a block L U factorization, which is a factorization A = LU ∈ ℝ^{n × n}, where L and U are block triangular and L has identity matrices on the diagonal: L= I L21 I⋮⋱ Lm1 ⋯ Lm,m−1 I ,U= U11 U12 ⋯ U1m U22 ⋮⋱ Um−1,m Umm .

It is amusing to remark that we were so involved with matrix inversion that we probably talked of nothing else for months. Just in this period Mrs. von Neumann acquired a big, rather wild but gentle Irish Setter puppy, which she called Inverse in honor of our work!

— HERMAN H. GOLDSTINE, The Computer: From Pascal to von Neumann (1972)

The most computationally intensive portion of the tasks assigned to the processors is integrating the KKR matrix inverse over the first Brillouin zone. To evaluate the integral, hundreds or possibly thousands of complex double precision matrices of order between 80 and 300 must be formed and inverted. Each matrix corresponds to a different vertex of the tetrahedrons into which the Brillouin zone has been subdivided.

— M. T. HEATH, G. A. GEIST, and J. B. DRAKE, Superconductivity Computations, in Early Experience with the Intel IPSC/860 at Oak Ridge National Laboratory (1990)

Press 1/x to invert the matrix. Note that matrix inversion can produce erroneous results if you are using ill-conditioned matrices.

— HEWLETT-PACKARD, HP 48G Series User's Guide (1993)

Almost anything you can do with A⁻¹ can be done without it.

— GEORGE E. FORSYTHE and CLEVE B. MOLER, Computer Solution of Linear Algebraic Systems (1967)

14.1. Use and Abuse of the Matrix Inverse

To most numerical analysts, matrix inversion is a sin. Forsythe, Malcolm, and Moler put it well when they say [430, 1977, p. 31] “In the vast majority of practical computational problems, it is unnecessary and inadvisable to actually compute A⁻¹.” The best example of a problem in which the matrix inverse should not be computed is the linear equations problem Ax = b. Computing the solution as x = A⁻¹ × b requires 2n³ flops, assuming A⁻¹ is computed by Gaussian elimination with partial pivoting (GEPP), whereas GEPP applied directly to the system costs only 2n³/3 flops.

Not only is the inversion approach three times more expensive, but it is much less stable. Suppose X = A⁻¹ is formed exactly, and that the only rounding errors are in forming x = fl(Xb). Then x ^ = (X + ΔX)b, where |ΔX| ≤ γ_n|X|, by (3.11). So A x ^ = A(X + ΔX)b = (I + AΔX)b, and the best possible residual bound is |b−A x ̑ | ≤ γn |A| | A−1 | |b| . For GEPP, Theorem 9.4 yields |b−A x ̑ | ≤ γ3n | L ̑ | | U ̑ | | x ̑ | . Since it is usually true that ‖ | L ^ ‖Û| ‖_∞ ≈ ‖A‖_∞ for GEPP, we see that the matrix inversion approach is likely to give a much larger residual than GEPP if A is ill conditioned and if ‖x‖_∞ ≪ ‖ | A⁻¹ | | b| ‖_∞. For example, we solved fifty 25 × 25 systems Ax = b in MATLAB, where the elements of x are taken from the normal N(0, 1) distribution and A is random with κ₂(A) = u^−1/2 ≈ 9 × 10⁷. As Table 14.1 shows, the inversion approach provided much larger backward errors than GEPP in this experiment.

Given the inexpedience of matrix inversion, why devote a chapter to it? The answer is twofold. First, there are situations in which a matrix inverse must be computed. Examples are in statistics [63, 1974, §7.5], [806, 1984, §2.3], [834, 1989, p. 342 ff], where the inverse can convey important statistical information, in certain matrix iterations arising in eigenvalue-related problems [47, 1997], [193, 1987], [621, 1994], in the computation of matrix functions such as the square root [609, 1997] and the logarithm [231, 2001], and in numerical integrations arising in superconductivity computations [555, 1990] (see the quotation at the start of the chapter). Second, methods for matrix inversion display a wide variety of stability properties, making for instructive and challenging error analysis. (Indeed, the first major rounding error analysis to be published, that of von Neumann and Goldstine, was for matrix inversion—see §9.13).

Most of LAPACK's condition numbers and error bounds are based on estimated condition numbers … The price one pays for using an estimated rather than an exact condition number is occasional (but very rare) underestimates of the true error; years of experience attest to the reliability of our estimators, although examples where they badly underestimate the error can be constructed.

— E. ANDERSON et al., LAPACK Users' Guide (1999)

The importance of the counter-examples is that they make clear that any effort toward proving that the algorithms always produce useful estimations is fruitless. It may be possible to prove that the algorithms produce useful estimations in certain situations, however, and this should be pursued. An effort simply to construct more complex algorithms is dangerous.

— A. K. CLINE and R. K. REW, A Set of Counter-Examples to Three Condition Number Estimators (1983)

Singularity is almost invariably a clue.

— SIR ARTHUR CONAN DOYLE, The Boscombe Valley Mystery (1892)

Condition number estimation is the process of computing an inexpensive estimate of a condition number, where “inexpensive” usually means that the cost is an order of magnitude less than would be required to compute the condition number exactly. Condition estimates are required when evaluating forward error bounds for computed solutions to many types of linear algebra problems. An estimate of the condition number that is correct to within a factor 10 is usually acceptable, because it is the magnitude of an error bound that is of interest, not its precise value. In this chapter we are concerned with condition numbers involving A⁻¹, where A is a given nonsingular n × n matrix; these arise in linear systems Ax = b but also more generally in, for example, least squares problems and eigenvalue problems.

An important question for any condition estimator is whether there exists a “counterexample” : a parametrized matrix for which the quotient “estimate divided by true condition number” can be made arbitrarily small (or large, depending on whether the estimate is a lower bound or an upper bound) by varying a parameter. Intuitively, it might be expected that any condition estimator has counterexamples or, equivalently, that estimating a condition number to within a constant factor independent of the matrix A must cost as much, asymptotically, as computing A⁻¹. This property was conjectured by Demmel [314, 1992] and strong evidence for it is provided by Demmel, Diament, and Malajovich [321, 2001]. In the latter paper it is shown that the cost of computing an estimate of ‖A⁻¹ ‖ of guaranteed quality is at least the cost of testing whether the product of two n × n matrices is zero, and performing this test is conjectured to cost as much as actually computing the product [185, 1997, Problem 16.3].

15.1. How to Estimate Componentwise Condition Numbers

In the perturbation theory of Chapter 7 for linear equations we obtained perturbation bounds that involve the condition numbers κE,ƒ (A,x)= ‖ A−1 ‖ ‖ƒ‖ ‖x‖ + ‖ A−1 ‖ ‖E‖, condE,ƒ (A,x)= ‖ | A−1 | ( |E|x | + ƒ) ‖∞ ‖x ‖∞ , and their variants. To compute these condition numbers exactly we need to compute A⁻¹, which requires O(n³) operations, even if we are given a factorization of A. Is it possible to produce reliable estimates of both condition numbers in O(n²) operations? The answer is yes, but to see why we first need to rewrite the expression for cond_{E, f}(A, x). Consider the general expression ‖ |Ais usually true that ‖ | L ^ |‾‖Û| ‖_∞⁻¹|d‖_∞, where d is a given nonnegative vector (thus, d = f + E|x| for cond_{E, f}(A, x)); note that the practical error bound (7.31) is of this form. Writing D = diag(d) and e = [1,1, …, 1]^T , we have ‖| A−1 | d ‖∞ =‖| A−1 | De ‖∞ =‖| A−1 D|e ‖∞ =‖| A−1 D| ‖∞ =‖ A−1 D ‖∞ . 15.1

We must commence, not with a square, but with an oblong arrangement of terms consisting, suppose, of m lines and n columns. This will not in itself represent a determinant, but is, as it were, a Matrix out of which we may form various systems of determinants by fixing upon a number p, and selecting at will p lines and p columns, the squares corresponding to which may be termed determinants of the pth order.

— J. J. SYLVESTER, Additions to the Articles, “On a New Class of Theorems,” and “On Pascal's Theorem” (1850)

I have in previous papers defined a “Matrix” as a rectangular array of terms, out of which different systems of determinants may be engendered, as from the womb of a common parent; these cognate determinants being by no means isolated in their relations to one another, but subject to certain simple laws of mutual dependence and simultaneous deperition.

— J. J. SYLVESTER, On the Relation Between the Minor Determinants of Linearly Equivalent Quadratic Functions (1851)

The linear matrix equation AX−XB=C, 16.1 where A ∈ ℝ^{m × m}, B ∈ ℝ^{n × n}, and C ∈ ℝ^{m × n} are given and X ∈ ℝ^{m × n} is to be determined, is called the Sylvester equation. It is of pedagogical interest because it includes as special cases several important linear equation problems:

1. linear system: Ax = c,

2. multiple right-hand side linear system: AX = C,

3. matrix inversion: AX = I,

4. eigenvector corresponding to given eigenvalue b: (A − bI)x = 0,

5. commuting matrices: AX − XA = 0.

The Sylvester equation arises in its full generality in various applications. For example, the equations [I−X0I] [A−C0B] [IX0I] =[AAX−XB−C0B] show that block-diagonalizing a block triangular matrix is equivalent to solving a Sylvester equation. The Sylvester equation can also be produced from finite difference discretization of a separable elliptic boundary value problem on a rectangular domain, where A and B represent application of a difference operator in the “y” and “x” directions, respectively [1059, 1991].

That (16.1) is merely a linear system is emphasized by writing it in the form ( In ⊗A− BT ⊗ Im )vec (X) =vec (C), 16.2 where A ⊗ B := (a_ijB) is a Kronecker product and the vec operator stacks the columns of a matrix into one long vector. For future reference, we note the useful relation vec (AXB) = ( BT ⊕ A) vec (X). (See Horn and Johnson [637, 1991, Chap. 4] for a detailed presentation of properties of the Kronecker product and the vec operator.) The mn × mn coefficient matrix in (16.2) has a very special structure, illustrated for n = 3 by [A− b11 I− b21 I− b31 I− b12 IA− b22 I− b32 I− b13 I− b23 IA− b33 I] . In dealing with the Sylvester equation it is vital to consider this structure and not treat (16.2) as a general linear system.

Since the mn eigenvalues of I_n ⊗ A − B^T ⊗ I_m are given by λij ( In ⊗A− BT ⊗ Im )= λi (A) − λj (B), i=1:m,j=1:n, 16.3 the Sylvester equation is nonsingular precisely when A and B have no eigenvalues in common.

I recommend this method to you for imitation. You will hardly ever again eliminate directly at least not when you have more than 2 unknowns. The indirect [iterative] procedure can be done while half asleep, or while thinking about other things.

— CARL FRIEDRICH GAUSS, Letter to C. L. Gerling (1823)

The iterative method is commonly called the “Seidel process,” or the “Gauss-Seidel process.” But, as Ostrowski (1952) points out, Seidel (1874) mentions the process but advocates not using it. Gauss nowhere mentions it.

— GEORGE E. FORSYTHE, Solving Linear Algebraic Equations Can Be Interesting (1953)

The spurious contributions in null(A) grow at worst linearly and if the rounding errors are small the scheme can be quite effective.

— HERBERT B. KELLER, On the Solution of Singular and Semidefinite Linear Systems by Iteration (1965)

Iterative methods for solving linear systems have a long history, going back at least to Gauss. Table 17.1 shows the dates of publication of selected methods. It is perhaps surprising, then, that rounding error analysis for iterative methods is not well developed. There are two main reasons for the paucity of error analysis. One is that in many applications accuracy requirements are modest and are satisfied without difficulty, resulting in little demand for error analysis. Certainly there is no point in computing an answer to greater accuracy than that determined by the data, and in scientific and engineering applications the data often has only a few correct digits. The second reason is that rounding error analysis for iterative methods is inherently more difficult than for direct methods, and the bounds that are obtained are harder to interpret.

In this chapter we consider a simple but important class of iterative methods, stationary iterative methods, for which a reasonably comprehensive error analysis can be given. The basic question that our analysis attempts to answer is, “What is the limiting accuracy of a method in floating point arithmetic?” Specifically, “How small can we guarantee that the backward or forward error will be over all iterations k = 1,2, …?” Without an answer to this question we cannot be sure that a convergence test of the form ‖b − A x^ _k‖ ≤ ϵ (say) will ever be satisfied, for any given value of ϵ < ‖b − Ax₀‖!

As an indication of the potentially devastating effects of rounding errors we present an example constructed and discussed by Hammarling and Wilkinson [541, 1976]. Here, A is the 100 × 100 lower bidiagonal matrix with a_ii = 1.5 and a_{i, i-1} ≡ 1, and b_i ≡ 2.5. The successive overrelaxation (SOR) method is applied in MATLAB with parameter ω = 1.5, starting with the rounded version of the exact solution x, given by x_i = 1 − (−2/3)ⁱ. The forward errors ‖ x^ _k − x‖_∞/‖x‖_∞ and the ∞-norm backward errors η_A,b( x^ _k) are plotted in Figure 17.1. The SOR method converges in exact arithmetic, since the iteration matrix has spectral radius 1/2, but in the presence of rounding errors it diverges. The iterate x^ ₂₃₈ has a largest element of order 10¹³, x^ _{k +2} ≡ x^ _k for k ≥ 238, and for k > 100, x^ _k(60: 100) ≈ (−1)^k x^ ₁₀₀(60: 100). The divergence is not a result of ill conditioning of A, since κ_∞(A) ≈ 5. The reason for the initial rapid growth of the errors in this example is that the iteration matrix is far from normal; this allows the norms of its powers to become very large before they ultimately decay by a factor ≈ 1/2 with each successive power.

Unfortunately, the roundoff errors in the mth power of a matrix, say B^m, are usually small relative to ‖B‖^m rather than ‖B^m‖.

— CLEVE B. MOLER and CHARLES F. VAN LOAN, Nineteen Dubious Ways to Compute the Exponential of a Matrix (1978)

It is the size of the hump that matters: the behavior of ‖p(AΔt)ⁿ‖ = ‖p(AΔt)^t/Δt‖ for small but nonzero t. The eigenvalues and the norm, by contrast, give sharp information only about the limits t → ∞ or t → 0.

— DESMOND J. HIGHAM and LLOYD N. TREFETHEN, Stiffness of ODEs (1993)

Many people will go through life powering matrices and never see anything as dramatic as J₂(0.99)^k.

— G. W. STEWART, Matrix Algorithms. Volume II: Eigensystems (2001)

Powers of matrices occur in many areas of numerical analysis. One approach to proving convergence of multistep methods for solving differential equations is to show that a certain parameter-dependent matrix is uniformly “power bounded” [537, 1991, §V.7], [974, 1992]. Stationary iterative methods for solving linear equations converge precisely when the powers of the iteration matrix converge to zero. And the power method for computing the largest eigenvalue of a matrix computes the action of powers of the matrix on a vector. It is therefore important to understand the behaviour of matrix powers, in both exact and finite precision arithmetic.

It is well known that the powers A^k of A ∈ ℂ^{n × n} tend to zero as k → ∞ if ρ(A) < 1, where ρ is the spectral radius. However, this simple statement does not tell the whole story. Figure 18.1 plots the 2-norms of the first 30 powers of a certain 3 × 3 matrix with ρ(A) = 0.75. The powers do eventually decay, but initially they grow rapidly. (In this and other similar plots, k on the x-axis is plotted against ‖fl(A^k) ‖₂ on the y-axis, and the norm values are joined by straight lines for plotting purposes.) Figure 18.2 plots the 2-norms of the first 250 powers of a 14 × 14 nilpotent matrix C₁₄ discussed by Trefethen and Trummer [1158, 1987] (see §18.2 for details). The plot illustrates the statement of these authors that the matrix is not power bounded in floating point arithmetic, even though its 14th power should be zero.

These examples suggest two important questions.

• For a matrix with p(A) < 1, how does the sequence {‖A^k ‖} behave? In particular, what is the size of the “hump” max_k ‖A^k‖?

• What conditions on A ensure that fl(A^k) → 0 as k → ∞?

We examine these questions in the next two sections.

18.1. Matrix Powers in Exact Arithmetic

In exact arithmetic the limiting behaviour of the powers of A ∈ ℂ^{n × n} is determined by A's eigenvalues. As already noted, if ρ(A) < 1 then A^k → 0 as k → ∞; if ρ(A) > 1, ‖A^k‖ → ∞ as k → ∞. In the remaining case of ρ(A) = 1, ‖A^k‖ → ∞ if A has a defective eigenvalue λ such that |λ| = 1; A^k does not converge if A has a nondefective eigenvalue λ ≠ 1 such that |λ| = 1 (although the norms of the powers may converge); otherwise, the only eigenvalue of modulus 1 is the nondefective eigenvalue 1 and A^k converges to a nonzero matrix. These statements are easily proved using the Jordan canonical form A=XJ X−1 ∈ ℂn×n , 18.1a where X is nonsingular and J=diag ( J1 , J2 ,… , Js ), Ji = λi 1⋱⋱ λi 1 λi ∈ ℂ ni × ni , 18.1b where n₁ + n₂ + ⋯ + n_s = n. We will call a matrix for which A^k → 0 as k → ∞ (or equivalently, ρ(A) < 1) a convergent matrix.

Any orthogonal matrix can be written as the product of reflector matrices. Thus the class of reflections is rich enough for all occasions and yet each member is characterized by a single vector which serves to describe its mirror.

— BERESFORD N. PARLETT, The Symmetric Eigenvalue Problem (1998)

A key observation for understanding the numerical properties of the modified Gram-Schmidt algorithm is that it can be interpreted as Householder QR factorization applied to the matrix A augmented with a square matrix of zero elements on top. These two algorithms are not only mathematically … but also numerically equivalent. This key observation, apparently by Charles Sheffield, was relayed to the author in 1968 by Gene Golub .

— ÅKE BJöRCK, Numerics of Gram-Schmidt Orthogonalization (1994)

The great stability of unitary transformations in numerical analysis springs from the fact that both the ℓ₂-norm and the Frobenius norm are unitarily invariant. This means in practice that even when rounding errors are made, no substantial growth takes place in the norms of the successive transformed matrices.

— J. H. WILKINSON, Error Analysis of Transformations Based on the Use of Matrices of the Form I − 2ww^H (1965)

The QR factorization is a versatile computational tool that finds use in linear equation, least squares and eigenvalue problems. It can be computed in three main ways. The Gram-Schmidt process, which sequentially orthogonalizes the columns of A, is the oldest method and is described in most linear algebra textbooks. Givens transformations are preferred when A has a special sparsity structure, such as band or Hessenberg structure. Householder transformations provide the most generally useful way to compute the QR factorization. We explore the numerical properties of all three methods in this chapter. We also examine the use of iterative refinement on a linear system solved with a QR factorization and consider the inherent sensitivity of the QR factorization.

19.1. Householder Transformations

A Householder matrix (also known as a Householder transformation, or Householder reflector) is a matrix of the form P=I− 2 υT υ υ υT ,0≠υ ∈ ℝn . It enjoys the properties of symmetry and orthogonality, and, consequently, is involutary (P² = I). The application of P to a vector yields Px=x− ( 2 υT x υT υ ) υ. Figure 19.1 illustrates this formula and makes it clear why P is sometimes called a Householder reflector: it reflects x about the hyperplane span(v)^⊥.

Householder matrices are powerful tools for introducing zeros into vectors. Consider the question, “Given x and y can we find a Householder matrix P such that Px = y?” Since P is orthogonal we clearly require that ‖x‖₂ = ‖y‖₂. Now Px=y ⇔x−2 ( υT x υT υ ) υ=y, and this last equation has the form αυ = x − y for some α. But P is independent of the scaling of υ, so we can set α = 1.

With υ = x − y we have υT υ= xT x+ yT y=2 xT y, and, since x^Tx = y^Ty, υT x= xT x− yT x= 1 2 υT υ. Therefore Px=x−υ=y, as required. We conclude that, provided ‖x‖₂ = ‖y‖₂, we can find a Householder matrix P such that Px = y. (Strictly speaking, we have to exclude the case x = y, which would require υ = 0, making P undefined.)

For some time it has been believed that orthogonalizing methods did not suffer this squaring of the condition number … It caused something of a shock, therefore, when in 1966 Golub and Wilkinson … asserted that already the multiplications QA and Qb may produce errors in the solution containing a factor χ²(A).

— A. VAN DER SLUIS, Stability of the Solutions of Linear Least Squares Problems (1975)

Most packaged regression problems do compute a cross-products matrix and solve the normal equations using a matrix inversion subroutine. All the programs … that disagreed (and some of those that agreed) with the unperturbed solution tried to solve the normal equations.

— ALBERT E. BEATON, DONALD B. RUBIN, and JOHN L. BARONE, The Acceptability of Regression Solutions: Another Look at Computational Accuracy (1976)

On January 1, 1801 Giuseppe Piazzi discovered the asteroid Ceres. Ceres was only visible for forty days before it was lost to view behind the sun … Gauss, using three observations, extensive analysis, and the method of least squares, was able to determine the orbit with such accuracy that Ceres was easily found when it reappeared in late 1801.

— DAVID K. KAHANER, CLEVE B. MOLER, and STEPHEN G. NASH, Numerical Methods and Software (1989)

In this chapter we consider the least squares (LS) problem min_x ‖b − Ax‖₂, where A ∈ ℝ^{m × n} (m ≥ n) has full rank. We begin by examining the sensitivity of the LS problem to perturbations. Then we examine the stability of methods for solving the LS problem, covering QR factorization methods, the normal equations and seminormal equations methods, and iterative refinement. Finally, we show how to compute the backward error of an approximate LS solution.

We do not develop the basic theory of the LS problem, which can be found in standard textbooks (see, for example, Golub and Van Loan [509, 1996, §5.3]). However, we recall the fundamental result that any solution of the LS problem satisfies the normal equations A^TAx = A^Tb (see Problem 20.1). Therefore if A has full rank there is a unique LS solution. More generally, whatever the rank of A the vector x_LS = A⁺b is an LS solution, and it is the solution of minimal 2-norm. Here, A⁺ is the pseudo-inverse of A (given by A⁺ = (A^TA)⁻¹A^T when A has full rank); see Problem 20.3. (For more on the pseudo-inverse see Stewart and Sun [1083, 1990, §3.1].)

20.1. Perturbation Theory

Perturbation theory for the LS problem is, not surprisingly, more complicated than for linear systems, and there are several forms in which bounds can be stated. We begin with a normwise perturbation theorem that is a restatement of a result of Wedin [1211, 1973, Thm. 5.1]. For an arbitrary rectangular matrix A we define the condition number κ₂(A) = ‖A‖₂‖A⁺‖₂. If A has r = rank(A) nonzero singular values, σ₁ ≥ ⋯ ≥ σ_r, then κ₂(A) = σ₁/σ_r.

Theorem 20.1 (Wedin). Let A ∈ ℝ^{m × n} (m ≥ n) and A + ΔA both be of full rank, and let ‖b−Ax‖ 2 =min ,r=b−Ax, ‖ (b+Δb)− (A+ΔA)y ‖ 2 =min ,s=b+Δb− (A+ΔA)y, ‖ΔA‖ 2 ≤ϵ ‖A‖2 , ‖Δb‖ 2 ≤ϵ ‖b‖2 . Then, provided that κ₂(A) ϵ < 1, ‖x−y‖2 ‖x‖2 ≤ κ2 (A)ϵ 1− κ2 (A)ϵ (2+ ( κ2 (A)+1 ) ‖r‖2 ‖A‖ 2 ‖x‖2 ) , 20.1 ‖r−s‖ 2 ‖b‖2 ≤ (1+ 2 κ2 (A) )ϵ . 20.2 These bounds are approximately attainable.

Proof. We defer a proof until §20.10, since the techniques used in the proof are not needed elsewhere in this chapter.

The bound (20.1) is usually interpreted as saying that the sensitivity of the LS problem is measured by κ₂(A) when the residual is small or zero and by κ₂(A)² otherwise. This means that the sensitivity of the LS problem depends strongly on b as well as A, unlike for a square linear system.

I'm thinking of two numbers. Their average is 3. What are the numbers?

— CLEVE B. MOLER, The World's Simplest Impossible Problem (1990)

This problem arises in important algorithms used in mathematical programming … In these cases, B is usually very large and sparse and, because of storage difficulties, it is often uneconomical to store and access Q₁ … Sometimes it has been thought that [the seminormal equations method] could be disastrously worse than [the Q method] … It is the purpose of this note to show that such algorithms are numerically quite satisfactory.

— C. C. PAIGE, An Error Analysis of a Method for Solving Matrix Equations (1973)

Having considered well-determined and over determined linear systems, we now turn to the remaining class of linear systems: those that are underdetermined.

21.1. Solution Methods

Consider the underdetermined system Ax = b, where A ∈ ℝ^{m × n} with m ≤ n. The system can be analysed using a QR factorization AT =Q [R0] , 21.1 where Q ∈ ℝ^{n × n} is orthogonal and R ∈ ℝ^{m × m} is upper triangular. (We could, alternatively, use an LQ factorization of A, but we will keep to the standard notation.) We have b=Ax= [ RT 0] QT x= RT y1 , 21.2 where y= y1 y2 = QT x, y1 ∈ ℝm . If A has full rank then y₁ = R^−Tb is uniquely determined and all solutions of Ax = b are given by x=Q y1 y2 , y2 ∈ ℝn−m arbitrary. The unique solution x_LS that minimizes ‖x‖₂ is obtained by setting y₂ = 0. We have xLS =Q [ R−T b0] 21.3 =Q R0 R−1 R−T b=Q R0 ( RT R)−1 b = AT (A AT )−1 b= A+ b, 21.4 where A⁺ = A^T(AA^T)⁻¹ is the pseudo-inverse of A. Hence x_LS can be characterized as x_ls = A^Ty, where y solves the normal equations AA^Ty = b.

Equation (21.3) defines one way to compute x_LS. We will refer to this method as the “Q method”. When A is large and sparse it is desirable to avoid storing and accessing Q, which can be expensive. An alternative method with this property uses the QR factorization (21.1) but computes x_LS as x_LS = A^Ty, where RT Ry=b 21.5 (cf. (21.4)). These latter equations are called the seminormal equations (SNE). As the “semi” denotes, however, this method does not explicitly form AA^T , which would be undesirable from the standpoint of numerical stability. Note that equations (21.5) are different from the equations R^TRx = A^Tb for an overdetermined least squares (LS) problem, where A = Q[R^T 0]^T ∈ ℝ^{m × n} with m ≥ n, which are also called seminormal equations (see §20.6).

We began, 25 years ago, to take up [the conditioning of] the class of Vandermonde matrices. The original motivation came from unpleasant experiences with the computation of Gauss type quadrature rules from the moments of the underlying weight function.

— WALTER GAUTSCHI, How (Un)stable Are Vandermonde Systems? (1990)

Extreme ill-conditioning of the [Vandermonde] linear systems will eventually manifest itself as n increases by yielding an error curve which is not sufficiently levelled on the current reference … or more seriously fails to have the correct number of sign changes.

— M. ALMACANY, C. B. DUNHAM, and J. WILLIAMS, Discrete Chebyshev Approximation by Interpolating Rationals (1984)

A Vandermonde matrix is defined in terms of scalars α₀, α₁, …, α_n ∈ ℂ by V=V ( α0 , α1 ,… , αn ) = [11…1 α0 α0 … αn ⋮⋮⋮ α 0 n α 1 n … α n n ] ∈ ℂ (n+1)× (n+1) . Vandermonde matrices play an important role in various problems, such as in polynomial interpolation. Suppose we wish to find the polynomial p_n(x) = a_nxⁿ + a_{n − 1}x^{n − 1} + ⋯ + a₀ that interpolates to the data ( αi , ƒi ) i=0 n , for distinct points α_i, that is, p_n(α_i) = f_i, i = 0:n. Then the desired coefficient vector a = [a₀, a₁, …, a_n]^T is the solution of the dual Vandermonde system VT a=ƒ (dual) . The primal system Vx=b (primal) represents a moment problem, which arises, for example, when determining the weights for a quadrature rule: given moments b_i find weights x_i such that ∑ j=0 n xj α j i = bi ,i=0:n .

Because a Vandermonde matrix depends on only n + 1 parameters and has a great deal of structure, it is possible to perform standard computations with reduced complexity. The easiest algorithm to derive is for matrix inversion.

22.1. Matrix Inversion

Assume that V is nonsingular and let V−1 =W= ( wij ) i,j=0 n . The ith row of the equation WV = I may be written ∑ j=0 n wij α k j = δik ,k=0:n. These equations specify a fundamental interpolation problem that is solved by the Lagrange basis polynomial: ∑ j=0 n wij xj = Π k=0 n k≠i ( x− αk αi − αk ) ≕ li (x) . 22.1 The inversion problem is now reduced to finding the coefficients of l_i(x). It is clear from (22.1) that V is nonsingular iff the α_i are distinct.

A simple but extremely valuable bit of equipment in matrix multiplication consists of two plain cards, with a re-entrant right angle cut out of one or both of them if symmetric matrices are to be multiplied. In getting the element of the ith row and jth column of the product, the ith row of the first factor and the jth column of the second should be marked by a card beside, above, or below it.

— HAROLD HOTELLING, Some New Methods in Matrix Calculation (1943)

It was found that multiplication of matrices using punched card storage could be a highly efficient process on the Pilot ACE, due to the relative speeds of the Hollerith card reader used for input (one number per 16 ms.) and the automatic multiplier (2 ms.). While a few rows of one matrix were held in the machine the matrix to be multiplied by it was passed through the card reader. The actual computing and selection of numbers from store occupied most of the time between the passage of successive rows of the cards through the reader, so that the overall time was but little longer than it would have been if the machine had been able to accommodate both matrices.

— MICHAEL WOODGER, The History and Present Use of Digital Computers at the National Physical Laboratory (1958)

23.1. Methods

A fast matrix multiplication method forms the product of two n × n matrices in O(n^ω) arithmetic operations, where ω < 3. Such a method is more efficient asymptotically than direct use of the definition (AB)ij = ∑ k=1 n aik bkj , 23.1 which requires O(n³) operations. For over a century after the development of matrix algebra in the 1850s by Cayley, Sylvester, and others, this definition provided the only known method for multiplying matrices. In 1967, however, to the surprise of many, Winograd found a way to exchange half the multiplications for additions in the basic formula [1249, 1968]. The method rests on the identity, for vectors of even dimension n, xT y= ∑ i=1 n/2 ( x2i−1 + y2i ) ( x2i + y2i−1 ) − ∑ i=1 n/2 x2i−1 x2i − ∑ i=1 n/2 y2i−1 y2i . 23.2 When this identity is applied to a matrix product AB, with x a row of A and y a column of B, the second and third summations are found to be common to the other inner products involving that row or column, so they can be computed once and reused. Winograd's paper generated immediate practical interest because on the computers of the 1960s floating point multiplication was typically two or three times slower than floating point addition. (On today's machines these two operations are usually similar in cost.)

Shortly after Winograd's discovery, Strassen astounded the computer science community by finding an O(n^log₂7) operations method for matrix multiplication (log₂7 ≈ 2.807). A variant of this technique can be used to compute A⁻¹ (see Problem 23.8) and thereby to solve Ax = b, both in O(n^log₂7) operations. Hence the title of Strassen's 1969 paper [1093, 1969], which refers to the question of whether Gaussian elimination is asymptotically optimal for solving linear systems.

Strassen's method is based on a circuitous way to form the product of a pair of 2 × 2 matrices in 7 multiplications and 18 additions, instead of the usual 8 multiplications and 4 additions. As a means of multiplying 2 × 2 matrices the formulae have nothing to recommend them, but they are valid more generally for block 2 × 2 matrices. Let A and B be matrices of dimensions m × n and n × p respectively, where all the dimensions are even, and partition each of A, B, and C = AB into four equally sized blocks: [ C11 C12 C21 C22 ] = [ A11 A12 A21 A22 ] [ B11 B12 B21 B22 ] . 23.3

Once the [FFT] method was established it became clear that it had a long and interesting prehistory going back as far as Gauss. But until the advent of computing machines it was a solution looking for a problem.

— T. W. KöRNER, Fourier Analysis (1988)

Life as we know it would be very different without the FFT.

— CHARLES F. VAN LOAN, Computational Frameworks for the Fast Fourier Transform (1992)

24.1. The Fast Fourier Transform

The matrix-vector product y = F_nx, where Fn = (exp (−2πi (r−1) (s−1)/n ) ) r,s=1 n , is the key computation in the numerical evaluation of Fourier transforms. If the product is formed in the obvious way then O(n²) operations are required. The fast Fourier transform (FFT) is a way to compute y in just O(nlogn) operations. This represents a dramatic reduction in complexity.

The FFT is best understood (at least by a numerical analyst!) by interpreting it as the application of a clever factorization of the discrete Fourier transform (DFT) matrix F_n.

Theorem 24.1 (Cooley-Tukey radix 2 factorization). If n = 2^t then the DFT matrix F_n may be factorized as Fn = At … A1 Pn , 24.1 where P_n is a permutation matrix and Ak = I 2t−k ⊗ B 2k , B 2k = Ir Ωr Ir − Ωr , r= 2k−1 , Ωr =diag (1, ωk ,…, ω k r−1 ), ωk =exp (−2πi/ 2k ) .

Proof. See Van Loan [1182, 1992, Thm. 1.3.3].

The theorem shows that we can write y = F_nx as y= At … A1 Pn x, which is formed as a sequence of matrix-vector products. It is the sparsity of the A_k (two nonzeros per row) that yields the O(n log n) operation count.

We will not consider the implementation of the FFT, and therefore we do not need to define the “bit reversing” permutation matrix P_n in (24.1). However, the way in which the weights ω k j are computed does affect the accuracy. We will assume that computed weights ω ^ k j are used that satisfy, for all j and k, w ̑ k j = w k j + ϵkj ,| ϵkj |≤μ. 24.2 Among the many methods for computing the weights are ones for which we can take μ = cu, μ = cu log j, and μ = cuj, where c is a constant that depends on the method; see Van Loan [1182, 1992, §1.4].

We are now ready to prove an error bound.

Not a single student showed up for Newton's second lecture, and throughout almost every lecture for the next seventeen years … Newton talked to an empty room.

— MICHAEL WHITE, Isaac Newton: The Last Sorcerer (1997)

[On Newton's method] If we start with an approximation to a zero which is appreciably more accurate than the limiting accuracy which we have just described, a single iteration will usually spoil this very good approximation and produce one with an error which is typical of the limiting accuracy

— J. H. WILKINSON, Rounding Errors in Algebraic Processes (1963)

25.1. Newton's Method

Newton's method is a key tool in scientific computing for solving nonlinear equations and optimization problems. Our interest in this chapter is in Newton's method for solving algebraic systems of nonlinear equations.

Let F : ℝⁿ → ℝⁿ be continuously differentiable and denote by J its Jacobian matrix (∂F_i/∂x_j). Given a starting vector x₀, Newton's method for finding a solution of F(x) = 0 generates a sequence {x_i} defined by J ( xi ) ( xi+1 − xi ) =−F ( xi ) , i≥0. 25.1 The attraction of the method stems from the fact that, under appropriate conditions, it converges rapidly to a solution from any starting vector x₀ sufficiently close to that solution. In particular, if the Jacobian is nonsingular at the solution then the rate of convergence is quadratic [333, 1983, Thm. 5.2.1].

Computationally, Newton's method is implemented as solve J ( xi ) di =−F ( xi ) , xi+1 = xi + di . In practical computation rounding and other errors are introduced, so the computed iterates x^ _i actually satisfy x ^ i+1 = x ^ i − (J ( x ^ i )+ Ei )−1 (F ( x ^ i )+ ei ) + εi , 25.2 where

• e_i is the error made when computing the residual F( x^ _i),

• E_i is the error incurred in forming J( x^ i) and solving the linear system for d_i, and

• ε_i is the error made when adding the correction d^ _i to x^ _i.

Note that the term E_i in (25.2) enables the model to cover modified Newton methods, in which the Jacobian is held constant for several iterations in order to reduce the cost of the method.

The question of interest in this chapter is how the various errors in (25.2) affect the convergence of Newton's method. In particular, it is important to know the limiting accuracy and limiting residual, that is, how small ‖x_* − x^ _i‖ and ‖F( x^ _i)‖ are guaranteed to become as n increases, where x_* is the solution to which the iteration would converge in the absence of errors. One might begin an analysis by arguing that (sufficiently small) errors made on a particular iteration are damped out by Newton's method, since the method is locally convergent. However, we are concerned in (25.2) with systematic errors, possibly much larger than the working precision, and how the errors from different iterations interact is not obvious.

In the next section we give general error analysis results for Newton's method in floating point arithmetic that cater for a wide variety of situations: they allow for extended precision in computation of the residual and the use of possibly unstable linear system solvers.

Given the pace of technology, I propose we leave math to the machines and go play outside.

— CALVIN, Calvin and Hobbes by Bill Watterson (1992)

To analyse a given numerical algorithm we proceed as follows. A number which measures the effect of roundoff error is assigned to each set of data. “Hill-climbing” procedures are then applied to search for values large enough to signal instability.

— WEBB MILLER, Software for Roundoff Analysis (1975)

The prospects for effective use of interval arithmetic look very good, so efforts should be made to increase its availability and to make it as user-friendly as possible.

— DONALD E. KNUTH, The Art of Computer Programming, Volume 2, Seminumerical Algorithms (1998)

Despite its wide use, until quite recently the Nelder-Mead method and its ilk have been deprecated, scorned, or ignored by almost all of the mainstream optimization community.

— MARGARET H. WRIGHT, Direct Search Methods: Once Scorned, Now Respectable (1996)

Automatic error analysis is any process in which we use the computer to help us analyse the accuracy or stability of a numerical computation. The idea of automatic error analysis goes back to the dawn of scientific computing. For example, running error analysis, described in §3.3, is a form of automatic error analysis; it was used extensively by Wilkinson on the ACE machine. Various forms of automatic error analysis have been developed. In this chapter we describe in detail the use of direct search optimization for investigating questions about the stability and accuracy of algorithms. We also describe interval analysis and survey other forms of automatic error analysis.

26.1. Exploiting Direct Search Optimization

Is Algorithm X numerically stable? How large can the growth factor be for Gaussian elimination (GE) with pivoting strategy P? By how much can condition estimator C underestimate the condition number of a matrix? These types of questions are common, as we have seen in this book. Usually, we attempt to answer such questions by a combination of theoretical analysis and numerical experiments with random and nonrandom data. But a third approach can be a valuable supplement to the first two: phrase the question as an optimization problem and apply a direct search method.

A direct search method for the problem max x ∈ ℝn ƒ (x), ƒ : ℝn →ℝ 26.1 is a numerical method that attempts to locate a maximizing point using function values only and does not use or approximate derivatives of ƒ. Such methods are usually based on heuristics that do not involve assumptions about the function ƒ. Various direct search methods have been developed; for surveys see Powell [948, 1970], [950, 1998], Swann [1112, 1972], [1113, 1974], Lewis, Torczon, and Trossett [783, 2000], and Wright [1260, 1996]. Most of the methods were developed in the 1960s, in the early years of numerical optimization. For problems in which ƒ is smooth, direct search methods have largely been supplanted by more sophisticated optimization methods that use derivatives (such as quasi-Newton methods and conjugate gradient methods), but they continue to find use in applications where ƒ is not differentiable, or even not continuous. These applications range from chemical analysis [994, 1977], where direct search methods have found considerable use, to the determination of drug doses in the treatment of cancer [105, 1991]; in both applications the evaluation of ƒ is affected by experimental errors. Lack of smoothness of ƒ, and the difficulty of obtaining derivatives when they exist, are characteristic of the optimization problems we consider here.

The use of direct search can be illustrated with the example of the growth factor for GE on A ∈ ℝ^{n × n}, ρn (A) = maxi,j,k | a ij (k) | maxi,j | aij | , where the a ij (k) are the intermediate elements generated during the elimination.

The first American Venus probe was lost due to a program fault caused by the inadvertent substitution of a statement of the form DO 3 I = 1.3 for one of the form DO 3 I = 1, 3.

— JIM HORNING, Note on Program Reliability (1979)

Numerical subroutines should deliver results that satisfy simple, useful mathematical laws whenever possible.

— DONALD E. KNUTH, The Art of Computer Programming, Volume 2, Seminumerical Algorithms (1998)

No method of solving a computational problem is really available to a user until it is completely described in an algebraic computing language and made completely reliable. Before that, there are indeterminate aspects in any algorithm.

— GEORGE E. FORSYTHE, Today's Computational Methods of Linear Algebra (1967)

The extended precision calculation of pi has substantial application as a test of the “global integrity” of a supercomputer … Such calculations … are apparently now used routinely to check supercomputers before they leave the factory. A large-scale calculation of pi is entirely unforgiving; it soaks into all parts of the machine and a single bit awry leaves detectable consequences.

— J. M. BORWEIN, P. B. BORWEIN, and D. H. BAILEY, Ramanujan, Modular Equations, and Approximations to Pi or How to Compute One Billion Digits of Pi (1989)

In this chapter we discuss some miscellaneous aspects of floating point arithmetic that have an impact on software development.

27.1. Exploiting IEEE Arithmetic

IEEE standard 754 and 854 arithmetic can be exploited in several ways in numerical software, provided that proper access to the arithmetic's features is available in the programming environment. Unfortunately, although virtually all modern floating point processors conform to the IEEE standards, language standards and compilers generally provide poor support (exceptions include the Standard Apple Numerics Environment (SANE) [260, 1988], Apple's PowerPC numerics environment [261, 1994], and Sun's SPARCstation compilers [1101, 1992], [1100, 1992]). We give four examples to show the benefits of IEEE arithmetic on software.

Suppose we wish to evaluate the dual of the vector p-norm (1 ≤ p ≤ ∞), that is, the q-norm, where p⁻¹ + q⁻¹ = 1. In MATLAB notation we simply write norm(x, 1/(1−1/p) ), and the extreme cases p = 1 and ∞ correctly yield q = 1/(1 − 1/p) = ∞ and 1 in IEEE arithmetic. (Note that the formula q = p/(p − 1) would not work, because inf/inf evaluates to a NaN.)

The following code finds the biggest element in an array a(1: n), and elegantly solves the problem of choosing a “sufficiently negative number” with which to initialize the variable max: max⁡=−inf for j=1:n if aj >max   max = aj  endend Any unknown or missing elements of the array could be assigned NaNs (a(j) = NaN) and the code would (presumably) still work, since a NaN compares as unordered with everything.

Consider the following continued fraction [691, 1981]: r (x)=7− 3 x−2− 1 x−7+ 10 x−2− 2 x−3 , which is plotted over the range [0,4] in Figure 27.1. Another way to write the rational function r(x) is in the form r (x)= (((7x−101)x+540)x−1240)x+958 (((x−14)x+72)x−151)x+112 , in which the polynomials in the numerator and denominator are written in the form in which they would be evaluated by Horner's rule.

Many tricks or treats associated with the Hilbert matrix may seem rather frightening or fascinating.

— MAN-DUEN CHOI, Tricks or Treats with the Hilbert Matrix (1983)

I start by looking at a 2 by 2 matrix. Sometimes I look at a 4 by 4 matrix. That's when things get out of control and too hard. Usually 2 by 2 or 3 by 3 is enough, and I look at them, and I compute with them, and I try to guess the facts. First, think of a question. Second, I look at examples, and then third, guess the facts.

— PAUL R. HALMOS (1991)

When people look down on matrices, remind them of great mathematicians such as Frobenius, Schur, C. L. Siegel, Ostrowski, Motzkin, Kac, etc., who made important contributions to the subject.

— OLGA TAUSSKY, How I Became a Torchbearer for Matrix Theory (1988)

Ever since the first computer programs for matrix computations were written in the 1940s, researchers have been devising matrices suitable for test purposes and investigating the properties of these matrices. In the 1950s and 1960s it was common for a whole paper to be devoted to a particular test matrix: typically its inverse or eigenvalues would be obtained in closed form.

Early collections of test matrices include those of Newman and Todd [891, 1958] and Rutishauser [1006, 1968]; most of Rutishauser's matrices come from continued fractions or moment problems. Two well-known books present collections of test matrices. Gregory and Karney [524, 1969] deal exclusively with the topic, while Westlake [1217, 1968] gives an appendix of test matrices.

In this chapter we present a gallery of matrices. We describe their properties and explain why they are useful (or not so useful, as the case may be) for test purposes. The coverage is limited, and all the matrices we consider are available in MATLAB. A comprehensive, well-documented collection of parametrized test matrices is contained in MATLAB's gallery function (see help gallery), and MATLAB has several other special matrices that can be used for test purposes (see help elmat).

In addition to the matrices discussed in this chapter, interesting matrices discussed elsewhere in this book include magic squares (Problem 6.4), the Kahan matrix (equation (8.11)), Hadamard matrices (§9.4), and Vandermonde matrices (Chapter 22).

Table 28.1 lists MATLAB functions that can be used to generate many of the test matrices considered in this book.

The matrices described here can be modified in various ways while still preserving some or all of their interesting properties. Among the many ways of constructing new test matrices from old are

• similarity transformations A ← X⁻¹AX,

• unitary transformations A ← UAV, where U*U = V*V = I,

• Kronecker products A ← A⊗B or B⊗ A (for which MATLAB has a function kron),

• powers A ← A^k.

28.1. The Hilbert and Cauchy Matrices

The Hilbert matrix H_n ∈ ℝ^{n × n}, with elements h_ij = l/(i + j − 1), is perhaps the most famous of all test matrices. It was widely used in the 1950s and 1960s for testing inversion algorithms and linear equation solvers. Its attractions were threefold: it is very ill conditioned for even moderate values of n, formulae are known for the elements of the inverse, and the matrix arises in a practical problem: least squares fitting by a polynomial expressed in the monomial basis.

Despite its past popularity and notoriety, the Hilbert matrix is not a good test matrix. It is too special. Not only is it symmetric positive definite, but it is totally positive.

1.1. Since E∼ _rel( x^ ) = E_rel( x^ )|x|/| x^ | and x^ = x(1 ± E_rel( x^ )), we have Erel ( x ^ ) 1+ Erel ( x ^ ) ≤ E∼ rel ( x ^ ) ≤ Erel ( x ^ ) 1− Erel ( x ^ ) . Hence if E_rel( x^ ) < 0.01, say, then there is no difference between E_rel and E∼ ∼_rel for practical purposes.

1.2. Nothing can be concluded about the last digit before the decimal point. Evaluating y to higher precision yields  t     y35262537412640768743.9999999999992500740262537412640768743.9999999999992500725972 This shows that the required digit is, in fact, 3. The interesting fact that y is so close to an integer was pointed out by Lehmer [779, 1943], who explains its connection with number theory. It is known that y is irrational [1090, 1991].

1.3. 1. x/ ( x+1 +1) .

2. 2 sin x−y 2 cos x+y 2 .

3. (x − y)(x + y). Cancellation has not been avoided, but it is now harmless if x and y are known exactly (see also Problem 3.8).

4. sin x/(l + cos x).

5. c= ((a − b)² + ab(2sinθ/2)²)^1/2.

1.4. a + ib = (x + iy)² = x² − y² + 2ixy, so b = 2xy and a = x² − y², giving x² − b²/(4x²) = a, or 4x⁴ − 4ax² − b² = 0. Hence x2 = 4a± 16 a2 +16 b2 8 = a± a2 + b2 2 . If a ≥ 0 we use this formula with the plus sign, since x² ≥ 0. If a < 0 this formula is potentially unstable, so we use the rewritten form x2 = b2 2 (−a+ a2 + b2 ) . In either case we get two values for x and recover y from y = b/(2x). Note that there are other issues to consider here, such as scaling to avoid overflow.

Caveat receptor … Anything free comes with no guarantee!

— JACK DONGARRA and ERIC GROSSE, Netlib mail header

In this appendix we provide information on how to acquire software mentioned in the book. First, we describe some basic aspects of the Internet.

B.1. Internet

A huge variety of information and software is available over the Internet, the worldwide combination of interconnected computer networks. The location of a particular object is specified by a URL, which stands for “Uniform Resource Locator”. Examples of URLs are

http://www.netlib.org/index.html ftp://ftp.netlib.org

The first example specifies a World Wide Web server (http = hypertext transfer protocol) together with a file in hypertext format (html = hypertext markup language), while the second specifies an anonymous ftp site. In any URL, the site address may, optionally, be followed by a filename that specifies a particular file. For more details about the Internet see on-line information, or one of the many books on the subject.

B.2. Netlib

Netlib is a repository of freely available mathematical software, documents, and databases of interest to the scientific computing community [350, 1987], [168, 1994]. It includes

• research codes,

• golden oldies (classic programs that are not available in standard libraries),

• the collected algorithms of the ACM,

• program libraries such as EISPACK, LINPACK, LAPACK, and MINPACK.

• back issues of NA-Digest, a weekly digest for the numerical analysis community,

• databases of conferences and performance data for a wide variety of machines.

Netlib also enables the user to download technical reports from certain institutions, to download software and errata for textbooks, and to search a “white pages” database.

Netlib can be accessed in several ways.

Since the programming is likely to be the main bottleneck in the use of an electronic computer we have given a good deal of thought to the preparation of standard routines of considerable generality for the more important processes involved in computation. By this means we hope to reduce the time taken to code up large-scale computing problems, by building them up, as it were, from prefabricated units.

— J. H. WILKINSON, The Automatic Computing Engine at the National Physical Laboratory (1948)

In spite of the self-contained nature of the linear algebra field, experience has shown that even here the preparation of a fully tested set of algorithms is a far greater task than had been anticipated.

— J. H. WILKINSON and C. REINSCH, Handbook for Automatic Computation: Linear Algebra (1971)

In this appendix we briefly describe some of the freely available program libraries that have been mentioned in this book. These packages are all available from netlib (see §B.2).

C. 1. Basic Linear Algebra Subprograms

The Basic Linear Algebra Subprograms (BLAS) are sets of Fortran primitives for matrix and vector operations. They cover all the common operations in linear algebra. There are three levels, corresponding to the types of object operated upon. In the examples below, x and y are vectors, A, B, C are rectangular matrices, and T is a square triangular matrix. Names of BLAS routines are given in parentheses. The leading “x” denotes the Fortran data type, whose possible values are some or all of

S real

D double precision

C complex

Z complex*16, or double complex

Level 1: [774, 1979] Vector operations. Inner product: x^Ty (xdot); y ← αx + y (xAXPY); vector 2-norm (y^Ty)^1/2 (xNRM2); swap vectors x ↔ y (xSWAP); scale a vector x ← αx (xSCAL); and other operations.

Level 2: [347, 1988], [348, 1988] Matrix-vector operations. Matrix times vector (gaxpy): y ← αAx + βy (xGEMV); rank-1 update: A ← A + αxy^T (xGERM); triangular solve: x ← T⁻¹x (xTRSV); and variations on these.

Level 3: [342, 1990], [343, 1990] Matrix-matrix operations. Matrix multiplication: C ← αAB + βC (xGEMM); multiple right-hand side triangular solve: A ← αT⁻¹A (xTRSM); rank-r and rank-2r updates of a symmetric matrix (xSYRK, xSYR2K); and variations on these.

The BLAS are intended to be called in innermost loops of linear algebra codes. Usually, most of the computation in a code that uses BLAS calls is done inside these calls. LINPACK [341, 1979] uses the level-1 BLAS throughout (model Fortran implementations of the level-1 BLAS are listed in [341, 1979, App. D]). LAPACK [20, 1999] exploits all three levels, using the highest possible level at all times.

Each set of BLAS comprises a set of subprogram specifications only. The specifications define the parameters to each routine and state what the routine must do, but not how it must do it. Thus the implementor has freedom over the precise implementation details (loop orderings, block algorithms, special code for special cases) and even the method (fast versus conventional matrix multiply), but the implementation is required to be numerically stable, and code that tests the numerical stability is provided with the model implementations [343, 1990], [348, 1988].

For more details on the BLAS and the advantages of using them, see the defining papers listed above, or, for example, [349, 1998] or [509, 1996, Chap. 1].

The Matrix Computation Toolbox is a collection of MATLAB M-files containing functions for constructing test matrices, computing matrix factorizations, visualizing matrices, and carrying out direct search optimization. Various other miscellaneous functions are also included.

This toolbox supersedes the Test Matrix Toolbox [606, 1995] that was described in the first edition of this book. Most of the test matrices in that toolbox have been incorporated into MATLAB in the gallery function. The new toolbox incorporates some of the other routines in the Test Matrix Toolbox and adds several new ones.

The back matter includes bibliography, name index, and subject index