Library of Congress CataloginginPublication Data Higham, Nicholas J., 1961Accuracy and stability of numerical algorithms / Nicholas J. Higham. p. cm. Includes bibliographical references (p.  ) and index. ISBN O8987 l3552 (pbk.) 1. Numerical analysisData processing. 2. Computer algorithms. I. Title. QA297.H53 1996 5 19.4’0285’5 1 dc20 9539903
o is a registered trademark.
Dedicated to Alan M. Turing and James H. Wilkinson
14 Condition Number Estimation 14.1 H o w t o E s t i m a t e C o m p o n e n t w i s e C o n d i t i o n Numbers . . . . . . . . . . . . . . . . . . . . . . . . . 14.2 The pNorm Power Method. . . . . . . . . . . . . . . 14.3 LAPACK lNorm Estimator . . . . . . . . . . . . . . 14.4 Other Condition Estimators . . . . . . . . . . . . . . 14.5 Condition Numbers of Tridiagonal Matrices . . . . . 14.6 Notes and References . . . . . . . . . . . . . . . . . . 14.6.1 LAPACK . . . . . . . . . . . . . . . . . . . . . Problems . . . . . . . . . . . . . . . . . . . . . . . . .
579 B S i n g u l a r V a l u e D e c o m p o s i t i o n , MMatrices B.1 Singular Value Decomposition . . . . . . . . . . . . . . . . . . 580 B.2 MMatrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 580 C Acquiring Software C.1 Internet. . . . . . C.2 Netlib . . . . . . . C.3 MATLAB . . . . . C.4 NAG Library and
13.1 Residuals for inverses computed by MATLAB'S INV function. . . 264 14.1 Underestimation ratio for Algorithm 14.4 for 5?5 matrix A( ?). 297 16.1 Forward and backward errors for SOR iteration. 17.1 17.2 17.3 17.4 17.5 17.6
. . . . . . . . 327
A typical hump for a convergent, nonnormal matrix. Diverging powers of a nilpotent matrix, C1 4 . . . . . Infinity norms of powers of 2 ? 2 matrix J in (17.2). Computed powers of chebspec matrices. . . . . . . . Pseudospectra for chebspec matrices. . . . . . . . . Pseudospectrum for SOR iteration matrix. . . . . . . xvii
. . . . . . . . . . . . . . 363 18.1 Householder matrix P times vector . . . . . . . . . . . . . . . . . . 372 18.2 Givens rotation, y = G(i,j,?) 22.1 Exponent versus time for matrix multiplication. . . . . . . . . . 449 22.2 Errors for Strassen’s method with two random matrices of dimension n = 1024. . . . . . . . . . . . . . . . . . . . . . . . . . 457 23.1 Error in FFT followed by inverse FFT. . . . . . . . . . . . . . . 468 24.1 The possible steps in one iteration of the MDS method when n=2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 478 25.1 Rational function r . . . . . . . . . . . . . . . . . . . . . . . . 493 25.2 Error in evaluating rational function r . . . . . . . . . . . . . . 494 26.1 spy(rem(pascal(32),2)). . . . . . . . . . . . . . . . . . . . . 524 26.2 Pseudospectra of compan(A). . . . . . . . . . . . . . . . . , . . 526 26.3 Pseudospectra of 32 ? 32 pentadiagonal Toeplitz matrices. . . . 528
List of Tables 1.1 Computed approximations = ?l((1+1/n) n ) to e = 2.71828 . . . . 16 1.2 Computed values of from Algorithms 1 and 2. . . . 23 1.3 Results from GE without pivoting on an upper Hessenberg mat r i x . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.1 2.2 2.3 2.4 2.5
Times for solution of a linear system of order n . . . . . . . . . 189 Records for largest dense linear systems solved. . . . . . . . . . 199
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
41 46 54 55 55
. . . . . . . . . . . . . .
99
. . . . . . . . 121 . . . . . 122
11.1 ?A,b values for A = orthog(25).
. . . . . . . . . . . . . . . . 240
11.2 ?A,b values for A = clement(50) . . . . . . . . . . . . . . . . . 241 11.3 ?A,b values for A = gfpp(50) . . . . . . . . . . . . . . . . . . . 241 12.1 Stability of block and point LU factorization. . . . . . . . . . . 256 13.1 13.2 13.3 13.4 13.5 13.6
Backward errors for the norm. . . . . . . . . . . . . Mflop rates for inverting a triangular matrix on a Cray 2. . . . Mflop rates for inverting a full matrix on a Cray 2. . . . . . . . Times (minutes and seconds) for inverting an n ? n matrix. . . Additional timings for inverting an n ? n matrix. . . . . . . . . GaussJordan elimination for U =6. . . . . . . . . . . . . . .
16.1 Dates of publication of selected iterative methods. xix
262 270 275 276 276 279
. . . . . . . 326
LIST OF TABLES 16.2 Result for Jacobi method, ? = ?8 j . . . . . . . . . . . . . 333 16.3 Results for Jacobi method, a = (?8 j ). . . . . . . . . . . 333 19.1 LS backward errors and residual for Vandermonde system. . . . 405 20.1 Backward errors for underdetermined Vandermonde system. . . 422 21.1 Bounds and estimates for . . . . . . . . . . . . . . . . 428 21.2 Parameters in the threeterm recurrence (21.6). . . . . . . . . . 433 21.3 Results for dual ChebyshevVandermondelike system. . . . . . 438 25.1 Results from Cholesky factorization. . . . . . . . . . . . . . . . 496 25.2 Effect of extended run time on Patriot missile operation. . . . . 507 26.1 Condition numbers of Hilbert and Pascal matrices. . . . . . . . 516
Preface It has been 30 years since the publication of Wilkinson’s books Rounding Errors in Algebraic Processes [1088, 1963] and The Algebraic Eigenvalue Problem [1089, 1965]. These books provided the first thorough analysis of the effects of rounding errors on numerical algorithms, and they rapidly became highly influential classics in numerical analysis. Although a number of more recent books have included analysis of rounding errors, none has treated the subject in the same depth as Wilkinson. This book gives a thorough, uptodate treatment of the behaviour of numerical algorithms in finite precision arithmetic. It combines algorithmic derivations, perturbation theory, and rounding error analysis. Software practicalities are emphasized throughout, with particular reference to LAPACK. The best available error bounds, some of them new, are presented in a unified format with a minimum of jargon. Historical perspective is given to provide insight into the development of the subject, and further information is provided in the many quotations. Perturbation theory is treated in detail, because of its central role in revealing problem sensitivity and providing error bounds. The book is unique in that algorithmic derivations and motivation are given succinctly, and implementation details minimized, so that attention can be concentrated on accuracy and stability results. The book was designed to be a comprehensive reference and contains extensive citations to the research literature. Although the book’s main audience is specialists in numerical analysis, it will be of use to all computational scientists and engineers who are concerned about the accuracy of their results. Much of the book can be understood with only a basic grounding in numerical analysis and linear algebra. This first two chapters are very general. Chapter 1 describes fundamental concepts of finite precision arithmetic, giving many examples for illustration and dispelling some misconceptions. Chapter 2 gives a thorough treatment of floating point arithmetic and may well be the single most useful chapter in the book. In addition to describing models of floating point arithmetic and the IEEE standard, it explains how to exploit “lowlevel” features not represented in the models and contains a large set of informative exercises. In the rest of the book the focus is, inevitably, on numerical linear algebra, because it is in this area that rounding errors are most influential and have xxi
xxii
PREFACE
been most extensively studied. However, I found that it was impossible to cover the whole of numerical linear algebra in a single volume. The main omission is the area of eigenvalue and singular value computations, which is still the subject of intensive research and requires a book of its own to summarize algorithms, perturbation theory, and error analysis. This book is therefore certainly not a replacement for The Algebraic Eigenvalue Problem. Two reasons why rounding error analysis can be hard to understand are that, first, there is no standard notation and, second, error analyses are often cluttered with rederivations of standard results. In this book I have used notation that I find nearly always to be the most convenient for error analysis: the key ingredient is the symbol ? n = nu/(1  nu), explained in §3.1. I have also summarized many basic error analysis results (for example, in Chapters 3 and 8) and made use of them throughout the book. I like to think of these basic results as analogues of the Fortran BLAS (Basic Linear Algebra Subprograms): once available in a standard form they can be used as black boxes and need not be reinvented. A number of the topics included here have not been treated in depth in previous numerical analysis textbooks. These include floating point summation, block LU factorization, condition number estimation. the Sylvester equation, powers of matrices. finite precision behaviour of stationary iterative methods, Vandermonde systems, and fast matrix multiplication, each of which has its own chapter. But there are also some notable omissions. I would have liked to include a chapter on Toeplitz systems, but this is an area in which stability and accuracy are incompletely understood and where knowledge of the underlying applications is required to guide the investigation. The important problems of updating and downdating matrix factorizations when the matrix undergoes a “small” change have also been omitted due to lack of time and space. A further omission is analysis of parallel algorithms for all the problems considered in the book (though blocked and partitioned algorithms and one particular parallel method for triangular systems are treated). Again, there are relatively few results and this is an area of active research. Throughout the history of numerical linear algebra, theoretical advances have gone hand in hand with software development. This tradition has continued with LAPACK (1987), a project to develop a stateoftheart Fortran package for solving linear equations and eigenvalue problems. LAPACK has enjoyed a synergy with research that has led to a number of important breakthroughs in the design and analysis of algorithms, from the standpoints of both performance and accuracy. A key feature of this book is that it provides the material needed to understand the numerical properties of many of the algorithms in LAPACK, the except ions being the routines for eigenvalue and singular value problems. In particular, the error bounds computed by the LAPACK linear equation solvers are explained, the LAPACK condition estimator is described in detail, and some of the software issues confronted by
XX111 the LAPACK developers are highlighted. Chapter 25 examines the influence of floating point arithmetic on general numerical software, offering salutary stories, useful techniques, and brief descriptions of relevant codes. This book has been written with numerical analysis courses in mind, although it is not designed specifically as a textbook. It would be a suitable reference for an advanced course (for example, for a graduate course on Numerical Linear Algebra following the syllabus recommended by the ILAS Education Committee [601, 1993]), and could be used by instructors at all levels as a supplementary text from which to draw examples, historical perspective, statements of results, and exercises. The exercises (actually labelled “problems”) are an important part of the book, and many of them have not, to my knowledge, appeared in textbooks before. Where appropriate I have indicated the source of an exercise; a name without a citation means that the exercise came from private communication or unpublished notes. Research problems given at the end of some sets of exercises emphasize that most of the areas covered are still active. In addition to surveying and unifying existing results (including some that have not appeared in the mainstream literature) and sometimes improving upon their presentation or proof, this book contains new results. Some of particular note are as follows. 1. The error analysis in §5.3 for evaluation of the Newton interpolating polynomial. 2. The forward error analysis for iterative refinement in §11.1. 3. The error analysis of GaussJordan elimination in §13.4. 4. The unified componentwise error analysis of QR factorization methods in Chapter 18, and the corresponding analysis of their use for solving the least squares problem in Chapter 19. 5. Theorem 20.3, which shows the backward stability of the QR factorization method for computing the minimum 2norm solution to an underdetermined system. The Notes and References are an integral part of each chapter. In addition to containing references, historical information, and further details, they include material not covered elsewhere in the chapter, and should always be consulted, in conjunction with the index, to obtain the complete picture. I have included relatively few numerical examples except in the first chapter. There are two reasons. One is to reduce the length of the book. The second reason is because today it so easy for the reader to perform experiments in MATLAB* or some other interactive system. To this end I have made *MATLAB is a registered trademark of The Math Works, Inc.
xxiv
PREFACE
available the Test Matrix Toolbox, which contains MATLAB Mfiles for many of the algorithms and special matrices described in the book; see Appendix E. This book has been designed to be as easy to use as possible. There are thorough name and subject indexes, page headings show chapter and section titles and numbers, and there is extensive crossreferencing. I have adopted the unusual policy of giving with (nearly) every citation not only its numerical location in the bibliography but also the names of the authors and the year of publication. This provides as much information as possible in a citation and reduces the need for the reader to turn to the bibliography. A database accstabnumalg.bib containing all the references in the bibliography is available over the Internet from the bibnet project (which can be accessed via netlib, described in §C.2). Special care has been taken to minimize the number of typographical and other errors, but no doubt, some remain. I will be happy to receive notification of errors, as well as comments and suggestions for improvement. Acknowledgements Three books, in addition to Wilkinson’s, have strongly influenced my research in numerical linear algebra and have provided inspiration for this book: Golub and Van Loan’s Matrix Computations [470, 1989] (first edition 1983), Parlett’s The Symmetric Eigenvalue Problem [820, 1980], and Stewart’s Introduction to Matrix Computations [941, 1973]. Knuth’s The Art of Computer Programming books [666, 19731981] have also influenced my style and presentation. Jim Demmel has contributed greatly to my understanding of the subject of this book and provided valuable technical help and suggestions. The first two chapters owe much to the work of Velvel Kahan; I am grateful to him for giving me access to unpublished notes and for suggesting improvements to early versions of Chapters 2 and 25. Des Higham read various drafts of the book, offering sound advice and finding improvements that had eluded me. Other people who have given valuable help, suggestions, or advice are Zhaojun Bai, Brad Baxter, ?ke Bj?rck, Martin CampbellKelly, Shivkumar Chandrasekaran, Alan Edelman, Warren Ferguson, Philip Gill, Gene Golub, George Hall, Sven Hammarling, Andrzej Kielbas, Philip Knight, Beresford Parlett, David Silvester, Michael Saunders, Ian Smith, Doron Swade, Nick Trefethen, Jack Williams, and Hongyuan Zha. David Carlisle provided invaluable help and advice concerning . Working with SIAM on the publication of this book was a pleasure. Special thanks go to Nancy Abbott (design), Susan Ciambrano (acquisition), Ed Cilurso (production), Beth Gallagher (copy editing), Corey Gray (production),
XXV Mary Rose Muccie (copy editing and indexing), Colleen Robishaw (design), and Sam Young (production). Research leading to this book has been supported by grants from the Engineering and Physical Sciences Research Council, by a Nuffield Science Research Fellowship from the Nuffield Foundation, and by a NATO Collaborative Research Grant held with J. W. Demmel. I was fortunate to be able to make extensive use of the libraries of the University of Manchester, the University of Dundee, Stanford University, and the University of California, Berkeley. This book was typeset in using the book document style. The references were prepared in and the index with MakeIndex. It is difficult to imagine how I could have written the book without these wonderful tools. I used the “big” software from the distribution, running on a 486DX workstation. I used text editors The Semware Editor (Semware Corporation) and GNU Emacs (Free Software Foundation) and checked spelling with PCWrite (Quicksoft).
Manchester April 1995
Nicholas J. Higham
About the Dedication This book is dedicated to the memory of two remarkable English mathematicians, James Hardy Wilkinson (19191986), FRS, and Alan Mathison Turing (19121954), FRS, both of whom made immense contributions to scientific computation. Turing’s achievements include his paper “On Computable Numbers, with an Application to the Entscheidungsproblem”, which answered Hilbert’s decidability question using the abstract device now known as a Turing machine [1025, 1936]; his work at Bletchley Park during World War II on breaking the ciphers of the Enigma machine; his 1945 report proposing a design for the Automatic Computing Engine (ACE) at the National Physical Laboratory [1026, 1945]; his 1948 paper on LU factorization and its rounding error analysis [1027, 1948]; his consideration of fundamental questions in artificial intelligence (including his proposal of the “Turing test”); and, during the last part of his life, spent at the University of Manchester, his work on morphogenesis (the development of structure and form in an organism). Turing is remembered through the Turing Award of the Association for Computing Machinery (ACM), which has been awarded yearly since 1966 [3, 1987]. For more about Turing, read the superb biography by Hodges [575, 1983], described by a reviewer as “one of the finest pieces of scholarship to appear in the history of computing” [182, 1984]. Wilkinson, like Turing a Cambridgetrained mathematician, was Turing’s assistant at the National Physical Laboratory. When Turing left, Wilkinson managed the group that built the Pilot ACE, contributing to the design and construction of the machine and its software. Subsequently, he used the machine to develop and study a variety of numerical methods. He developed backward error analysis in the 1950s and 1960s publishing the books Rounding Errors in Algebraic Processes [1088, 1963]† (REAP) and The Algebraic Eigenvalue Problem [1089, 1965]‡ (AEP), both of which rapidly achieved the status of classics. (AEP was reprinted in paperback in 1988 and, after being out of print for many years, REAP is now also available in paperback.) The AEP was described by the late Professor Leslie Fox as “almost certainly the most important and widely read title in numerical analysis”. Wilkinson also † ‡
REAP has been translated into Polish [1091, AEP has been translated into Russian [1094,
xxvii
19 6 7 ] 1970 ].
and German [1093,
19 6 9 ].
XXV111
ABOUT THE DEDICATION
contributed greatly to the development of mathematical software. The volume Handbook for Automatic Computation, Volume II: Linear Algebra [1102, 1971 ], coedited with Reinsch, contains highquality, properly documented software and has strongly influenced subsequent software projects such as the NAG Library, EISPACK, LINPACK, and LAPACK. Wilkinson received the 1970 Turing Award. In his Turing Award lecture he described life with Turing at the National Physical Laboratory in the 1940s [1096, 1971]. Wilkinson is remembered through SIAM’s James H. Wilkinson Prize in Numerical Analysis and Scientific Computing, awarded every 4 years; the Wilkinson Prize for Numerical Software, awarded by Argonne National Laboratory, the National Physical Laboratory, and the Numerical Algorithms Group; and the Wilkinson Fellowship in Scientific Computing at Argonne National Laboratory. For more about Wilkinson see the biographical memoir by Fox [403, 1987], Fox’s article [402, 1978], Parlett’s essay [821, 1990], the prologue and epilogue of the proceedings [252, 1990] of a conference held in honour of Wilkinson at the National Physical Laboratory in 1987. and the tributes in [23, 1987]. Lists of Wilkinson’s publications are given in [403, 1987] and in the special volume of the journal Linear Algebra and its Applications (88/89, April 1987) published in his memory.
Previous Home
Chapter 1 Principles of Finite Precision Computation
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 erroranalysts 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. HEWLETTPACKARD, HP15C Advanced Functions Handbook (1982)
1
Cited in Goldstine [461,
1977 ,
p. 258].
1
Next
PRINCIPLES OF FINITE PRECISION COMPUTATION
2
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 highlevel 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 subscripted lower case lower case lower case Greek
A ,B ,C ?,? letters letters ?i j , bij, cij , ?i j , letters , y, z, c, g, h letters ?, ?, ?,? , ?
ij
for for for for
matrices, matrix elements, vectors, scalars,
following the widely used convention originally introduced by Householder [587, 19 6 4 ]. The vector space of all real m ? n matrices is denoted by IR m ? n and the vector space of real nvectors by IRn . Similarly, denotes the vector space of complex m ? n matrices. Algorithms are expressed using a pseudocode based on the MATLAB language [232, 1988], [735, 1992]. Comments begin with the % symbol. Submatrices are specified with the colon notation, as used in MATLAB and Fortran 90: 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 2
For the purposes of this book an algorithm is a MATLAB program; cf. Smale [924,
1990].
1.1 NOTATION AND BACKGROUND
3
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 ?l(·), and we assume that the basic arithmetic operations op=+,,*,/ satisfy ?  < u .
(1.1)
Here, u is the unit roundoff (or machine precision), which is typically of order 108 or 1016 in single and double precision computer arithmetic, respectively, and between 1010 and 1012 on pocket calculators. For more on floating point arithmetic see Chapter 2. Computed quantities (and, in this chapter only, arbitrary approximations) wear a hat. Thus denotes the computed approximation to . Definitions are often (but not always) indicated by “:=” or “=:”, with the colon next to the object being defined. We make use of the floor and ceiling functions: is the largest integer less than or equal to , and is the smallest integer greater than or equal to . 2 The normal distribution with mean µ and variance is denoted by We measure the cost of algorithms in flops. A flop is an elementary floating point operation: +,,/, or *. We normally state only the highestorder terms of flop counts. Thus, when we say that an algorithm for n ? n matrices requires 2n 3/3 flops, we really mean 2n 3/3+O(n 2) flops. Other definitions and notation are introduced when needed. Two top ics, however, do not fit comfortably into the main text and are described in Appendix B: the singular value decomposition (SVD) and Mmatrices. All our numerical experiments were carried out either in MATLAB 4.2 [735, 1992], sometimes in conjunction with the Symbolic Math Toolbox [204, 1993], or with the Salford Software/Numerical Algorithms Group FTN903 Fortran 90 compiler, Version 1.2 [888, 1993]. Whenever we say a computation was “done in Fortran 90” we are referring to the use of this compiler. All the results quoted were obtained on a 486DX workstation, unless otherwise stated, but many of the experiments were repeated on a Sun SPARCstation, using the NAGWare4 FTN90 compiler [785, 1992]. Both machines use IEEE standard floating point arithmetic and the unit roundoff is u = 2 53 1.1 ? 10 16 in M ATLAB and in double precision in Fortran 90. (Strictly speaking, in Fortran 90 we should not use the terminology single and double precision but should refer to the appropriate KIND parameters; see, e.g., Metcalf and Reid [749, 1990, §2.6]. However, these terms are vivid and unambiguous in IEEE arithmetic, so we use them throughout the book.) 3 FTN90 is a joint trademark of Salford Software Ltd. and The Numerical Algorithms Group Ltd. 4 NAGWare is a trademark of The Numerical Algorithms Group Ltd.
PRINCIPLES OF FINITE PRECISION COMPUTATION
4
1.2. Relative Error and Significant Digits Let be an approximation to a real number the accuracy of are its absolute error
. The most useful measures of
and its relative error
(which is undefined if = 0). An equivalent definition of relative error is = (1+p). Some authors omit the absolute values = ?, where from these definitions. When the sign is important we will simply talk about “the error . In scientific computation, where answers to problems can vary enormously in magnitude, it is usually the relative error that is of interest, because it is scale independent: scaling and leaves unchanged. Relative error is connected with the notion of correct significant digits (or correct significant figures). The significant digits in a number are the first nonzero digit and all succeeding digits. Thus 1.7320 has five significant digits, while 0.0491 has only three. What is meant by correct significant digits in a number that approximates another seems intuitively clear, but a precise definition is problematic, as we explain in a moment. First, note that for a number with p significant digits there are only p +1 possible answers to the question “how many correct significant digits does have?” (assuming is not a constant such as 2.0 that is known exactly). Therefore the number of correct significant digits is a fairly crude measure of accuracy in comparison with the relative error. For example, in the following two cases agrees with to three but not four significant digits by any reasonable definition, yet the relative errors differ by a factor of about 44: = 1.00000, = 9.00000,
= 1.00499, = 8.99899,
= 4.99 ? 103, = 1.12 ? 104.
Here is a possible definition of correct significant digits: an approximation to has p correct significant digits if and round to the same number to p significant digits. Rounding is the act of replacing a given number by the nearest p significant digit number, with some rule for breaking ties when there are two nearest. This definition of correct significant digits is mathematically elegant and agrees with intuition most of the time. But consider the numbers = 0.9949,
= 0.9951.
According to the definition does not have two correct significant digits but does have one and three correct significant digits!
6
PRINCIPLES OF FINITE PRECISION COMPUTATION
or, if the data is itself the solution to another problem, it may be the result of errors in an earlier computation. The effects of errors in the data are generally easier to understand than the effects of rounding errors committed during a computation, because data errors can be analysed using perturbation theory for the problem at hand, while intermediate rounding errors require an analysis specific to the given method. This book contains perturbation theory for most of the problems considered, for example, in Chapters 7 (linear systems), 19 (the least squares problem), and 20 (underdetermined systems). Analysing truncation errors, or discretization errors, is one of the major tasks of the numerical analyst. Many standard numerical methods (for example, the trapezium rule for quadrature, Euler’s method for differential equations, and Newton’s method for nonlinear equations) can be derived by taking finitely many terms of a Taylor series. The terms omitted constitute the truncation error, and for many methods the size of this error depends on a parameter (often called h, “the stepsize”) whose appropriate value is a compromise between obtaining a small error and a fast computation. Because the emphasis of this book is on finite precision computation, with virtually no mention of truncation errors, it would be easy for the reader to gain the impression that the study of numerical methods is dominated by the study of rounding errors. This is certainly not the case. Trefethen explains it well when he discusses how to define numerical analysis [1016, 1992]: Rounding errors and instability are important, and numerical analysts will always be the experts in these subjects and at pains to ensure that the unwary are not tripped up by them. But our central mission is to compute quantities that are typically uncomputable, from an analytic point of view, and to do it with lightning speed. In this quotation “uncomputable” means that approximations are necessary, and thus Trefethen’s point is that developing good approximations is a more fundamental task than analysing the effects of rounding errors on those approximations. A possible way to avoid rounding and truncation errors (but not data errors) is to try to solve a problem using a symbolic manipulation package, such as Maple5 [199, 1991] or Mathematica6 [1109, 1991]. Indeed, we have used this approach to compute “exact answers” in some of our numerical experiments, While we acknowledge the value of symbolic manipulation as part of the toolkit of the scientific problem solver, we do not study it in this book. 5 6
Maple is a registered trademark of Waterloo Maple Software. Mathematica is a registered trademark of Wolfram Research Inc.
1.3 SOURCES OF ERRORS
5
A definition of correct significant digits that does not suffer from the latter anomaly states that agrees with to p significant digits if is less than half a unit in the pth significant digit of . However, this definition implies that 0.123 and 0.127 agree to two significant digits, whereas many people would say that they agree to less than two significant digits. In summary, while the number of correct significant digits provides a useful way in which to think about the accuracy of an approximation, the relative error is a more precise measure (and is base independent). Whenever we give an approximate answer to a problem we should aim to state an estimate or bound for the relative error. When and are vectors the relative error is most often defined with a norm, as For the commonly used norms := maxi := the inequality < ? ? 10p implies that components with have about p correct significant decimal digits, but for the smaller components the inequality merely bounds the absolute error. A relative error that puts the individual relative errors on an equal footing is the componentwise relative error
which is widely used in error analysis and perturbation analysis (see Chapter 7, for example). As an interesting aside we mention the “tablemaker’s dilemma”. Suppose you are tabulating the values of a transcendental function such as the sine function and a particular entry is evaluated as 0.124500000000 correct to a few digits in the last place shown, where the vertical bar follows the final significant digit to be tabulated. Should the final significant digit be 4 or 5? The answer depends on whether there is a nonzero trailing digit and, in principle, we may never be able answer the question by computing only a finite number of digits.
1.3. Sources of Errors There are three main sources of errors in numerical computation: rounding, data uncertainty, and truncation. Rounding errors, which are an unavoidable consequence of working in finite precision arithmetic, are largely what this book is about. The remainder of this chapter gives basic insight into rounding errors and their effects. Uncertainty in the data is always a possibility when we are solving practical problems. It may arise in several ways: from errors in measuring physical quantities, from errors in storing the data on the computer (rounding errors),
1.4 PRECISION VERSUS ACCURACY
7
1.4. Precision Versus Accuracy The terms accuracy and precision are often confused or used interchangeably, but it is worth making a distinction between them. Accuracy refers to the absolute or relative error of an approximate quantity. Precision is the accuracy with which the basic arithmetic operations +,,*,/ are performed, and for floating point arithmetic is measured by the unit roundoff u (see (1.1)). Accuracy and precision are the same for the scalar computation c = a*b, but accuracy can be much worse than precision in the solution of a linear system of equations, for example. It is important to realize that accuracy is not limited by precision, at least in theory. This may seem surprising, and may even appear to contradict many of the results in this book. However, arithmetic of a given precision can be used to simulate arithmetic of arbitrarily high precision, as explained in §25.9. (The catch is that such simulation is too expensive to be of practical use for routine computation.) In all our error analyses there is an implicit assumption that the given arithmetic is not being used to simulate arithmetic of a higher precision.
1.5. Backward and Forward Errors Suppose that an approximation to y = is computed in an arithmetic of precision u, where f is a real scalar function of a real scalar variable. How should we measure the “quality” of ? In most computations we would be happy with a tiny relative error, but this cannot always be achieved. Instead of focusing on the relative error of we can ask “for what set of data have we actually solved our problem?“, that is, for what do we have = ? In general, there may be many such so we should ask for the smallest one. The value of (or min ), possibly divided by , is called the backward error. The absolute and relative errors of are called forward errors, to distinguish them from the backward error. Figure 1.1 illustrates these concepts. The process of bounding the backward error of a computed solution is called backward error analysis, and its motivation is twofold. First, it interprets rounding errors as being equivalent to perturbations in the data. The data frequently contains uncertainties due to previous computations or errors committed in storing numbers on the computer. If the backward error is no larger than these uncertainties then the computed solution can hardly be criticizedit may be the solution we are seeking, for all we know. The second attraction of backward error analysis is that it reduces the question of bounding or estimating the forward error to perturbation theory, which for many problems is well understood (and only has to be developed once, for the
8
PRINCIPLES OF FINITE PRECISION COMPUTATION Input space
Figure 1.1. Backward and forward errors for y = line = computed.
Output space
Solid line = exact; dotted
given problem, and not for each method). We discuss perturbation theory in the next section. A method for computing y = is called backward stable if, for any it produces a computed with a small backward error, that is, = for some small . The definition of “small” will be context dependent. In general, a given problem has several methods of solution, some of which are backward stable and some not. As an example, assumption (1.1) says that the computed result of the operation ± y is the exact result for perturbed data (1 + ?) and y(1 + ?) with ? < u; thus addition and subtraction are, by assumption, backward stable operations. Most routines for computing the cosine function do not satisfy = cos( + ) with a relatively small but only the weaker relation + ?y = cos( + ), with relatively small ?y and . A result of the form (1.2) is known as a mixed forwardbackward error result and is illustrated in Figure 1.2. Provided that and are sufficiently small, (1.2) says that the computed value scarcely differs from the value + that would have been produced by an input + scarcely different from the actual input Even more simply, is almost the right answer for almost the right data. In general, an algorithm is called numerically stable if it is stable in the mixed forwardbackward error sense of (1.2) (hence a backward stable algorithm can certainly be called numerically stable). Note that this definition is specific to problems where rounding errors are the dominant form of errors. The term stability has different meanings in other areas of numerical analysis.
1.6 CONDITIONING Input space
Figure 1.2. Mixed forwardbackward error for y = line = computed.
Output space
Solid line = exact; dotted
1.6. Conditioning The relationship between forward and backward error for a problem is governed by the conditioning of the problem, that is, the sensitivity of the solution to perturbations in the data. Continuing the y = example of the previous section, let an approximate solution satisfy = Then, assuming for simplicity that f is twice continuously differentiable,
and we can bound or estimate the righthand side. This expansion leads to the notion of condition number. Since
the quantity
measures, for small the relative change in the output for a given relative change in the input, and it is called the (relative) condition number of f . If or f is a vector then the condition number is defined in a similar way using norms and it measures the maximum relative change, which is attained for some, but not all, vectors As an example, consider the function = log The condition number is = which is large for 1. This means that a small relative change in can produce a much larger relative change in log for 1. The
10
PRINCIPLES OF FINITE PRECISION COMPUTATION
reason is that a small relative change in produces a small absolute change in = log (since and that change in log may be large in a relative sense. When backward error, forward error, and the condition number are defined in a consistent fashion we have the useful rule of thumb that forward error
condition number ? backward error,
with approximate equality possible. One way to interpret this rule of thumb is to say that the computed solution to an illconditioned problem can have a large forward error. For even if the computed solution has a small backward error, this error can be amplified by a factor as large as the condition number when passing to the forward error. One further definition is useful. If a method produces answers with forward errors of similar magnitude to those produced by a backward stable method, then it is called forward stable. Such a method need not be backward stable itself. Backward stability implies forward stability, but not vice versa. An example of a method that is forward stable but not backward stable is Cramer’s rule for solving a 2 ? 2 linear system, which is discussed in §1.10.1.
1.7. Cancellation Cancellation is what happens when two nearly equal numbers are subtracted. It is often, but not always, a bad thing. Consider the function = (1 With = 1.2?105 the value of rounded to 10 significant figures is c = 0.9999 9999 99, so that 1  c = 0.0000 0000 01. 10
Then (1  c) = 10 /1.44?1010 = 0.6944. . . , which is clearly wrong given the fact that 0 < < l/2 for all 0. A 10 significant figure approximation to cos is therefore not sufficient to yield a value of with even one correct figure. The problem is that 1  c has only 1 significant figure. The subtraction 1  c is exact, but this subtraction produces a result of the same size as the error in c. In other words, the subtraction elevates the importance of the earlier error. In this particular example it is easy to rewrite to avoid the cancellation. Since cos = 1  2 sin2
Evaluating this second formula for with a 10 significant figure approximation to sin yields = 0.5, which is correct to 10 significant figures.
1.8 SOLVING A QUADRATIC EQUATION
11
To gain more insight into the cancellation phenomenon consider the subtraction (in exact arithmetic) where and The terms ?a and ?b are relative errors or uncertainties in the data, perhaps attributable to previous computations. With x = a  b we have
The relative error bound for is large when a  b > y z > 0 then the cancellation in the evaluation of x + (y  z) is harmless.
1.8. Solving a Quadratic Equation Mathematically, the problem of solving the (real) quadratic equation ax 2 bx + c = 0 is trivial: there are two roots (if a 0)) given by
+
(1.3) Numerically, the problem is more challenging, as neither the successful evaluation of (1.3) nor the accuracy of the computed roots can be taken for granted. The easiest issue to deal with is the choice of formula for computing the roots. If b2 >> 4ac then and so for one choice of sign the formula (1.3) suffers massive cancellation. This is damaging cancellation because one of the arguments, is inexact, so the subtraction brings into prominence the earlier rounding errors. How to avoid the cancellation is well known: obtain the larger root (in absolute value), x 1, from
and the other from the equation x 1 x 2 = c/a.
12
PRINCIPLES OF FINITE PRECISION COMPUTATION
Unfortunately, there is a more pernicious source of cancellation: the subtraction b 2  4ac. Accuracy is lost here when b 2 4ac (the case of nearly equal roots), and no algebraic rearrangement can avoid the cancellation. The only way to guarantee accurate computed roots is to use double precision (or some trick tantamount to the use of double precision) in the evaluation of b 2  4ac. Another potential difficulty is underflow and overflow. If we apply the formula (1.3) in IEEE single precision arithmetic (described in §2.3) to the equation 1020 x2  3·1020 x+2·1020 = 0 then overflow occurs, since the maximum floating point number is of order 10 the roots, however, are innocuous: x = 1 and x = 2. Dividing through the equation by max( a, b, c ) = 1020 cures the problem, but this strategy is ineffective for the equation 1020x 2 3x+2·1020 = 0, whose roots are 1020 and 2·1020. In the latter equation we need to scale the variable: defining x = 1020 y gives 1020 y2 3·10 20 y+ 2 · 1 0 2 0 = 0 , which is the first equation we considered. These ideas can be built into a general scaling strategy (see the Notes and References), but the details are nontrivial. As this discussion indicates, not only is it difficult to devise an accurate and robust algorithm for solving a quadratic equation, but it is a nontrivial task to prepare specifications that define precisely what “accurate” and “robust” mean for a given system of floating point arithmetic.
1.9. Computing the Sample Variance In statistics the sample variance of n numbers x 1, . . . , xn is defined as (1.4) where the sample mean
Computing from this formula requires two passes through the data, one to compute and the other to accumulate the sum of squares. A twopass computation is undesirable for large data sets or when the sample variance is to be computed as the data is generated. An alternative formula, found in many statistics textbooks, uses about the same number of operations but requires only one pass through the data: (1.5)
1.10 SOLVING LINEAR EQUATIONS
13
This formula is very poor in the presence of rounding errors because it computes the sample variance as the difference of two positive numbers, and therefore can suffer severe cancellation that leaves the computed answer dominated by roundoff. In fact, the computed answer can be negative, an event aptly described by Chan, Golub, and LeVeque [194, 1983] as “a blessing in disguise since this at least alerts the programmer that disastrous cancellation has occurred”. In contrast, the original formula (1.4) always yields a very accurate (and nonnegative) answer, unless n is large (see Problem 1.10). Surprisingly, current calculators from more than one manufacturer (but not HewlettPackard) appear to use the onepass formula, and they list it in their manuals. As an example, if x = [10000, 10001, 10002]T then, in single precision arithmetic (u 6 ? 108), the sample variance is computed as 1.0 by the twopass formula (relative error 0) but 0.0 by the onepass formula (relative error 1). It might be argued that this data should be shifted by some estimate of the mean before applying the onepass formula which does not change ), but a good estimate is not always available and there are alternative onepass algorithms that will always produce an acceptably accurate answer. For example, instead of accumulating and we can accumulate
which can be done via the updating formulae (1.6a) (1.6b) after which = Q n / (n  1). Note that the only subtractions in these recurrences are relatively harmless ones that involve the data xi . For the numerical example above, (1.6) produces the exact answer. The updating formulae (1.6) are numerically stable, though their error bound is not as small as the one for the twopass formula (it is proportional to the condition number KN in Problem 1.7). The problem of computing the sample variance illustrates well how mathematically equivalent formulae can have different numerical stability properties.
1.10. Solving Linear Equations For an approximate solution y to a linear system Ax = b ( A IRn ? n , b n IR ) the forward error is defined as xy/x, for some appropriate norm.
14
PRINCIPLES OF FINITE PRECISION COMPUTATION
Another measure of the quality of y, more or less important depending on the circumstances, is the size of the residual r = b  Ay . When the linear system comes from an interpolation problem, for example, we are probably more interested in how closely Ay represents b than in the accuracy of y. The residual is scale dependent: multiply A and b by a and r is multiplied by a. One way to obtain a scaleindependent quantity is to divide by  A y, yielding the relative residual
The importance of the relative residual is explained by the following resuit, which was probably first proved by Wilkinson (see the Notes and References). We use the 2norm, defined by x2 = (xT x)? and A2 =
Lemma 1.1. With the notation above, and for the 2norm,
Proof. If (A+?A)y = b then r := bAy = ?Ay, so r  2 <  ?A  2  y   2 , giving (l.7) On the other hand, (A +?A)y = b for ?A = ryT /(yT y) and ?A2 =  r2 /y2, so the bound (1.7) is attainable. Lemma 1.1 says that p(y) measures how much A (but not b) must be perturbed in order for y to be the exact solution to the perturbed system, that is, p(y) equals a normwise relative backward error. If the data A and b are uncertain and p(y) is no larger than this uncertainty (e.g., p(y) = O(u)) then the approximate solution y must be regarded as very satisfactory. For other problems the backward error may not be as easy to compute as it is for a general linear system, as we will see for the Sylvester equation (§15.2) and the least squares problem (§19.7). To illustrate these concepts we consider two specific linear equation solvers: Gaussian elimination with partial pivoting (GEPP) and Cramer’s rule.
1.10.1. GEPP Versus Cramer’s Rule Cramer’s rule says that the components of the solution to a linear system Ax = b are given by x i = det(A i (b))/det(A), where Ai ( b) denotes A with its ith column replaced by b. These formulae are a prime example of a method
1.10 SOLVING LINEAR EQUATIONS
15
that is mathematically elegant, but useless for solving practical problems. The two flaws in Cramer’s rule are its computational expense and its numerical instability. The computational expense needs little comment, and is, fortunately, explained in most modern linear algebra textbooks (for example, Strang [961, 1993] cautions the student “it would be crazy to solve equations that way”). The numerical instability is less well known, but not surprising. It is present even for n = 2, as a numerical example shows. We formed a 2 ? 2 system Ax = b with condition number K2 (A) = A 2 A 12 1013, and solved the system by both Cramer’s rule and GEPP in MATLAB (unit roundoff u 1.1 ? 10 16 ). The results were as follows, where r = b  A : Cramer’s rule 1.0000 2.0001
1.5075 ? 107 1.9285 ? 107
GEPP 1.0002 2.0004
4.5689 ? 10 17 2.1931 ? 10 17
The scaled residual for GEPP is pleasantly smallof order the unit roundoff. That for Cramer’s rule is ten orders of magnitude larger, showing that the computed solution from Cramer’s rule does not closely satisfy the equations, or, equivalently, does not solve a nearby system. The solutions themselves are similar, both being accurate to three significant figures in each component but incorrect in the fourth significant figure. This is the accuracy we would expect from GEPP because of the rule of thumb “forward error backward error ? condition number”. That Cramer’s rule is as accurate as GEPP in this example, despite its large residual, is perhaps surprising, but it is explained by the fact that Cramer’s rule is forward stable for n = 2; see Problem 1.9. For general n, the accuracy and stability of Cramer’s rule depend on the method used to evaluate the determinants, and satisfactory bounds are not known even for the case where the determinants are evaluated by GEPP. The small residual produced by GEPP in this example is typical: error analysis shows that GEPP is guaranteed to produce a relative residual of order u when n = 2 (see §9.2). To see how remarkable a property this is, consider the rounded version of the exact solution: z = fl(x) = x + ?x, where ?x 2 1,
which rounds to 1, since the HP 48G works to about 12 decimal digits. Thus for x > 1, the repeated square roots reduce x to 1.0, which the squarings leave unchanged. For 0 < x < 1 we have
on a 12digit calculator, so we would expect the square root to satisfy
This upper bound rounds to the 12 significant digit number 0.99. . .9. Hence after the 60 square roots we have on the calculator a number x < 0.99. . .9. The 60 squarings are represented by s(x) = and
Because it is smaller than the smallest positive representable number, this result is set to zero on the calculatora process known as underflow. (The converse situation, in which a result exceeds the largest representable number, is called overflow.) The conclusion is that there is nothing wrong with the calculator. This innocuouslooking calculation simply exhausts the precision and range of a machine with 12 digits of precision and a 3digit exponent. 1.12.3. An Infinite Sum It is well known that = 1.6449 3406 6848. . . . Suppose we were not aware of this identity and wished to approximate the sum numerically. The most obvious strategy is to evaluate the sum for increasing k until
1.13 INCREASING THE PRECISION
19
the computed sum does not change. In Fortran 90 in single precision this yields the value 1.6447 2532, which is first attained at k = 4096. This agrees with the exact infinite sum to just four significant digits out of a possible nine. The explanation for the poor accuracy is that we are summing the numbers from largest to smallest, and the small numbers are unable to contribute to the sum. For k = 4096 we are forming s + 40962 = s + 224, where s 1.6. Single precision corresponds to a 24bit mantissa, so the term we are adding to s “drops off the end” of the computer word, as do all successive terms. The simplest cure for this inaccuracy is to sum in the opposite order: from smallest to largest. Unfortunately, this requires knowledge of how many terms to take before the summation begins. With 109 terms we obtain the computed sum 1.6449 3406, which is correct to eight significant digits. For much more on summation, see Chapter 4.
1.13. Increasing the Precision When the only source of errors is rounding, a common technique for estimating the accuracy of an answer is to recompute it at a higher precision and to see how many digits of the original and the (presumably) more accurate answer agree. We would intuitively expect any desired accuracy to be achievable by computing at a high enough precision. This is certainly the case for algorithms possessing an error bound proportional to the precision, which includes all the algorithms described in the subsequent chapters of this book. However, since an error bound is not necessarily attained, there is no guarantee that a result computed in t digit precision will be more accurate than one computed in s digit precision, for a given t > s; in particular, for a very ill conditioned problem both results could have no correct digits. For illustration, consider the system Ax = b, where A is the inverse of the 5?5 Hilbert matrix and b i = (l)i i. (For details of the matrices used in this experiment see Chapter 26.) We solved the system in varying precisions with unit roundoffs u = 2 t , t = 15:40, corresponding to about 4 to 12 decimal places of accuracy. (This was accomplished in MATLAB by using the function chop from the Test Matrix Toolbox to round the result of every arithmetic operation to t bits; see Appendix E.) The algorithm used was Gaussian elimination (without pivoting), which is perfectly stable for this symmetric positive definite matrix. The upper plot of Figure 1.3 shows t against the relative errors and the relative residuals b The lower plot of Figure 1.3 gives corresponding results for A = P5 + 5I, where P5 is the Pascal matrix of order 5. The condition numbers (A) are 1.62?102 for the inverse Hilbert matrix and 9.55?105 for the shifted Pascal matrix. In both cases the general trend is that increasing the precision decreases the residual and relative error, but the behaviour is
20
P RINCIPLES
OF
FINITE P RECISION C OMPUTATION
Figure 1.3. Forward errors and relative residuals b versus precision t =  log2 u on the x axis. not monotonic. The reason for the pronounced oscillating behaviour of the relative error (but not the residual) for the inverse Hilbert matrix is not clear. An example in which increasing the precision by several bits does not improve the accuracy is the evaluation of y = x + a sin(bx),
a = 108,
b = 224.
(1.8)
Figure 1.4 plots t versus the absolute error, for precisions u = 2t, t = 10:40. Since a sin(bx) 8.55?109, for t less than about 20 the error is dominated by the error in representing x = l/7. For 22 < t < 31 the accuracy is (exactly) constant! The plateau over the range 22 < t < 31 is caused by a fortuitous rounding error in the addition: in the binary representation of the exact answer the 23rd to 32nd digits are 1s, and in the range of t of interest the final rounding produces a number with a 1 in the 22nd bit and zeros beyond, yielding an unexpectedly small error that affects only bits 33rd onwards. A more contrived example in which increasing the precision has no beneficial effect on the accuracy is the following evaluation of z = f(x): y = abs(3(x0.5)0.5)/25 if y = 0 z = l
1.14 C ANCELLATION
OF
R OUNDING E RRORS
21
Figure 1.4. Absolute error versus precision, t = log2 u. else z = ey % Store to inhibit extended precision evaluation. z = (z  1)/y end In exact arithmetic, z = f(2/3) = 1, but in Fortran 90 on a Sun SPARCstation and on a 486DX workstation, = fl(f(2/3)) = 0.0 in both single and double precision arithmetic. A further example is provided by the “innocuous calculation” of §1.12.2, in which a step function is computed in place of f(x) = x for a wide range of precisions. It is worth stressing that how precision is increased can greatly affect the results obtained. Increasing the precision without preserving important properties such as monotonicity of rounding can vitiate an otherwise reliable algorithm. Increasing the precision without maintaining a correct relationship among the precisions in different parts of an algorithm can also be harmful to the accuracy.
1.14. Cancellation of Rounding Errors It is not unusual for rounding errors to cancel in stable algorithms, with the result that the final computed answer is much more accurate than the inter
P RINCIPLES
22
OF
FINITE P RECISION C OMPUTATION
mediate quantities. This phenomenon is not universally appreciated, perhaps because we tend to look at the intermediate numbers in an algorithm only when something is wrong, not when the computed answer is sat isfactory. We describe two examples. The first is a very short and rather unusual computation, while the second involves a wellknown algorithm for computing a standard matrix decomposition. 1.14.1. Computing (ex

1)/x
Consider the function f(x) = (ex  1)/x = which arises in various applications. The obvious way to evaluate f is via the algorithm % Algorithm 1. if x = 0 f = l else f = (ex  1)/x end This algorithm suffers severe cancellation for x 