《数值分析》课程教学资源（课外阅读）Numerical Analysis（Trefethen, Princeton Companion to Mathematics, 2008）

604 IV.Branches of Mathematics really goes on isa far more interesting process ofx be done even in principle by formulas,for most mathe. Press/ wnat na ppens instead i answers that are accurate to three or ten digits of pre aundhedorasctentfncoremgineing tion.Provi approximate solutions by an eler entary example.Sup. pery te pose we have one polynomial of degree 4. p(z)=C0+C1z+C222+C3z3+C4z hetineplherconstr and another of degree5, 767-9 4(2=d0+d1z+d2z2+d3z3+d4z+d5z5 on to the th It is well-known that there is an explicit formula that V.1990:Algebraic c expresses the roots of p in terms of radicals (discov bridge,MA:MIT Press/Elsevier more than 250 years later;see THE OF THE QUINTIC I for r )Thus e in a ce Pand are.Yetn practice they hardly A N.Trefethen differ at all.If a scientist or a mathematician ants 1 The Need for Numerical Computation digits of precision in less than a millisecond.Did the nth they turn to computers. Nev ertheless,there is a wide maybe not.Most of the time.the user neither knows The It is otion w made it a dtredcou th and others princip Here are three more examples of problems that can ments led to physical laws exp ssed mathematically nce of elementary and,in the remarkable cycle 05 fruits are all aroun operat i)Lin ments.The day has long sinc tions in n unknowns. :solve a system of n linear equa passed when an advanc in the physical da signif (Linear minimize a linear function roblem:find the shortest tour mathematics. between n cities. nle imagine that mathemati And here are five that,like finding the roots of q,cannot generally be solved in this manner cians generate formulas,and then,by inserting num bers intot essary results h

✐ 604 IV. Branches of Mathematics Karp, R. M., and V. Ramachandran. 1990. Parallel algorithms for shared-memory machines. In Handbook of Theoret￾ical Computer Science, volume A, Algorithms and Com￾plexity, edited by J. van Leeuwen. Cambridge, MA: MIT Press/Elsevier. Kearns, M. J., and U. V. Vazirani. 1994. An Introduction to Computational Learning Theory. Cambridge, MA: MIT Press. Kitaev, A., A. Shen, and M. Vyalyi. 2002. Classical and Quan￾tum Computation. Providence, RI: American Mathematical Society. Kushilevitz, E., and N. Nisan. 1996. Communication Com￾plexity. Cambridge: Cambridge University Press. Ron, D. 2001. Property testing (a tutorial). In Handbook on Randomized Computing, volume II. Dordrecht: Kluwer. Shaltiel, R. 2002. Recent developments in explicit construc￾tions of extractors. Bulletin of the European Association for Theoretical Computer Science 77:67–95. Sipser, M. 1997. Introduction to the Theory of Computation. Boston, MA: PWS. Strassen, V. 1990: Algebraic complexity theory. In Hand￾book of Theoretical Computer Science, volume A, Algo￾rithms and Complexity, edited by J. van Leeuwen. Cam￾bridge, MA: MIT Press/Elsevier. IV.21 Numerical Analysis Lloyd N. Trefethen 1 The Need for Numerical Computation Everyone knows that when scientists and engineers need numerical answers to mathematical problems, they turn to computers. Nevertheless, there is a wide￾spread misconception about this process. The power of numbers has been extraordinary. It is often noted that the scientific revolution was set in motion when Galileo and others made it a principle that everything must be measured. Numerical measure￾ments led to physical laws expressed mathematically, and, in the remarkable cycle whose fruits are all around us, finer measurements led to refined laws, which in turn led to better technology and still finer measure￾ments. The day has long since passed when an advance in the physical sciences could be achieved, or a signifi- cant engineering product developed, without numerical mathematics. Computers certainly play a part in this story, yet there is a misunderstanding about what their role is. Many people imagine that scientists and mathemati￾cians generate formulas, and then, by inserting num￾bers into these formulas, computers grind out the nec￾essary results. The reality is nothing like this. What really goes on is a far more interesting process of exe￾cution of algorithms. In most cases the job could not be done even in principle by formulas, for most mathe￾matical problems cannot be solved by a finite sequence of elementary operations. What happens instead is that fast algorithms quickly converge to “approximate” answers that are accurate to three or ten digits of pre￾cision, or a hundred. For a scientific or engineering application, such an answer may be as good as exact. We can illustrate the complexities of exact versus approximate solutions by an elementary example. Sup￾pose we have one polynomial of degree 4, p(z) = c0 + c1z + c2z2 + c3z3 + c4z4, and another of degree 5, q(z) = d0 + d1z + d2z2 + d3z3 + d4z4 + d5z5. It is well-known that there is an explicit formula that expresses the roots of p in terms of radicals (discov￾ered by Ferrari around 1540), but no such formula for the roots of q (as shown by Ruffini and abel [VI.33] more than 250 years later; see the insolubility of the quintic [V.21] for more details). Thus, in a cer￾tain philosophical sense the root-finding problems for p and q are utterly different. Yet in practice they hardly differ at all. If a scientist or a mathematician wants to know the roots of one of these polynomials, he or she will turn to a computer and get an answer to sixteen digits of precision in less than a millisecond. Did the computer use an explicit formula? In the case of q, the answer is certainly no, but what about p? Maybe, maybe not. Most of the time, the user neither knows nor cares, and probably not one mathematician in a hundred could write down formulas for the roots of p from memory. Here are three more examples of problems that can be solved in principle by a finite sequence of elementary operations, like finding the roots of p. (i) Linear equations: solve a system of n linear equa￾tions in n unknowns. (ii) Linear programming: minimize a linear function of n variables subject to m linear constraints. (iii) Traveling salesman problem: find the shortest tour between n cities. And here are five that, like finding the roots of q, cannot generally be solved in this manner. (iv) Find an eigenvalue [I.3 §4.3] of an n × n matrix. (v) Minimize a function of several variables

IV.21.Numerical Analysis 605 (vi)Evaluate an integral ics began to explode,but now mainly in the hands of specialists.New journals ed suc matik (959).The revolution was sparked by hardwar Can we conclude that (1)-(l)willbe easier than (v)-(vii but it included mathemati l and algorithmic develo in practice?Absolutely not.Problem (iii)is usually very nn is,say,m the nundr if the integral is in one dimension.Problems()and (v) of around 10,but so did the best algorithms known are of al 100.and easy wh increase in ik 1000000.m fact.n these matters Half a century on,numerical analysis has grown in guide to practice that,for each of the three to one of the largest branches of mathemati cs,the resea ignores the exact solution and uses appr ions journals across the sciences and engine fast!)methods instead g the pysis hanks to the rs of these back man of ers.we have reached a point where most of the clas ariables.(Thi defin lems of the physical sciences n inc osed over the real numbers.but not their discrete ana this ogues.)In the remain ona strong foundation rends This field encomnasses questions of inter 2 A Brief History polatio series expansi ons,and HARMO ANALY a5s00 wit NE Throughout history.leading mathematicians have been with,and in many case polyn omial and rational minimax appro imation asso sul in use today ed with n mes such as CHEE [VL45]and B standing example Among many other contributions basis functions.and wAVELETS.31.We shall not 179 in lea ata fitting have space to address these subjects,bu in almost ical quadrature (1814),as well as inventing THE FAST HOURIER TRANS 6()though the Cooley and Tukey in 1965. Around 1900,the numerical side of mathematics 3 Machine Arithmetic and Rounding Errors I is well-kno n that natics generally and of great uated on a computer,fo hines to work in hase f we desg xample.many advan of the early twentieth cen eal numbers by a system of floating-point arithmetic ne in whi each mber is in a digita numerical calculation matter unless the numberis s huge or tin as to cause or underflow Floating-point arithmetic wa invented by Konrac Zuse in Berlin in the 1930s,and

✐ IV.21. Numerical Analysis 605 (vi) Evaluate an integral. (vii) Solve an ordinary differential equation (ODE). (viii) Solve a partial differential equation (PDE). Can we conclude that (i)–(iii) will be easier than (iv)–(viii) in practice? Absolutely not. Problem (iii) is usually very hard indeed if n is, say, in the hundreds or thousands. Problems (vi) and (vii) are usually rather easy, at least if the integral is in one dimension. Problems (i) and (iv) are of almost exactly the same difficulty: easy when n is small, like 100, and often very hard when n is large, like 1 000 000. In fact, in these matters philosophy is such a poor guide to practice that, for each of the three problems (i)–(iii), when n and m are large one often ignores the exact solution and uses approximate (but fast!) methods instead. Numerical analysis is the study of algorithms for solving the problems of continuous mathematics, by which we mean problems involving real or complex variables. (This definition includes problems like lin￾ear programming and the traveling salesman problem posed over the real numbers, but not their discrete ana￾logues.) In the remainder of this article we shall review some of its main branches, past accomplishments, and possible future trends. 2 A Brief History Throughout history, leading mathematicians have been involved with scientific applications, and in many cases this has led to the discovery of numerical algorithms still in use today. gauss [VI.26], as usual, is an out￾standing example. Among many other contributions, he made crucial advances in least-squares data fitting (1795), systems of linear equations (1809), and numer￾ical quadrature (1814), as well as inventing the fast fourier transform [III.26] (1805), though the last did not become widely known until its rediscovery by Cooley and Tukey in 1965. Around 1900, the numerical side of mathematics started to become less conspicuous in the activities of research mathematicians. This was a consequence of the growth of mathematics generally and of great advances in fields in which, for technical reasons, math￾ematical rigor had to be the heart of the matter. For example, many advances of the early twentieth cen￾tury sprang from mathematicians’ new ability to reason rigorously about infinity, a subject relatively far from numerical calculation. A generation passed, and in the 1940s the computer was invented. From this moment numerical mathemat￾ics began to explode, but now mainly in the hands of specialists. New journals were founded such as Math￾ematics of Computation (1943) and Numerische Mathe￾matik (1959). The revolution was sparked by hardware, but it included mathematical and algorithmic develop￾ments that had nothing to do with hardware. In the half￾century from the 1950s, machines sped up by a factor of around 109, but so did the best algorithms known for some problems, generating a combined increase in speed of almost incomprehensible scale. Half a century on, numerical analysis has grown in￾to one of the largest branches of mathematics, the specialty of thousands of researchers who publish in dozens of mathematical journals as well as applica￾tions journals across the sciences and engineering. Thanks to the efforts of these people going back many decades, and thanks to ever more powerful comput￾ers, we have reached a point where most of the clas￾sical mathematical problems of the physical sciences can be solved numerically to high accuracy. Most of the algorithms that make this possible were invented since 1950. Numerical analysis is built on a strong foundation: the mathematical subject of approximation theory. This field encompasses classical questions of inter￾polation, series expansions, and harmonic analy￾sis [IV.11] associated with newton [VI.14], fourier [VI.25], Gauss, and others; semiclassical problems of polynomial and rational minimax approximation asso￾ciated with names such as chebyshev [VI.45] and Bern￾stein; and major newer topics, including splines, radial basis functions, and wavelets [VII.3]. We shall not have space to address these subjects, but in almost every area of numerical analysis it is a fact that, sooner or later, the discussion comes down to approximation theory. 3 Machine Arithmetic and Rounding Errors It is well-known that computers cannot represent real or complex numbers exactly. A quotient like 1/7 eval￾uated on a computer, for example, will normally yield an inexact result. (It would be different if we designed machines to work in base 7!) Computers approximate real numbers by a system of floating-point arithmetic, in which each number is represented in a digital equiv￾alent of scientific notation, so that the scale does not matter unless the number is so huge or tiny as to cause overflow or underflow. Floating-point arithmetic was invented by Konrad Zuse in Berlin in the 1930s, and

606 IV.Branches of Mathematics by the end of the 1950s it was standard across the the magnetic moment of the to the Bohr mag Until the 1980s.different computers had widely dif- to more than 12 or 13 digits of a .Thus IEEE ferent arithmetic properties.Then,in 1985.after y ears number on.th was adopted.or EEE arithmetic for short in two sense then,floating-point arithmetic is far his standard has su er to its ideal than is physics.It is a curiou mber coists of a 64-bit word divided widely regarded as an ugly and dangerous compromis or a signed expon ts the the 10 fathers of the field discovered that inexact arithmeti can be a source of danger causing errors in resu Comnuters do not merely renr ent numbers.of per orm operatio on t em such of rounding errors from microscopic to macroscopic d division,and e by certain mo se of these elementary opera s.In floating-point arith and Henric took great pains to publicize theisk ts of computed result mentary oper careless reliance on machine arithmetic.Thes is one of these four operations in its ideal formand to the current widespread impre is the same operation as realized on the computer sion that the main business of nu al analysis is ,assuming ith e quickly:rounding-e r analysis,while ofter x©y=(x*y)1+E discu y the cen ute valu ould remain. computer.n the IEEE system. 4 Numerical Linear Algebra 10 ing to compare the fineness of this discretization with ever since.There are several that of the discretizations of physics.nahandful is at the bottom of atoms or molecules in a line from o noint to anothe arrival of computers is on the order of 10(the cube root of Avogadro's num The starting point of this subject is nation.apro edure tha density.pressure.stress.strain.and temperature.Com perations.Equivalently,it solvesequations of the form puter arith etic,however,is more than a million time b,whe is an nxn r and x and ms the precision to which fundamental constantsar on computers around the world almost every time a known.such as(roughly)4 digits for the gravitationa ystem of lir ea equations is solved.Even if n is as 2008 The de

✐ 606 IV. Branches of Mathematics by the end of the 1950s it was standard across the computer industry. Until the 1980s, different computers had widely dif￾ferent arithmetic properties. Then, in 1985, after years of discussion, the IEEE (Institute of Electrical and Elec￾tronics Engineers) standard for binary floating-point arithmetic was adopted, or IEEE arithmetic for short. This standard has subsequently become nearly univer￾sal on processors of many kinds. An IEEE (double pre￾cision) real number consists of a 64-bit word divided into 53 bits for a signed fraction in base 2 and 11 bits for a signed exponent. Since 2−53 ≈ 1.1 × 10−16, IEEE numbers represent the numbers of the real line to a rel￾ative accuracy of about 16 digits. Since 2±210 ≈ 10±308, this system works for numbers up to about 10308 and down to about 10−308. Computers do not merely represent numbers, of course; they perform operations on them such as addition, subtraction, multiplication, and division, and more complicated results are obtained from sequences of these elementary operations. In floating-point arith￾metic, the computed result of each elementary opera￾tion is almost exactly correct in the following sense: if “∗” is one of these four operations in its ideal form and “-* ” is the same operation as realized on the computer, then for any floating-point numbers x and y, assuming that there is no underflow or overflow, x -* y = (x ∗ y)(1 + ε). Here ε is a very small quantity, no greater in abso￾lute value than a number known as machine epsilon, denoted by εmach, that measures the accuracy of the computer. In the IEEE system, εmach = 2−53 ≈ 1.1 × 10−16. Thus, on a computer, the interval [1, 2], for example, is approximated by about 1016 numbers. It is interest￾ing to compare the fineness of this discretization with that of the discretizations of physics. In a handful of solid or liquid or a balloonful of gas, the number of atoms or molecules in a line from one point to another is on the order of 108 (the cube root of Avogadro’s num￾ber). Such a system behaves enough like a continuum to justify our definitions of physical quantities such as density, pressure, stress, strain, and temperature. Com￾puter arithmetic, however, is more than a million times finer than this. Another comparison with physics con￾cerns the precision to which fundamental constants are known, such as (roughly) 4 digits for the gravitational constant G, 7 digits for Planck’s constant h and the ele￾mentary charge e, and 12 digits for the ratio µe/µB of the magnetic moment of the electron to the Bohr mag￾neton. At present, almost nothing in physics is known to more than 12 or 13 digits of accuracy. Thus IEEE numbers are orders of magnitude more precise than any number in science. (Of course, purely mathematical quantities like π are another matter.) In two senses, then, floating-point arithmetic is far closer to its ideal than is physics. It is a curious phenomenon that, nevertheless, it is floating-point arithmetic rather than the laws of physics that is widely regarded as an ugly and dangerous compromise. Numerical analysts themselves are partly to blame for this perception. In the 1950s and 1960s, the founding fathers of the field discovered that inexact arithmetic can be a source of danger, causing errors in results that “ought” to be right. The source of such prob￾lems is numerical instability: that is, the amplification of rounding errors from microscopic to macroscopic scale by certain modes of computation. These men, including von neumann [VI.91], Wilkinson, Forsythe, and Henrici, took great pains to publicize the risks of careless reliance on machine arithmetic. These risks are very real, but the message was communicated all too successfully, leading to the current widespread impres￾sion that the main business of numerical analysis is coping with rounding errors. In fact, the main busi￾ness of numerical analysis is designing algorithms that converge quickly; rounding-error analysis, while often a part of the discussion, is rarely the central issue. If rounding errors vanished, 90% of numerical analysis would remain. 4 Numerical Linear Algebra Linear algebra became a standard topic in undergradu￾ate mathematics curriculums in the 1950s and 1960s, and has remained there ever since. There are several reasons for this, but I think one is at the bottom of it: the importance of linear algebra has exploded since the arrival of computers. The starting point of this subject is Gaussian elimi￾nation, a procedure that can solve n linear equations in n unknowns using on the order of n3 arithmetic operations. Equivalently, it solves equations of the form Ax = b, where A is an n×n matrix and x and b are col￾umn vectors of size n. Gaussian elimination is invoked on computers around the world almost every time a system of linear equations is solved. Even if n is as large as 1000, the time required is well under a sec￾ond on a typical 2008 desktop machine. The idea of

IV.21.Numerical Analysis 607 elimination was first discovered by Chinese scholars with permuted rows.If the permutations are chosen to about 20 ago and more recent con argest entry b ow t in column The modern way of describing such algorithms how then has the additional property for ali ever was apparently introd as late as the 193 and J. racted from the ond row This ation can he interpreted as the multiplication of Aon the left by the practice,pivoting makes Gaussian elimination almost nsisting of the perfectly stable,and it is routinely done by almost all analogous row opera ions correspond to further multi- ations on th by lower-tri matrice Wilkinson and others that for certain exceptional matri on is still setting L=M-1, ancy represents an embarrassing gap at the heart of A-LU Here L isunit lower-triangular,that is,o wer-triangular d entr 04 ince o rep nd I which Gaussian eliminati n amplifies rounding errors can say that Gau in a certain sense exp tion is a process of lower-triangule the dimension.but a theorem to this effect has nev Nany of numerical lnear algebra been proved are also ba ased on writing a matrix as a product of matri Meanw e,beg nning in the lat 195( may s that this neld has a central the use of algorithms ha dogma: or UNITARY [Il.50s3]matrices,th at Is, real matrices algorithms-matrix factorizations Q: otesth orcom x ones w The st of such developments is the idea of QR factorization.I sanm matrx with mn a QR factorization a pro matrix A-QR eQ has orthonormal columns and R is upper-tri Soon after computers came into use it was obs that even for matrices that do have LU factoriza pression of the familar idea ot gram-schmidt ortho form i of o are ope amounts.Stability can be achieved by interchanging by elementary upp triangular matrices one could sav ord bring ma at the Gram-Schmid algorithm aims for Q and get: oin.Since pivoting acts on rows.it again corres A hig howed in 1958 that a dual strategy of orthogonalri- elimination with nivoting is ularization is more effective for mar PA-LU operations cach of which reflects across a ngular,I is unit lower-t hyperplane,one reduces A to upper-triangular form a orthogonal operations:one aims at R and gets Q

✐ IV.21. Numerical Analysis 607 elimination was first discovered by Chinese scholars about 2000 years ago, and more recent contributors include lagrange [VI.22], Gauss, and jacobi [VI.35]. The modern way of describing such algorithms, how￾ever, was apparently introduced as late as the 1930s. Suppose that, say, α times the first row of A is sub￾tracted from the second row. This operation can be interpreted as the multiplication of A on the left by the lower-triangular matrix M1 consisting of the identity with the additional nonzero entry m21 = −α. Further analogous row operations correspond to further multi￾plications on the left by lower-triangular matrices Mj . If k steps convert A to an upper-triangular matrix U, then we have MA = U with M = Mk ··· M2M1, or, upon setting L = M−1, A = LU. Here L is unit lower-triangular, that is, lower-triangular with all its diagonal entries equal to 1. Since U rep￾resents the target structure and L encodes the opera￾tions carried out to get there, we can say that Gaussian elimination is a process of lower-triangular upper￾triangularization. Many other algorithms of numerical linear algebra are also based on writing a matrix as a product of matri￾ces that have special properties. To borrow a phrase from biology, we may say that this field has a central dogma: algorithms ←→ matrix factorizations. In this framework we can quickly describe the next algorithm that needs to be considered. Not every matrix has an LU factorization; a 2 × 2 counterexample is the matrix A =  0 1 1 0 . Soon after computers came into use it was observed that even for matrices that do have LU factoriza￾tions, the pure form of Gaussian elimination is unsta￾ble, amplifying rounding errors by potentially large amounts. Stability can be achieved by interchanging rows during the elimination in order to bring maxi￾mal entries to the diagonal, a process known as piv￾oting. Since pivoting acts on rows, it again corresponds to a multiplication of A by other matrices on the left. The matrix factorization corresponding to Gaussian elimination with pivoting is PA = LU, where U is upper-triangular, L is unit lower-triangular, and P is a permutation matrix, i.e., the identity matrix with permuted rows. If the permutations are chosen to bring the largest entry below the diagonal in column k to the (k, k) position before the kth elimination step, then L has the additional property |ij |
1 for all i and j. The discovery of pivoting came quickly, but its theoretical analysis has proved astonishingly hard. In practice, pivoting makes Gaussian elimination almost perfectly stable, and it is routinely done by almost all computer programs that need to solve linear systems of equations. Yet it was realized in around 1960 by Wilkinson and others that for certain exceptional matri￾ces, Gaussian elimination is still unstable, even with pivoting. The lack of an explanation of this discrep￾ancy represents an embarrassing gap at the heart of numerical analysis. Experiments suggest that the frac￾tion of matrices (for example, among random matri￾ces with independent normally distributed entries) for which Gaussian elimination amplifies rounding errors by a factor greater than ρn1/2 is in a certain sense expo￾nentially small as a function of ρ as ρ → ∞, where n is the dimension, but a theorem to this effect has never been proved. Meanwhile, beginning in the late 1950s, the field of numerical linear algebra expanded in another direction: the use of algorithms based on orthogonal [III.50 §3] or unitary [III.50 §3] matrices, that is, real matrices with Q−1 = QT or complex ones with Q−1 = Q∗, where Q∗ denotes the conjugate transpose. The starting point of such developments is the idea of QR factorization. If A is an m × n matrix with m n, a QR factorization of A is a product A = QR, where Q has orthonormal columns and R is upper-tri￾angular. One can interpret this formula as a matrix ex￾pression of the familiar idea of Gram–Schmidt orthog￾onalization, in which the columns q1, q2,... of Q are determined one after another. These column opera￾tions correspond to multiplication of A on the right by elementary upper-triangular matrices. One could say that the Gram–Schmidt algorithm aims for Q and gets R as a by-product, and is thus a process of triangular orthogonalization. A big event was when Householder showed in 1958 that a dual strategy of orthogonal tri￾angularization is more effective for many purposes. In this approach, by applying a succession of elementary matrix operations each of which reflects Rm across a hyperplane, one reduces A to upper-triangular form via orthogonal operations: one aims at R and gets Q

608 IV.Branches of Mathematics ations preserve norms and thus do not amplify the for virtually every scentist with nobody but a few sp rounding errors introduced at each step. ialists knowing the details of how it is do ous re Istory is that QR factorization can be used by itself to solve least squares pro orthonormal b chmarks that all computer ma ufacturers run to In particular,one of the central problems of numerical is lucky ough to make the TOP500 list. dated twice linear algebra is the determinati ion or the a year since 1993.it is because of its prowess in solving vectors,then by forming a matrix x er trix prooml s Ax"b of dimensions ranging whose columns are thes e eigenvectors and a diagonal The eigenvalue decomposition is familiar to all math ematicians,but the development of numerical lin a AX -XD cene the and hence,since X is nonsingular SVD was discov AN [VI.52],anc eteenth c an A=XDX-1 beginning in around 1965.If A is an mx n matrix with m>n,an SVD of A is a factorization orthonormal eigenvectors always exists.giving A=USV* A=ODO where U ismx Visn×i where The 2 >o.One could compute the 1960s by Francis,Kublanovsk aya,and Wilkinson:the ating it to the eigenvalue problems for that doesnot square A.Computing the SVD is th stan ing the the in principle infinite.Nevertheless.its convergence is orm of the inverse lA-1l =la.in the case where A ctTraordnarvapidnthe condition number in the sense that at each step the number of cor K(A)=用A1MA-1‖=/On rect digits in one of the eigenvalue-eigenvector pairs t is also one of the great trlumphs of numerical analysis,and its impact through widely used least-squa res,computation of ranges and nullspaces n it l d in the and dete and Fortran and later to the software library gensystem )and its desce the dis on above cerns"clas ted in al. e n nsuing quarter-century brought in a whole new set IMSL and Numerical Recipes collecti of tools:methods for large-s LAB.Mro CI s MA se ne

✐ 608 IV. Branches of Mathematics as a by-product. The Householder method turns out to be more stable numerically, because orthogonal oper￾ations preserve norms and thus do not amplify the rounding errors introduced at each step. From the QR factorization sprang a rich collec￾tion of linear algebra algorithms in the 1960s. The QR factorization can be used by itself to solve least￾squares problems and construct orthonormal bases. More remarkable is its use as a step in other algorithms. In particular, one of the central problems of numerical linear algebra is the determination of the eigenvalues and eigenvectors of a square matrix A. If A has a com￾plete set of eigenvectors, then by forming a matrix X whose columns are these eigenvectors and a diagonal matrix D whose diagonal entries are the corresponding eigenvalues, we obtain AX = XD, and hence, since X is nonsingular, A = XDX−1, the eigenvalue decomposition. In the special case in which A is hermitian [III.50 §3], a complete set of orthonormal eigenvectors always exists, giving A = QDQ∗, where Q is unitary. The standard algorithm for com￾puting these factorizations was developed in the early 1960s by Francis, Kublanovskaya, and Wilkinson: the QR algorithm. Because polynomials of degree 5 or more cannot be solved by a formula, we know that eigen￾values cannot generally be computed in closed form. The QR algorithm is therefore necessarily an iterative one, involving a sequence of QR factorizations that is in principle infinite. Nevertheless, its convergence is extraordinarily rapid. In the symmetric case, for a typ￾ical matrix A, the QR algorithm converges cubically, in the sense that at each step the number of cor￾rect digits in one of the eigenvalue–eigenvector pairs approximately triples. The QR algorithm is one of the great triumphs of numerical analysis, and its impact through widely used software products has been enormous. Algorithms and analysis based on it led in the 1960s to computer codes in Algol and Fortran and later to the software library EISPACK (“Eigensystem Package”) and its descendant LAPACK. The same methods have also been incor￾porated in general-purpose numerical libraries such as the NAG, IMSL, and Numerical Recipes collections, and in problem-solving environments such as MAT￾LAB, Maple, and Mathematica. These developments have been so successful that the computation of matrix eigenvalues long ago became a “black box” operation for virtually every scientist, with nobody but a few spe￾cialists knowing the details of how it is done. A curi￾ous related story is that EISPACK’s relative LINPACK for solving linear systems of equations took on an unex￾pected function: it became the original basis for the benchmarks that all computer manufacturers run to test the speed of their computers. If a supercomputer is lucky enough to make the TOP500 list, updated twice a year since 1993, it is because of its prowess in solving certain matrix problems Ax = b of dimensions ranging from 100 into the millions. The eigenvalue decomposition is familiar to all math￾ematicians, but the development of numerical linear algebra has also brought its younger cousin onto the scene: the singular value decomposition (SVD). The SVD was discovered by Beltrami, jordan [VI.52], and sylvester [VI.42] in the late nineteenth century, and made famous by Golub and other numerical analysts beginning in around 1965. If A is an m× n matrix with m n, an SVD of A is a factorization A = UΣV ∗, where U is m×n with orthonormal columns, V is n×n and unitary, and Σ is diagonal with diagonal entries σ1 σ2 ··· σn 0. One could compute the SVD by relating it to the eigenvalue problems for AA∗ and A∗A, but this proves numerically unstable; a bet￾ter approach is to use a variant of the QR algorithm that does not square A. Computing the SVD is the stan￾dard route to determining the norm [III.62] A = σ1 (here · is the hilbert space [III.37] or “2” norm), the norm of the inverse A−1 = 1/σn in the case where A is square and nonsingular, or their product, known as the condition number, κ(A) = A A−1 = σ1/σn. It is also a step in an extraordinary variety of fur￾ther computational problems including rank-deficient least-squares, computation of ranges and nullspaces, determination of ranks, “total least-squares,” low-rank approximation, and determination of angles between subspaces. All the discussion above concerns “classical” numer￾ical linear algebra, born in the period 1950–75. The ensuing quarter-century brought in a whole new set of tools: methods for large-scale problems based on Krylov subspace iterations. The idea of these iterations is as follows. Suppose a linear algebra problem is given