604 IV.Branches of Mathematics really goes on is a far more interesting process of be done even in principle by formulas.for most mathe. matical problems cannot be solved by a finite sequence Press/Elsevier. of elementary operations.What happens instead is Ke.Shen,and M.Vyaly 2002.Classical and Quan rision,or a hundred.For a scientific or engineering um Computation.Providence,RI:American Mathematical application,such an answer may be as good as exact. Kushilevitz,E,and N.Nisan.1996.Communication Com- We can illustra a tutorial).In Handbook on p(2)-C0+C1z+C222+C323+C42 s.Bulletin of the E an Association and another of degree 5, ions of extractors q(z)-do+diz+d2z2+d3z+dz+dsz It is well-known that there is an explicit formula that Algebralc complexity theory.In Hand- expresses the roots of p in terms of radicals (discov book of Theoretical C red by Ferrari around bridge.MA-MIT Press/Elsevier. more than 250 vears later see THe O THE QUINTIC [V.21]for more details).Thus,in a cer IV.21 Numerical Analysis tain philosophic se the root-findir Lloyd N.Trefethen know the roots of one of these polynomials,he or she 1 The Need for Numerical Computation vill turn to a computer and get an answer to sixt Everyone knows that when scientists and engineers of precs ha a need numerical answers to mathematical problems. In thend.Did they turn to computers there is a wide the answer is certainly no,but what about p?Maybe. naybe not.Most of the time,the user neither knows The power of numbers has been extraordinary.It is proba y not one ofe noted that the scientific revolution was setn made it a principle from memory en G others Here are three more examples of problems that can be solved in principle by a finite sequence of elementary and,in the remarkable cycle whose fruits are all around operations,like fnding the roots of p. us,finer measurements led to refine d laws,which ir yand st ()Linear equations:solve a system of n linear equ in the physical sciences could be achieved,or a signifi- (Linear programming:minimize a linear function cant engineering product developed,without numerical of n variables subject to m linear constraints. )Traveling sale man problem:find the shortest tou ertainly play etween i citle nc d he Many people imagine that scientists and mathemati gen Canhtkctmitheosof4.cano formulas,and th h,b essary results.The reality is nothing like this.What
✐ 604 IV. Branches of Mathematics Karp, R. M., and V. Ramachandran. 1990. Parallel algorithms for shared-memory machines. In Handbook of Theoretical Computer Science, volume A, Algorithms and Complexity, 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 Quantum Computation. Providence, RI: American Mathematical Society. Kushilevitz, E., and N. Nisan. 1996. Communication Complexity. 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 constructions 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 Handbook of Theoretical Computer Science, volume A, Algorithms and Complexity, edited by J. van Leeuwen. Cambridge, 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 widespread 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 measurements 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 measurements. 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 mathematicians generate formulas, and then, by inserting numbers into these formulas, computers grind out the necessary results. The reality is nothing like this. What really goes on is a far more interesting process of execution of algorithms. In most cases the job could not be done even in principle by formulas, for most mathematical 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 precision, 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. Suppose 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 (discovered 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 certain 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 equations 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
IV21.Numerical Analysis 605 5 r 1 ation (ODE) ics began but now mainly in the hands of (v)Solve a partial differential equation (PDE). Can we conclude that (i)-(iii)will be easier than (iv)-(viii ly not.Pro em (is usually ver thad ma Problems (vi)and (vii)are usually rather easy.at leas century from the 1950s,machines sped up by a factor of around 10 if the dimension.Problems (and (v) but so did the best algorithms known is small like 100 and often very hard when n is larg like 1000000.In fact,in the matters philosophy is Half a century on.umericl analysis has grownn the largest branches gnores the exact solution and uses approximate(but dozens of mathematical journals as well as applica- ons journals across the sciences and ast!)meth is the study of algorithms for solving the problems of continuous mathematics.by ers,we have reached a point wh which we e most of the clas cal p ear progra ences ming and the traveling salesman problem the algorithms that make this possible were invented ed o wer the real numbers,but not their dis rete ana since 1950 0 possible future trends This field encompass classical questions of inter 2 A Brief History po ANALY Throughout history,leading mathematicians have been [VL251.Gauss.and othe semiclassical problems of to the di c applic ions,and i amln山 thm still in use todav.GAuss IVL261. as usual,is an out stein and major newer topics. snlines radial standing example Among n mnany othe basis functions.and WAVELETS IVII.31 We shall t have sp 1795)systems of lin ar equations (809)and nume (0 dre o than almos ical quadrature (1814), well as inventing THE FAS or later.the discussion comes down to approximation theory. Cooley and Tukey in 1965. 3 Machine Arithmetic and Rounding Errors the numerical sid of mather d t It is well-known that computers cannot represent real or complex numbers exactly.A quoti ent like 1/7 eval. the grow s genera and of grea ated on a cor er,fo orm matical rigor had to be the heart of the matter.For chines to work in base 7!) Computers appr imate xample,many advances of the entieth cer bers by a system of float s about infinity,a subject relativel y far from alent of scientific notation.so that the scale docs no merical calculation. nd in the 1940s nted by ad the 1930s
✐ 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 linear programming and the traveling salesman problem posed over the real numbers, but not their discrete analogues.) 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 outstanding example. Among many other contributions, he made crucial advances in least-squares data fitting (1795), systems of linear equations (1809), and numerical 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, mathematical rigor had to be the heart of the matter. For example, many advances of the early twentieth century 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 mathematics began to explode, but now mainly in the hands of specialists. New journals were founded such as Mathematics of Computation (1943) and Numerische Mathematik (1959). The revolution was sparked by hardware, but it included mathematical and algorithmic developments that had nothing to do with hardware. In the halfcentury 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 into one of the largest branches of mathematics, the specialty of thousands of researchers who publish in dozens of mathematical journals as well as applications journals across the sciences and engineering. Thanks to the efforts of these people going back many decades, and thanks to ever more powerful computers, we have reached a point where most of the classical 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 interpolation, series expansions, and harmonic analysis [IV.11] associated with newton [VI.14], fourier [VI.25], Gauss, and others; semiclassical problems of polynomial and rational minimax approximation associated with names such as chebyshev [VI.45] and Bernstein; 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 evaluated 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 equivalent 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 electron to the Bohr mag in rs had widely dif neton.At pres st nothing in physics ise thar of discussion,the IEEE(nstitute of Electrical and Ele any number in science.(Of course,purely mathematical ues lke n are an rithmetic is fa This standard has subsequently become nearly univer closer to its ideal than is physics.It is a curious al on pro ors of many kir ·An IEEE(do phenomenon that. floating-poin into bits for a signed fraction in base2 and 11 bits widely regarded as an ugly and dangerous compromise for a signed exponent.Since ≈1.1×10-16，E Numerical analysts themselves are partly to blame for Geld dis this system work an be a source of danger causing errors in results down to abou that ought to e source of course thev of rounding errors from microscopic to ma addition,sub action,multipli ation,and division,and These me s are o NN IV int arith metic,the computedr ault of e careless rellance on machine arithmetic.These risks are ℃ ery real.but the message was communicated all to is the same operation as realized on the computer ion that the main business of numerical analysis is then for any floating-point numbe oping with rounding errors In fact,the main b hat there is no underflo ess of n sisis desig X©v=(X*Y)(1+E) Here s is a very small quantity,no greater in abso arely the central issue.If 90%of numerical analysis lute value than a number known as machine epsilon denoted by 10gutc .In the 4 Numerical Linear Algebra Thus,on a computer,thei erval [1,2],for example the of th fphysics.In a handful of sons for this,but I think one is at the bottom of i d or of gas,th number ance or lnear algebra has exploded since the The starting point of this subject is Gaussian elimi ber).Such a system ehaves enough like a continuum nation,a procedure that can solve n linear equations ord 0 onsof the for puter arithmetic,however,is more than a million times xb,where Ais matri and x and b are co than thi Another nparison with physics co imn vectors of siz nd the known.such as (roughly)4 digits for the gravitationa ystem of linear equations is solved.Even if n is as constant G.7digits for Planck's constant h and the el mentary charge e,and 12 digits for the ratio ue/uB o on a typical 2008
✐ 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 different arithmetic properties. Then, in 1985, after years of discussion, the IEEE (Institute of Electrical and Electronics Engineers) standard for binary floating-point arithmetic was adopted, or IEEE arithmetic for short. This standard has subsequently become nearly universal on processors of many kinds. An IEEE (double precision) 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 relative 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 arithmetic, the computed result of each elementary operation 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 absolute 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 interesting 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 number). Such a system behaves enough like a continuum to justify our definitions of physical quantities such as density, pressure, stress, strain, and temperature. Computer arithmetic, however, is more than a million times finer than this. Another comparison with physics concerns 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 elementary charge e, and 12 digits for the ratio µe/µB of the magnetic moment of the electron to the Bohr magneton. 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 problems 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 impression that the main business of numerical analysis is coping with rounding errors. In fact, the main business 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 undergraduate 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 elimination, 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 column 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 second on a typical 2008 desktop machine. The idea of
IV21.Numerical Analysis 607 If the are chosen to include [V1221 Gauss and Iacori [VL351 to the (kk) osition before the kth elimination step The modern way of describing such algorith ms,how then L has the additional property Ifiil s 1 for all i and ts cted from the second row.This operation can be ly hard.In s the multip left by pe with the additional nonzero entry m Furthe ed to solve correspond to further mult If k ster rt A to per-trianaular matrix son and ot then we have.or.upon pivoting.The lack of an explanation of this discrep setting L M ancy represents an embarrassing gap at the heart of A-LU ents sugges Here ngular,tha 1.er-tri the ces with independent norr ally distributed entries)for ents the target st ture and i encodes the opera ch Gaussian e amplines rounding error inCamiedouttogetth 1 ecan s say tha uppe n is the dimension,but a theorem to this effect has never Many other algorithms of numerical linear algebra n prove in the late 1950s that hat edon writing a matrix as a product of matr numerical linear algebrae xpanded in another direction from biology,we may say that this field has a central he use of algorithms ba d on ORTHOGONAL.50s3 dogma: 0s3)m algorithms-matrix factorizations notes the conjugate transpose.The starting point n this framework we ribe th of A is a product A-OR A- al colu int matrice idea of Gram-Schmidt orthog nat even I actor a th 41.q or Q ar rounding erre b potentially large correspond to multiplication of A on the right eve by ele mentary upper-t riangula r matrices.One could say mal entries to the dia nal.a process known as piv oting.Since pivoting acts on rows,it again corresponds lizatio A big event was when Householder facton 958 that a d elimination with pivoting is this approach.by applying a succession of elementary PA=LU matrix upper-a operations each of which reflects R gular,L ng one redu uppe is a perm
✐ 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, however, was apparently introduced as late as the 1930s. Suppose that, say, α times the first row of A is subtracted 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 multiplications 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 represents the target structure and L encodes the operations carried out to get there, we can say that Gaussian elimination is a process of lower-triangular uppertriangularization. Many other algorithms of numerical linear algebra are also based on writing a matrix as a product of matrices 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 factorizations, the pure form of Gaussian elimination is unstable, amplifying rounding errors by potentially large amounts. Stability can be achieved by interchanging rows during the elimination in order to bring maximal entries to the diagonal, a process known as pivoting. 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 matrices, Gaussian elimination is still unstable, even with pivoting. The lack of an explanation of this discrepancy represents an embarrassing gap at the heart of numerical analysis. Experiments suggest that the fraction of matrices (for example, among random matrices with independent normally distributed entries) for which Gaussian elimination amplifies rounding errors by a factor greater than ρn1/2 is in a certain sense exponentially 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-triangular. One can interpret this formula as a matrix expression of the familiar idea of Gram–Schmidt orthogonalization, in which the columns q1, q2,... of Q are determined one after another. These column operations 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 triangularization 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 havebeen so successful that the computation of mati num orthogon s long ago be peratio rounding errors introduced at each step. cialists kn wing the details of how it is done.A curi From the O factoriz ation sprang a rich ous related story is that EISPACK's relative LINPACK for OR factorization can be ems and construct orthonormal bases benchmarks that all computer manufacturers run to se as a st the spee the computers.a super linear algebra is the determination of the a vear since 1993.it is because of its prowess in solving and eigen ctors of squa ertain matrix probler ms Ax "b of dimensions ranging ors and a diap ematicians,but the development of numerical linear algebra has rought it ont the AX XD. SVD was dise and hence,since X is nonsingular A-XDX-1, SYLVESTER [VL s by m n th specia m≥，an SVD of A is a factorization orthonormal eigenvectors alway s exists,giving A=UEV* A=QDQ*, thonormal columns,Visnxn algo for com and AA,b values cannot generally be computed in closed form that does not square A.Computing the SVD is the stan. The QR algorithm is therefore necessarily an iterative ermining the [11.62]IAll the HI A-I orm)the is square and nonsingular,or their product,known as .the the condition number. rect digits in one of the eigenvalue K(A)=AllA-=/om approx triples t isalso a step numerical analys and its impact h east-squares,computation of ranges and nullspaces software products has bee letermination of s,to squa low-ran on it le p EISPACK ystem Package")and its descendant All the discu ssion above co erms“lass Canm LAPA The sar method en inc cal linear algebra,borr 5.Th of tools:methods for large- cale problems based on subspace iterations.The idea of these iterations Maple,and s as foll ows.Suppose a linear algebra problem is giver
✐ 608 IV. Branches of Mathematics as a by-product. The Householder method turns out to be more stable numerically, because orthogonal operations preserve norms and thus do not amplify the rounding errors introduced at each step. From the QR factorization sprang a rich collection of linear algebra algorithms in the 1960s. The QR factorization can be used by itself to solve leastsquares 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 complete 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 computing 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 eigenvalues 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 typical matrix A, the QR algorithm converges cubically, in the sense that at each step the number of correct 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 incorporated in general-purpose numerical libraries such as the NAG, IMSL, and Numerical Recipes collections, and in problem-solving environments such as MATLAB, 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 specialists knowing the details of how it is done. A curious related story is that EISPACK’s relative LINPACK for solving linear systems of equations took on an unexpected 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 mathematicians, 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 better approach is to use a variant of the QR algorithm that does not square A. Computing the SVD is the standard 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 further 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” numerical 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
IV21.Numerical Analysis 609 that invoves a matrix of large dimension,say n xthat satisfies a certain variational propert 1990.Is there a"fast matrix inverse"in store for us? as minimizing xTAx -x'b (for solving Ax =b if Ais Numerical Solution of Differential Equations ric).Now if K is a k-dimensional subspace of with Long before much attention was pald to linear alge n,th n itmay be po var The magical choice of K&is a Krylov subspace cal in pro The ppack to NewToN,and even to ARCHIMEDes [V131.Th Kx(A.q)=span(q.Aq. 4k-1a sic qu las are d he ide for an intial q.Fo re ons th t have tas in these subspaces often converge very rapidly to the Newto R”aski values o diverge at a rate as high as2 as the gund If the points are chosen optimaly,then to ten-digit precision in just a few hundred iterations algorithm nd is stable.It tur that t mal points are roots of Legendre polyn omials rylo subspace iteratic s originated with th ed near th sketc e gr ient and Lanczos iterations pub poses is Clenshaw-Curtis quadrature,where the inte olve problems of a large enough s ale for the meth ation point ne cos(jm/n). ≤n.Th ds to b ompetitiv They took off i the 1970s with and unlike Gau quadrature can he ted in (nlog n)ope rations by the fast Fourier transform ditioning.In precondit oning a system Ax b,one eplaces it by a mathemat cally equivalent system such of potential theory MAX Mb 85 of ngular matrix M.If M is well ch as are based on polynomial interpolation in equally and a Krylov subspace iteration a p ally number Since the 1970s.preconditioned matrix iterations multistep methods for the numerical solution of ODEs. The idea here is that for an initial value prob em u' ave emerged as an indispensable tool of computa time most heavily cited article in all of mathematicsn the n=n,n≥0. s the 198 aper by We then replac us to nonsymmetric matrices approximate values Finally,we must mentio n the bigg gest unsolved pro v"u(tn). n≥0 he inverted in oina)one rations for ,2 (The problems of solving system Ax b or com EULER IVL19L is vale v+1-vn+△tftn,v
✐ IV.21. Numerical Analysis 609 that involves a matrix of large dimension, say n . 1000. The solution may be characterized as the vector x ∈ Rn that satisfies a certain variational property such as minimizing 1 2 xTAx −xTb (for solving Ax = b if A is symmetric positive definite) or being a stationary point of (xTAx)/(xTx) (for solving Ax = λx if A is symmetric). Now if Kk is a k-dimensional subspace of Rn with k / n, then it may be possible to solve the same variational problem much more quickly in that subspace. The magical choice of Kk is a Krylov subspace Kk(A, q) = span(q, Aq, . . . , Ak−1q) for an initial vector q. For reasons that have fascinating connections with approximation theory, solutions in these subspaces often converge very rapidly to the exact solution in Rn as k increases, if the eigenvalues of A are favorably distributed. For example, it is often possible to solve a matrix problem involving 105 unknowns to ten-digit precision in just a few hundred iterations. The speedup compared with the classical algorithms may be a factor of thousands. Krylov subspace iterations originated with the conjugate gradient and Lanczos iterations published in 1952, but in those years computers were not powerful enough to solve problems of a large enough scale for the methods to be competitive. They took off in the 1970s with the work of Reid and Paige and especially van der Vorst and Meijerink, who made famous the idea of preconditioning. In preconditioning a system Ax = b, one replaces it by a mathematically equivalent system such as MAx = Mb for some nonsingular matrix M. If M is well chosen, the new problem involving MA may have favorably distributed eigenvalues and a Krylov subspace iteration may solve it quickly. Since the 1970s, preconditioned matrix iterations have emerged as an indispensable tool of computational science. As one indication of their prominence we may note that in 2001, Thomson ISI announced that the most heavily cited article in all of mathematics in the 1990s was the 1989 paper by van der Vorst introducing Bi-CGStab, a generalization of conjugate gradients for nonsymmetric matrices. Finally, we must mention the biggest unsolved problem in numerical analysis. Can an arbitrary n×n matrix A be inverted in O(nα) operations for every α > 2? (The problems of solving a system Ax = b or computing a matrix product AB are equivalent.) Gaussian elimination has α = 3, and the exponent shrinks as far as 2.376 for certain recursive (though impractical) algorithms published by Coppersmith and Winograd in 1990. Is there a “fast matrix inverse” in store for us? 5 Numerical Solution of Differential Equations Long before much attention was paid to linear algebra, mathematicians developed numerical methods to solve problems of analysis. The problem of numerical integration or quadrature goes back to Gauss and newton [VI.14], and even to archimedes [VI.3]. The classic quadrature formulas are derived from the idea of interpolating data at n + 1 points by a polynomial of degree n, then integrating the polynomial exactly. Equally spaced interpolation points give the Newton– Cotes formulas, which are useful for small degrees but diverge at a rate as high as 2n as n → ∞: the Runge phenomenon. If the points are chosen optimally, then the result is Gauss quadrature, which converges rapidly and is numerically stable. It turns out that these optimal points are roots of Legendre polynomials, which are clustered near the endpoints. (A proof is sketched in special functions [III.85].) Equally good for most purposes is Clenshaw–Curtis quadrature, where the interpolation points become cos(jπ /n), 0 j n. This quadrature method is also stable and rapidly convergent, and unlike Gauss quadrature can be executed in O(n log n) operations by the fast Fourier transform. The explanation of why clustered points are necessary for effective quadrature rules is related to the subject of potential theory. Around 1850 another problem of analysis began to get attention: the solution of ODEs. The Adams formulas are based on polynomial interpolation in equally spaced points, which in practice typically number fewer than ten. These were the first of what are now called multistep methods for the numerical solution of ODEs. The idea here is that for an initial value problem u = f (t, u) with independent variable t > 0, we pick a small time step ∆t > 0 and consider a finite set of time values tn = n∆t, n 0. We then replace the ODE by an algebraic approximation that enables us to calculate a succession of approximate values vn ≈ u(tn), n 0. (The superscript here is just a superscript, not a power.) The simplest such approximate formula, going back to euler [VI.19], is vn+1 = vn + ∆tf (tn, vn),
610 IV.Branches of Mathematics or,using the abbreviation fmf(tn.). ship between the number of"stages"(- vn+1=yn+△tfn per step)and th dof accuracy Both the ODE sel f and its num Kutta in 1901 and have order but it was not unti ut.and y become vectors of an appropriate 1963 that it was proved that s 6 stages are required o achieve ord The analysis of su dimension.The Adams formulas are higher-order gen s or that other areas.and a key fgure in this area since the 1960s fourth-order Adams-Bashforth formula is vn+1=v+△(55m-59P-1+37m-2-9fm-3) es are orders are rarely needed for practical purposes. reflects a new element in When computers began to be used to solve differ ence as ial equat at the stahility as hefore this phrase refers will no orde o the unb unded amplification of local errors by a enabling excellent accuracy for all kinds of computa ally those e ant ounding.Instability typically manifests itself as an lator r in the computed solution that b Most unfortunately,the habi in the numerical a s1 rature I em speak no ofe of the mund Dahlqust Dahlquist saw that the phenom more precisely their analyzed with great power and g 10 as distinct f rom rounding error.This ubiquitous lan e of the t r ing the hirth of mod n tone,bu ern num ical analysis. This land eoren f Kutta.For example.here are the forn consistency+stability=convergence nulas of the famous fourth-order Runge-Kutta method,which advance a perty that the discrete formula has locally positive function f: a=△tf(t,v, is the proper e an b=△tftn+△t,v+a). ence is the property that as t 0,in the absence c=△tf(tn+号△t,vn+号b). of rounding error the num erical solution co d=△tf(tn+△t.vn+cl, n+1=vn+言(a+2b+2c+d and haps in the air in the sense that practitioners realized Runge-Kutta methods tend to be easier to implement tha schem was not unstab the der to analy. nswer His theor derive the of dams-Bashforth a wide class of numerical methods as co mputer methods for ODEs were being devel oped, the same was happening for the much bigge
✐ 610 IV. Branches of Mathematics or, using the abbreviation f n = f (tn, vn), vn+1 = vn + ∆tf n. Both the ODE itself and its numerical approximation may involve one equation or many, in which case u(t, x) and vn become vectors of an appropriate dimension. The Adams formulas are higher-order generalizations of Euler’s formula that are much more effi- cient at generating accurate solutions. For example, the fourth-order Adams–Bashforth formula is vn+1 = vn + 1 24∆t(55f n − 59f n−1 + 37f n−2 − 9f n−3). The term “fourth-order” reflects a new element in the numerical treatment of problems of analysis: the appearance of questions of convergence as ∆t → 0. The formula above is of fourth order in the sense that it will normally converge at the rate O((∆t)4). The orders employed in practice are most often in the range 3–6, enabling excellent accuracy for all kinds of computations, typically in the range of 3–10 digits, and higherorder formulas are occasionally used when still more accuracy is needed. Most unfortunately, the habit in the numerical analysis literature is to speak not of the convergence of these magnificently efficient methods, but of their error, or more precisely their discretization or truncation error as distinct from rounding error. This ubiquitous language of error analysis is dismal in tone, but seems ineradicable. At the turn of the twentieth century, the second great class of ODE algorithms, known as Runge–Kutta or one-step methods, was developed by Runge, Heun, and Kutta. For example, here are the formulas of the famous fourth-order Runge–Kutta method, which advance a numerical solution (again scalar or system) from time step tn to tn+1 with the aid of four evaluations of the function f : a = ∆tf (tn, vn), b = ∆tf (tn + 1 2∆t,vn + 1 2 a), c = ∆tf (tn + 1 2∆t,vn + 1 2b), d = ∆tf (tn + ∆t,vn + c), vn+1 = vn + 1 6 (a + 2b + 2c + d). Runge–Kutta methods tend to be easier to implement but sometimes harder to analyze than multistep formulas. For example, for any s, it is a trivial matter to derive the coefficients of the s-step Adams–Bashforth formula, which has order of accuracy p = s. For Runge– Kutta methods, by contrast, there is no simple relationship between the number of “stages” (i.e., function evaluations per step) and the attainable order of accuracy. The classical methods with s = 1, 2, 3, 4 were known to Kutta in 1901 and have order p = s, but it was not until 1963 that it was proved that s = 6 stages are required to achieve order p = 5. The analysis of such problems involves beautiful mathematics from graph theory and other areas, and a key figure in this area since the 1960s has been John Butcher. For orders p = 6, 7, 8 the minimal numbers of stages are s = 7, 9, 11, while for p > 8 exact minima are not known. Fortunately, these higher orders are rarely needed for practical purposes. When computers began to be used to solve differential equations after World War II, a phenomenon of the greatest practical importance appeared: once again, numerical instability. As before, this phrase refers to the unbounded amplification of local errors by a computational process, but now the dominant local errors are usually those of discretization rather than rounding. Instability typically manifests itself as an oscillatory error in the computed solution that blows up exponentially as more numerical steps are taken. One mathematician concerned with this effect was Germund Dahlquist. Dahlquist saw that the phenomenon could be analyzed with great power and generality, and some people regard the appearance of his 1956 paper as one of the events marking the birth of modern numerical analysis. This landmark paper introduced what might be called the fundamental theorem of numerical analysis: consistency + stability = convergence. The theory is based on precise definitions of these three notions along the following lines. Consistency is the property that the discrete formula has locally positive order of accuracy and thus models the right ODE. Stability is the property that errors introduced at one time step cannot grow unboundedly at later times. Convergence is the property that as ∆t → 0, in the absence of rounding errors, the numerical solution converges to the correct result. Before Dahlquist’s paper, the idea of an equivalence of stability and convergence was perhaps in the air in the sense that practitioners realized that if a numerical scheme was not unstable, then it would probably give a good approximation to the right answer. His theory gave rigorous form to that idea for a wide class of numerical methods. As computer methods for ODEs were being developed, the same was happening for the much bigger
IV.21.Numerical Analysis 611 In the haif-century since von Neumann died,the for applications in st and m into a breathtakingly powerful subject known asom rology,and further developed by Southwell:in 1928 ational fluid dynamics.Early treatments of linear er on h enc But although the Courant-Friedrichs-Lewy work later is now a routine matter to solve problems involving ecame Tamous,the impact of thes e idea before con th hur t developed ionsare linear or the grids areniform the early years was the group of researchers around at the Los Alamos laboratory.including the st as for ODEs on Neun and hi discovered that some numerical methods for PDEswere use wino For examp t they Many of these succ esses have been facilitated by pick space and time steps Ax and At for a regular anoth ca for solving PDEs tha xi=j△x.tn=n△t,jin≥0 imating a differential operator by a difference quo nd replace the PDE by algebraic fo pute a succession of approximate values: mulas that com tself by functi mef tha =u(,x..n≥0 For instance,one migh partition the suc A well-known discretization for this purpose is the Lax- insist Wendroff formula: each piece is a polynomial of small degree The solu u*1=+A(1-1)+A2(1-2+%, ed by solving a al form of t where A=△t/A which can be generalized to nor and th re is often a aua ntee that the camnuted opti mal within that subspace.Fin emen uncti 06 ods are known for theirf andling compl Von Neumann and others realized that the presence or ted gec y are entir absence of such instabil ivil engineering.The number of books and a es that linea have been published about finite-element methods is in sis."Experience indicated that,as a practical matter,a mature field of numerical salution of A theory PDEs,what aspect of the current state of the art would se Rich son or Courant, myer in 1956,the same year as Dahlquist's paper.Many algorithms of linear algebra.The solution of a large letails were di as re scale PDE pro applied to nonlinear ones but broadly speaking the rix iteration result followed the same pattern of quating co that utilizes a finite-difference preconditioner imple ency athen nla
✐ IV.21. Numerical Analysis 611 subject of PDEs. Discrete numerical methods for solving PDEs had been invented around 1910 by Richardson for applications in stress analysis and meteorology, and further developed by Southwell; in 1928 there was also a theoretical paper on finite-difference methods by courant [VI.83], Friedrichs, and Lewy. But although the Courant–Friedrichs–Lewy work later became famous, the impact of these ideas before computers came along was limited. After that point the subject developed quickly. Particularly influential in the early years was the group of researchers around von Neumann at the Los Alamos laboratory, including the young Peter Lax. Just as for ODEs, von Neumann and his colleagues discovered that some numerical methods for PDEs were subject to catastrophic instabilities. For example, to solve the linear wave equation ut = ux numerically we may pick space and time steps ∆x and ∆t for a regular grid, xj = j∆x, tn = n∆t, j,n 0, and replace the PDE by algebraic formulas that compute a succession of approximate values: vn j ≈ u(tn, xj ), j, n 0. A well-known discretization for this purpose is the Lax– Wendroff formula: vn+1 j = vn j + 1 2λ(vn j+1 −vn j−1)+ 1 2λ2(vn j+1 −2vn j +vn j−1), where λ = ∆t/∆x, which can be generalized to nonlinear systems of hyperbolic conservation laws in one dimension. For ut = ux, if λ is held fixed at a value less than or equal to 1, the method will converge to the correct solution as ∆x, ∆t → 0 (ignoring rounding errors). If λ is greater than 1, on the other hand, it will explode. Von Neumann and others realized that the presence or absence of such instabilities could be tested, at least for linear constant-coefficient problems, by discrete fourier analysis [III.27] in x: “von Neumann analysis.” Experience indicated that, as a practical matter, a method would succeed if it was not unstable. A theory soon appeared that gave rigor to this observation: the Lax equivalence theorem, published by Lax and Richtmyer in 1956, the same year as Dahlquist’s paper. Many details were different—this theory was restricted to linear equations whereas Dahlquist’s theory for ODEs also applied to nonlinear ones—but broadly speaking the new result followed the same pattern of equating convergence to consistency plus stability. Mathematically, the key point was the uniform boundedness principle. In the half-century since von Neumann died, the Lax–Wendroff formula and its relatives have grown into a breathtakingly powerful subject known as computational fluid dynamics. Early treatments of linear and nonlinear equations in one space dimension soon moved to two dimensions and eventually to three. It is now a routine matter to solve problems involving millions of variables on computational grids with hundreds of points in each of three directions. The equations are linear or nonlinear; the grids are uniform or nonuniform, often adaptively refined to give special attention to boundary layers and other fast-changing features; the applications are everywhere. Numerical methods were used first to model airfoils, then whole wings, then whole aircraft. Engineers still use wind tunnels, but they rely more on computations. Many of these successes have been facilitated by another numerical technology for solving PDEs that emerged in the 1960s from diverse roots in engineering and mathematics: finite elements. Instead of approximating a differential operator by a difference quotient, finite-element methods approximate the solution itself by functions f that can be broken up into simple pieces. For instance, one might partition the domain of f into elementary sets such as triangles or tetrahedra and insist that the restriction of f to each piece is a polynomial of small degree. The solution is obtained by solving a variational form of the PDE within the corresponding finite-dimensional subspace, and there is often a guarantee that the computed solution is optimal within that subspace. Finite-element methods have taken advantage of tools of functional analysis to develop to a very mature state. These methods are known for their flexibility in handling complicated geometries, and in particular they are entirely dominant in applications in structural mechanics and civil engineering. The number of books and articles that have been published about finite-element methods is in excess of 10 000. In the vast and mature field of numerical solution of PDEs, what aspect of the current state of the art would most surprise Richardson or Courant, Friedrichs, and Lewy? I think it is the universal dependence on exotic algorithms of linear algebra. The solution of a largescale PDE problem in three dimensions may require a system of a million equations to be solved at each time step. This may be achieved by a GMRES matrix iteration that utilizes a finite-difference preconditioner implemented by a Bi-CGStab iteration relying on another multigrid preconditioner. Such stacking of tools was
612 IV.Branches of Mathematics surely not imagined by the early computer pioneers. their Taylor series;indeed,Arnol'd has argued that Tay to numertn 1947 vton's aa math the ucial tool for combating instability is the use o veryone knows the idea of Newton's method at the kth stp. imate x derive a required. x(k+1)=x(-F(x())/F(x()) Here are some examples that illustrate the suc. today's SCHRODINGER EQUATION:structural mechan other functions F and to systems of two equations cda languag NS I at a point c matr [IV.13$1.1D);cosmology(the Einstein equationsk oil dis diffusion equationsk:tsu ami modeling (the shallo water equa optical fib x)=x()-(J(x())-IF(x(k) Malik equation):metall rgy ithe Cahn-Hilliard eaua practice that to get .9s2D Jx®)xk+-x)=-Fx) 6 Numerical Optimization As long as is Lipschi convergence of this iteration is quadratic: The xk+1》-x41=O("x)-x,2) several variables and the closely related problem of Students en think it might b solution of nonlinear systems of equations.The devel formulas to enhance the e onent in this estimate to 3 Ho this is an llusio Taking two steps a in effi forward in part by a community of scholars with close links to operations research and economics. iency betweer n quadratic and quartic is at best acon- Calculus students lea n that a smooth fun on ma a houndary.The same two no ssibilities characterize the hich two big strands of the fi d of optimization.At on of pro rom the point related to multivariate calculus.At the other are prob of multivariate calculus,it is ems of line programming where the imizing a scalar function f of a variable xE":to ind a d minimum, a zero of th h Unconstrained nonlinear optimization is an old sub rix known as the es of f with o 22f ons by tne H6)=ax,1n
✐ 612 IV. Branches of Mathematics surely not imagined by the early computer pioneers. The need for it ultimately traces to numerical instability, for as Crank and Nicolson first noted in 1947, the crucial tool for combating instability is the use of implicit formulas, which couple together unknowns at the new time step tn+1, and it is in implementing this coupling that solutions of systems of equations are required. Here are some examples that illustrate the successful reliance of today’s science and engineering on the numerical solution of PDEs: chemistry (the schrödinger equation [III.83]); structural mechanics (the equations of elasticity); weather prediction (the geostrophic equations); turbine design (the navier– stokes equations [III.23]); acoustics (the Helmholtz equation); telecommunications (maxwell’s equations [IV.13 §1.1]); cosmology (the Einstein equations); oil discovery (the migration equations); groundwater remediation (Darcy’s law); integrated circuit design (the drift diffusion equations); tsunami modeling (the shallowwater equations); optical fibers (the nonlinear wave equations [III.49]); image enhancement (the Perona– Malik equation); metallurgy (the Cahn–Hilliard equation); pricing financial options (the black–scholes equation [VII.9 §2]). 6 Numerical Optimization The third great branch of numerical analysis is optimization, that is, the minimization of functions of several variables and the closely related problem of solution of nonlinear systems of equations. The development of optimization has been somewhat independent of that of the rest of numerical analysis, carried forward in part by a community of scholars with close links to operations research and economics. Calculus students learn that a smooth function may achieve an extremum at a point of zero derivative, or at a boundary. The same two possibilities characterize the two big strands of the field of optimization. At one end there are problems of finding interior zeros and minima of unconstrained nonlinear functions by methods related to multivariate calculus. At the other are problems of linear programming, where the function to be minimized is linear and therefore easy to understand, and all the challenge is in the boundary constraints. Unconstrained nonlinear optimization is an old subject. Newton introduced the idea of approximating functions by the first few terms of what we now call their Taylor series; indeed, Arnol’d has argued that Taylor series were Newton’s “main mathematical discovery.” To find a zero x∗ of a function F of a real variable x, everyone knows the idea of Newton’s method: at the kth step, given an estimate x(k) ≈ x∗, use the derivative F (x(k)) to define a linear approximation from which to derive a better estimate x(k+1): x(k+1) = x(k) − F(x(k))/F (x(k)). Newton (1669) and Raphson (1690) applied this idea to polynomials, and Simpson (1740) generalized it to other functions F and to systems of two equations. In today’s language, for a system of n equations in n unknowns, we regard F as an n-vector whose derivative at a point x(k) ∈ Rn is the n × n Jacobian matrix with entries Jij (x(k)) = ∂Fi ∂xj (x(k)), 1 i, j n. This matrix defines a linear approximation to F(x) that is accurate for x ≈ x(k). Newton’s method then takes the matrix form x(k+1) = x(k) − (J(x(k)))−1F(x(k)), which in practice means that to get x(k+1) from x(k), we solve a linear system of equations: J(x(k))(x(k+1) − x(k)) = −F(x(k)). As long as J is Lipschitz continuous and nonsingular at x∗ and the initial guess is good enough, the convergence of this iteration is quadratic: x(k+1) − x∗ = O(x(k) − x∗2). (1) Students often think it might be a good idea to develop formulas to enhance the exponent in this estimate to 3 or 4. However, this is an illusion. Taking two steps at a time of a quadratically convergent algorithm yields a quartically convergent one, so the difference in effi- ciency between quadratic and quartic is at best a constant factor. The same goes if the exponent 2, 3, or 4 is replaced by any other number greater than 1. The true distinction is between all of these algorithms that converge superlinearly, of which Newton’s method is the prototype, and those that converge linearly or geometrically, where the exponent is just 1. From the point of view of multivariate calculus, it is a small step from solving a system of equations to minimizing a scalar function f of a variable x ∈ Rn: to find a (local) minimum, we seek a zero of the gradient g(x) = ∇f (x), an n-vector. The derivative of g is the Jacobian matrix known as the Hessian of f , with entries Hij (x(k)) = ∂2f ∂xi∂xj (x(k)), 1 i, j n
IV.21.Numerical Analysis 613 to Even the problem of stating local optimality a Hessian is always symmetric a matter involving LAGRANGE MULTIPLIERS [IIL64]and Though the Nev wton formulas for minimization and distinction between active and inactive constraints ion. One of the obstacles quickly encountered was 1951 and also twelve years earlier,it was subsequently meth ort ils if the initial gu addre sed both practically and theoretically by the active resea rch topic today. algorithmic technologles known as line searches and The problem of constraints brings us to the other rablems with mo e than a few variables.it also This subject was born in the 1930s and 1940s with quickly became clear that the cost of evaluating jace Kantorovich in the Soviet Union and Dantzig in the ns at every ste be e ir Fo inexact Jacobians or Hessians and/or inexact solutions the famous SIMPLEX ALGORITHM for solving lin of the still achie ear p program is nothing m e tha meth ables subicct to m ear equality and/or inegualin ods in the 196 0s by Broyden,Davidon,Flet her,and constraints.How can this be a challenge?On answer is part hat m and may b e pr or Hessian factor illustration and also in thei own right.A famous early mple h t at th time eon utput n ec hich lefinite quasi-N on updating formula was p the 1970s the Soviet Union used an inpu utout com fo ur dit nt autho puter mode l involving thousands of variables as a tool y has hoon knoun over since as the REcs made medium. and In subs the scale of linear progra ing probl tractabl ,ne .the )t a technology enables derivati tof satisfy all the const aints.Fo ute progra n th region is a p well as producing numerical outputs it also nal value of f is guaranteed to be attained at one of erivative he idea of automatic different rtex points. (A point is called vertex if it is s ne ts.)The ent of“reverse mode "form lations,it did not becom m one ver ully pra f.Carle,anc an optimal poin d.Al trained optimization problems are relatively 1984,an upheaval occurred in this field.trig easy.but they are the true depth of thi endra k f Bell Laborat es.K d fo dealin a functio than the s vorking in the f:RH- Ris to be e minimized subject o certain equa nterior of the fe asible region inst d.Once a conne and ida d hy
✐ IV.21. Numerical Analysis 613 and one may utilize it just as before in a Newton iteration to find a zero of g(x), the new feature being that a Hessian is always symmetric. Though the Newton formulas for minimization and finding zeros were already established, the arrival of computers created a new field of numerical optimization. One of the obstacles quickly encountered was that Newton’s method often fails if the initial guess is not good. This problem has been comprehensively addressed both practically and theoretically by the algorithmic technologies known as line searches and trust regions. For problems with more than a few variables, it also quickly became clear that the cost of evaluating Jacobians or Hessians at every step could be exorbitant. Faster methods were needed that might make use of inexact Jacobians or Hessians and/or inexact solutions of the associated linear equations, while still achieving superlinear convergence. An early breakthrough of this kind was the discovery of quasi-Newton methods in the 1960s by Broyden, Davidon, Fletcher, and Powell, in which partial information is used to generate steadily improving estimates of the true Jacobian or Hessian or its matrix factors. An illustration of the urgency of this subject at the time is the fact that in 1970 the optimal rank-two symmetric positivedefinite quasi-Newton updating formula was published independently by no fewer than four different authors, namely Broyden, Fletcher, Goldfarb, and Shanno; their discovery has been known ever since as the BFGS formula. In subsequent years, as the scale of tractable problems has increased exponentially, new ideas have also become important, including automatic differentiation, a technology that enables derivatives of computed functions to be determined automatically: the computer program itself is “differentiated,” so that as well as producing numerical outputs it also produces their derivatives. The idea of automatic differentiation is an old one, but for various reasons, partly related to advances in sparse linear algebra and to the development of “reverse mode” formulations, it did not become fully practical until the work of Bischof, Carle, and Griewank in the 1990s. Unconstrained optimization problems are relatively easy, but they are not typical; the true depth of this field is revealed by the methods that have been developed for dealing with constraints. Suppose a function f : Rn → R is to be minimized subject to certain equality constraints cj (x) = 0 and inequality constraints dj (x) 0, where {cj} and {dj} are also functions from Rn to R. Even the problem of stating local optimality conditions for solutions to such problems is nontrivial, a matter involving lagrange multipliers [III.64] and a distinction between active and inactive constraints. This problem was solved by what are now known as the KKT conditions, introduced by Kuhn and Tucker in 1951 and also twelve years earlier, it was subsequently realized, by Karush. Development of algorithms for constrained nonlinear optimization continues to be an active research topic today. The problem of constraints brings us to the other strand of numerical optimization, linear programming. This subject was born in the 1930s and 1940s with Kantorovich in the Soviet Union and Dantzig in the United States. As an outgrowth of his work for the U.S. Air Force in the war, Dantzig invented in 1947 the famous simplex algorithm [III.84] for solving linear programs. A linear program is nothing more than a problem of minimizing a linear function of n variables subject to m linear equality and/or inequality constraints. How can this be a challenge? One answer is that m and n may be large. Large-scale problems may arise through discretization of continuous problems and also in their own right. A famous early example was Leontiev’s theory of input–output models in economics, which won him the Nobel Prize in 1973. Even in the 1970s the Soviet Union used an input–output computer model involving thousands of variables as a tool for planning the economy. The simplex algorithm made medium- and largescale linear programming problems tractable. Such a problem is defined by its objective function, the function f (x) to be minimized, and its feasible region, the set of vectors x ∈ Rn that satisfy all the constraints. For a linear program the feasible region is a polyhedron, a closed domain bounded by hyperplanes, and the optimal value of f is guaranteed to be attained at one of the vertex points. (A point is called a vertex if it is the unique solution of some subset of the equations that define the constraints.) The simplex algorithm proceeds by moving systematically downhill from one vertex to another until an optimal point is reached. All of the iterates lie on the boundary of the feasible region. In 1984, an upheaval occurred in this field, triggered by Narendra Karmarkar at AT&T Bell Laboratories. Karmarkar showed that one could sometimes do much better than the simplex algorithm by working in the interior of the feasible region instead. Once a connection was shown between Karmarkar’s method and the logarithmic barrier methods popularized by Fiacco and