Numerics with automatic result verification

Mathematics and Computers in Simulation 35 (1993) 435-450 North-Holland

435

MATCOM 948

Numerics with automatic Ulrich

result verification

Kulisch

Universitiit Karlsruhe, Germany

L.B. Rall Mathematics Department, University of Wisconsin, Madison, WI, United States

Abstract Kulisch, U. and L.B. Rall, Numerics Simulation 35 ( 1993) 435-450.

with automatic

result verification,

Mathematics

and Computers

in

Floating-point arithmetic is the fast way to perform scientific and engineering calculations. Today, individual floating-point operations are maximally accurate as a rule. However, after only two or just a few operations, the result can be completely wrong. Computers now carry out up to 10” floating-point operations in a second. Thus, particular attention must be paid to the reliability of the computed results. In recent years, techniques have been developed in numerical analysis which make it possible for the computer itself to verify the correctness of computed results for numerous problems and applications. Moreover, the computer frequently establishes the existence and uniqueness of the solution in this way. For example, a verified solution of a system of ordinary differential equations is just as valid as a solution obtained by a computer algebra system, which as a rule still requires a valid formula evaluation. Furthermore, the numerical routine remains applicable even if the problem does not have a closed-form solution. This paper first illustrates some problems with ordinary floating-point arithmetic by means of some simple examples. Then, necessary extensions of computer arithmetic are entered into. Accumulation must be performed in fixed-point arithmetic. The continuum will be brought into the computer in the form of intervals. In order to have simple access to the new techniques, programming languages themselves have to be extended. ACRITH-XSC, Pascal-XSC and C-XSC are examples of such language extensions, which have been implemented with varying amounts of compiler effort. Result verification is carried out by means of mathematical fixed-point theorems such as Brouwer’s fixed-point theorem and its generalizations. Numerical computation with result verification is of fundamental significance for many applications, for example, for simulation or mathematical modeling. Models which are frequently developed by means of heuristic methods can only be refined systematically if computational errors can be essentially excluded. This paper gives an introduction to the whole field of numerical computation with automatic result verification, written with nonspecialists in mind. It extends from the design of hardware for arithmetic operations in VLSI to compilers, algorithms, numerical analysis, and deep into applications. A symbiosis of these new techniques with a computer system for symbolic computation would be highly desirable.

Correspondence to: Prof. L.B. Rall, Mathematics Department, University of Wisconsin-Madison, Vleck Hall, 480 Lincoln Drive, Madison, WI 53706-1313, United States. e-mail: [email protected].

0378-4754/93/S

06.00 @ 1993 -

Elsevier Science Publishers

B.V. All rights reserved

309 Van

436

U. Kulisch, L.B. Rail /Numerics with automatic result verification

1. Introduction Without question, methods of calculation have put their stamp on mathematics throughout history, and progress in connection with numbers has meant progress in mathematics also. In prehistoric times, people already counted and calculated. However, the beginning of real techniques for computation can be ascribed to the introduction of the Hindu-Arabic method of calculation and its dissemination through the famous arithmetic books of the thirteenth to the sixteenth centuries. Many mathematical results, such as the famous theory of solution of polynomial equations, are direct consequences of this development. The invention of mechanical calculating machines by Pascal, Leibniz and others, of logarithms and the slide rule, and also meanwhile of many other mechanical and analog computational aids, represent without doubt milestones not only in mathematics, but also in the history of the human race. A great leap in development took place about fifty years ago with the development of the first electronic computers. This technology immediately made it possible to perform arithmetic operations faster by a factor of about 1000 as compared to its mechanical or electromechanical predecessors. The great technical gains of this century would have been impossible without modern computational means. Today’s automobile, airplane, space travel, modern radio and television, and not least the swift further development of the computer itself were only possible on the basis of enormous computational capacity and bear visible witness to this development. On the other hand, this gave massive stimulation to further development of algorithmic and numerical mathematics. Further development of circuit technology, which again would not have been possible without the computer, has increased the computational capacity of a processor by a factor of about lo9 as compared to the first electronic computers of the fifties. Comparison of the numbers lo3 and lo9 shows that the real computer revolution took place after the development of the first electronic computers. Remarkably, there was nothing equivalent to this development on the arithmeticalmathematical side. The enormous advances in computer technology really suggested an attempt to make the computer also more powerful arithmetically. On this, however, mathematicians exerted hardly any influence. In the area of scientific and engineering applications of concern here, the computer is basically still used today in the same way it was in the middle fifties. Just the four floating-point operations of addition, subtraction, multiplication, and division are still being performed, only today many times more often and much faster. One has apparently become used to the known inadequacy and incompleteness of floating-point arithmetic, or else one supposes that it is a necessary evil. The consequences which this entails can be illustrated by two simple examples, which coincidentally were disseminated by the e-mail network during the same week of June, 1992. l In the recent Gulf War, a scud missile hit an American army barracks in Dhahran, killing 28 American servicemen and women. A Patriot missile was fired at the incoming scud, but did not hit it because of a numerical error [ 9 1. l In a recent flight of the space shuttle Endeavour, the robot arm could not recover a satellite automatically because of a numerical error. The astronauts had to climb out and haul the satellite into the shuttle by hand. In both cases, the manufacturers certified that the algorithms used were completely tested and approved, that is to say, were verified. Nevertheless, they failed numerically when used in actual situations. These examples show that it is certainly necessary and also compulsory to make numerics and the computer more trustworthy than in the case of exclusive use of ordinary floating-point arithmetic.

U. Kulisch, L.B. Rall /Numerics with automatic result verification

431

This development began about 25 years ago and has gained more and more ground in the meantime. However, there still remains a problem of acceptance to be considered. Even with only ordinary floating-point arithmetic, the computer in the form of a high-performance workstation of a vector processor or a supercomputer, or a massively parallel installation, is still an extremely powerful tool, and working with this tool permeates ordinary thinking. This creates a chasm which one prefers not to cross. “The cobbler sticks to his last.” Of course, a still more powerful but unknown tool on the other side cannot immediately be mastered or even set up and used effectively. Thus, the chasm is widened still more. This paper attempts to build a bridge to help reduce the difficulty of crossing it. To make it more accessible, it has been composed with as little mathematical formalism and as few technical details as possible.

2. Integer arithmetic As a rule, electronic computers now have two different number systems built into their hardware: integers and floating-point numbers. Integer arithmetic operates over a limited range of integers. In principle, this range can be extended indefinitely by software. As a rule, however, software slows down the computational process considerably. Inasmuch as the hardware (and also the software) are intact, integer arithmetic and its software extensions operate without error. For example, a rational arithmetic or a multiple-precision arithmetic can be constructed using integer arithmetic. In this way, the computer carries along as many digits as necessary to represent the result of an operation exactly. In number theory or computer algebra, this is the preferred method for arithmetic. Also, most applications in computer science can be carried out error-free with integer arithmetic. Examples are compilation or other forms of translation, and algorithms for searching or sorting.

3. Floating-point arithmetic On the other hand, in numerical or so-called scientific computation, one has to work with real numbers. As a rule, a real number is represented by an infinite decimal or binary fraction. This has to be approximated by a finite fraction in a computer. Ordinarily, this is done by the so-called floating-point numbers. Floating-point arithmetic is built into the hardware, and is consequently very fast. However, each floating-point operation is subject to error as a rule. Today, each floatingpoint operation is, as one says, maximally accurate. That is, the result of each operation differs at most by 1 or i unit in the last place according to the choice of the rounding function. However, after only two or just a few floating-point operations, the result can be completely wrong. We consider two simple examples. l Let x and y be two vectors with six components: x = (102’, 1223, IO’*, 10i5,3,-10i2)

and

y = (1020,2,-1022,

Let the scalar product be denoted by x .y. Here, the corresponding all the products are summed. This gives

1013,2111, 1016).

components

x .y = 104’ + 2446 - 104’ + 102* + 6333 - 102* = 8779.

are multiplied,

and

438

U. Kulisch, L.B. Rail /Numerics with automatic result verification

On the contrary, every computer (including those with IEEE arithmetic) gives the value zero for this scalar product. The reason for this is that the summands are of such different orders of magnitude that they cannot be processed correctly in ordinary floating-point format. This catastrophic error occurs even though the data (the components of the vectors) use up less than 4% of the exponent range of small and medium size computers! l For the second example, we consider a floating-point system with base 10 and a mantissa length of five digits. The difference of the two numbers x = 0.100 05 - lo5 and y = 0.999 73 . lo4 is to be computed. This time, both operands are of the same order of magnitude. The computer now gives the completely correct result x - y = 0.7700. 10’. Now, suppose that the two numbers x and y are the results of two previous multiplications. Since we are dealing with five-digit arithmetic, these products of course have ten digits. Suppose the unrounded products are x1 ‘x2 = 0.100054824

1 . 105,

yl . y2 = 0.999 7342213.

104.

Subtracting these two numbers, after normalization and rounding to five places, one obtains the result 0 (xi . x2 - y1 . ~2) = 0.8 14 02 . 10’) which differs in every digit from the result computed above. In the second case, the value of the expression xi . x2 - yl . y2 was computed to the closest five-digit floating-point number. On the contrary, pure floating-point arithmetic with rounding after each individual operation gave a completely wrong result. It is plainly an almost inexplicable failure that this inadequacy has not yet been done away with after forty years of floating-point arithmetic and enormous advances in technology. One should be very, very careful with the argument “nothing like that happens in my calculations”. It is very difficult to have a view of all the data that enter as intermediate results into an hour-long computation on a workstation or a supercomputer. Many computer users appear to be like a wood carver who is forced to carve a work of art with a dull knife, even though it is so easy to sharpen a knife! Computers become ever faster. In the middle fifties, computers were able to carry out about 100 floating-point operations per second. Today, fast computers are able to execute some billion floating-point operations per second, and this in situations where some computations take hours. In the meantime, the data format used has hardly changed. One computes with about seventeen decimal digits. Sometimes, a format twice as long is used in cases of extremely critical applications. In classical error analysis of numerical algorithms, the error of each individual floating-point operation is estimated. It is evident that this is no longer practically possible for an algorithm for which 1014 floating-point operations are performed in an hour. Thus, an error analysis is no longer carried out as a rule. Indeed, the fact that the computed result could be wrong is often not taken into consideration. Once in a while, the user tries to justify the computed result by heuristic methods, or make it plausible. However, there are also users who are not aware of the fact that the computed result could also be completely wrong, perhaps due to confusion with the completely different properties of computation with integers. From the mathematical point of view, the problem of correctness of computed results is particularly of central importance because of the high computational speeds attainable today. The determination of the correctness of computed results is essential in many now classical applications such as simulation or mathematical modeling. This enables the user to distinguish between inaccuracies in computation and actual properties of the model. A mathematical model can only be developed systematically if errors entering into the computation can be essentially excluded.

U. Kulisch, L.B. Rail/Numerics

with automatic result verification

439

4. Necessary extensions of ordinary floating-point arithmetic In the past, what were regarded as scientific computations were carried out using a slide rule or tables of logarithms of numbers and standard functions. The use of logarithms avoids the necessity of keeping track of exponents, as when using a slide rule, but in either case precision is limited. Furthermore, additions and subtractions had to be carried out as side calculations in which terms of relatively small magnitude were neglected, and subtraction of almost equal terms could render the results useless. (The use of Gaussian logarithms for addition and subtraction still results in limited precision, neglect of terms of small magnitude, and the possibility of cancellation errors.) When done by hand, the number of possible computations was so small that the effects of neglecting small terms and cancellation could actually be observed and appropriate measures taken. A modern mathematics coprocessor chip can be viewed as automating computation using logarithm tables with a large (but fixed) number of places, including addition and subtraction logarithms. It is so fast, however, that there is no way to observe the cumulative effect of neglect of small quantities, or cancellation errors. By contrast with early scientific computation, which can be characterized as multiplicationdivision intensive, accountants and bookkeepers labored with addition and subtraction of long columns of numbers of various magnitudes. Errors in these addition-subtraction intensive computations were not tolerated, and various methods of validation (such as double-entry bookkeeping) were developed to ensure that answers came out correct to the penny. In order to handle potentially millions of additions and subtractions accurately, it is evident that the “slide-rule whiz” capabilities of ordinary floating-point arithmetic have to be augmented by the abilities of a patient accountant. Given the state of the art, there is no reason why this cannot be done on a single chip. Computers were invented to take on complicated work for people. The evident discrepancy between computational power and control of computational errors suggests also turning over the process of error estimation itself to the computer. This has been successfully done in the meantime for practically all the basic problems of numerical analysis and many applications. To achieve this, the computer has to be made arithmetically more powerful than is ordinarily the case. One thus proceeds from the above observations. In floating-point arithmetic, most errors occur in accumulations, that is, by execution of a sequence of additions and subtractions. On the other hand, multiplication and division in floating-point arithmetic are relatively stable operations. In fixedpoint arithmetic, however, accumulation is performed without errors. Thus, most errors encountered in floating-point arithmetic can be avoided if accumulation is carried out in fixed-point arithmetic. The present state of technology allows an easy realization of fixed-point accumulation. One only has to provide a fixed-point register in the arithmetic unit which covers the whole floating-point range. For classical /370_architecture, in the data format long it is sufficient to employ a fixed-point register of exactly 568 bits or 71 bytes. If such a hardware register is not available, then it can be simulated in the main memory by software. This naturally results in loss of speed as a consequence, which in many cases is outweighed by the gain in certainty. If this register is made twice as long, which is likewise easily possible, then dot products of vectors of any finite dimension can likewise always be computed exactly. (The products, that is, the summands in the dot product, have double mantissa lengths, and the exponent range is likewise doubled.) In the case of the frequently used data format of the so-called IEEE arithmetic standards, about four times as many bits or bytes are necessary because of its substantially greater exponent range, likewise presenting no problem to modern technology. The possibility of computing the

440

U. Kulisch, L.B. RaN/Numerics with automatic result verification

dot product of floating-point vectors of any finite dimension with complete exactitude opens a new dimension in numerical analysis. In particular, the optimal dot product proves itself to be an essential instrument in attaining higher computational accuracy. Now, it is also evident that with the dot product, all operations between floating-point numbers, as well as those for complex numbers 1 and, in particular, all operations on vectors or matrices with real or complex components can be carried out exactly as well as with maximal accuracy, that is, with only a single rounding in each component. Only a few years ago, just the error-free computation of the product of two floating-point matrices was considered to be harder than the calculation of the eigenvalues of a symmetric matrix! In the computer proposed here, this will be a basic arithmetic operation. In addition, higher precision floating-point arithmetic can be based on the optimal dot product. In order to adapt the computer also for automatic error control, its arithmetic must be extended with still another element. All operations on floating-point numbers, which are addition, subtraction, multiplication and division, and the dot product of floating-point vectors must be supplied with directed rounding, that is, rounding to the nearest smaller and larger floating-point numbers. An interval arithmetic for real and complex floating-point numbers, as well as for vectors and matrices with real and complex floating-point components can be built with these operations. Intervals bring the continuum into the computer and likewise open a completely new dimension in numerical analysis. An interval is represented in the computer by a pair of floating-point numbers. It describes the continuum of real numbers bounded between these two floating-point numbers. Operations on two intervals in the computer result from operations on two appropriately chosen endpoints of the operand intervals. In this, the computation of the lower endpoint of the result interval is rounded downward, and the computation of the upper endpoint is rounded upward. The result is certain to contain all results of the operation considered applied individually to elements of the first and second operand intervals. The interval evaluation of an arithmetic expression such as a polynomial costs about twice as much as the evaluation of the expression in simple floating-point arithmetic. It provides a superset of the range of the expression or function in the given interval of arguments.

5. Computation with automatic result verification for algebraic problems With these two arithmetic aids, an essential basic principle of computation with automatic result verification can already be described for algebraic problems. We will illustrate this by two simple examples. We first address the question of whether a given real function defined by an arithmetic expression has no zeros in a given interval X. This question cannot be answered with mathematical certainty if only floating-point arithmetic is available. One can evaluate the function say 1000 times in the interval X. If all computed function values turn out to be positive, then the probability is high that the function does not have any zeros in X. However, this conclusion is certainly not reliable; because of roundoff errors, a positive result could be computed for a negative function value. The function could also descend to a negative value which was not detected due to the choice of the evaluation points. On the other hand, a single interval evaluation of the function may suffice to solve this problem with complete mathematical certainty. If the computed interval result

‘For complex division, an additional consideration is necessary.

U. Ku&h,

L.B. Rail /Numerics with automatic result verification

441

does not contain zero, then the range of values also cannot contain zero, because the computed interval result is a superset of the range. As a consequence, the function does not have any zeros in the interval X. The interval evaluation of the function is only about twice as costly as a single floating-point evaluation of the function. As a second example, we sketch a method by which one can verify the correctness of the computed solution in the case of a linear system of equations. First, an approximation for the solution is computed by some ordinary method, for example by Gaussian elimination. This approximation is expanded by augmenting it with a small number & E in each component direction. This results in an n-dimensional cube X which has the computed approximation as its midpoint. This cube is just an interval vector consisting of intervals in each coordinate direction. Now, one transforms the system of equations in an appropriate way into fixed-point form x = A. x + b, where one naturally makes use of the fact that one has already computed a presumed approximate solution. Then the image Y : = A . X + b is computed. If Y is contained in X, which can be verified simply by comparison of endpoints, then Y contains a fixed point of the equation x = A. x + b by Brouwer’s fixed-point theorem2. If Y is strictly contained in X, so that the endpoints do not touch, then the fixed point is unique, that is, it has been verified computationally that the original matrix of the system of equations is nonsingular. By iteration, using the optimal dot product, the fixed point can then be determined to arbitrary accuracy in practice. In this way, the user obtains a mathematically certain and reliable result. The two examples of numerical computation with automatic result verification described here leave open the question of what to do if the verification step fails. In this case, refined methods have of course been developed which allow valid assertions. However, these will not be gone into here. A great advantage of automatic result verification is that the computer itself can quickly establish that the computation performed has not led to a correct and usable result. In this case, the program can choose an alternative algorithm or repeat the computation using higher precision. Similar techniques of automatic result verification can be applied to many other algebraic problem areas, such as the solution of nonlinear systems of equations, the computation of roots, the calculation of eigenvalues and eigenvectors of matrices, optimization problems, etc. In particular, valid and highly accurate evaluation of arithmetic expressions and functions on the computer is included. These routines also work for problems with interval data. The evaluation of a function is reduced to the solution of a simple system of equations, an approximate solution of which is corrected if necessary by defect correction and the optimal dot product to a guaranteed, maximal accuracy. This technique of enclosing function values with high accuracy is an essential aid for numerics. Newton’s method for the computation of the roots of a system of equations is frequently unstable in a neighborhood of a root. That is to say, a root results from the mutual cancellation of positive and negative terms. In the immediate vicinity of a root, these terms are only approximately equal and cancellation occurs in ordinary floating-point arithmetic. However, in the case of insufficiently accurate function evaluation, Newton’s method can easily overshoot its goal, possibly into the domain of attraction of another root. An iteration process which converges to a solution x in the real numbers may thus not converge to the value x in ordinary floating-point arithmetic. It can swing over to an entirely different root. A conclusive

2Brouwer’s fixed-point theorem asserts the following. Let CJ~ : R” -+ R” be a continuous mapping, X C U!”be closed, convex and bounded, and Y : = c$(X). If Y S X, then 4 has at least one fixed point in Y.

442

U, Kulisch, L.B. RaN /Numerics with automatic result verification

verification step is thus also indicated for iteration processes. A useful and often applied method to strive for high accuracy consists of enclosing the error of an approximation instead of the solution itself. The techniques sketched here for numerics with automatic result verification use the optimal dot product and interval arithmetic as essential tools in addition to ordinary floating-point arithmetic. Interval arithmetic permits computation of certain bounds for the solutions of a problem. One obtains high accuracy by means of the optimal dot product. The combination of both instruments is indeed a breakthrough in numerical analysis. Naive application of interval arithmetic alone certainly always leads to reliable bounds. However, these can be so wide that they are really useless. This effect was observed twenty years ago by many numerical analysts, leading frequently to the rejection once and for all of this useful, necessary and blameless tool.

6. Necessary

extensions of ordinary programming languages

In putting the new techniques and possibilities into practice, a larger difficulty is encountered at once. The usual programming languages such as ALGOL, FORTRAN, BASIC, Pascal, MODULA or C only provide the four basic operations of floating-point arithmetic for computing with real numbers. They do make it possible to write down simple expressions, formulas or functions of reals in the customary way. However, it is not practicable to call a maximally accurate dot product or a maximally accurate matrix multiplication in these languages. If a matrix multiplication, say, is programmed in the usual way by three nested loops, then the dot products which occur in them are again evaluated in ordinary floating-point arithmetic. We have already pointed out above that this is to be unconditionally avoided. Operations in product spaces, as in the case of vectors and matrices, can only be realized awkwardly in the languages cited by procedures, in which the programming practically descends to assembly level. In mathematics and science, the concept of arithmetic operations or of functions is by no means limited to the real numbers. Among others, the operations in ordinary vector spaces and also vector-valued functions are the most clumsy to realize and it is really not sensible always to have to reduce these operations to the elementary floating-point operations and then use awkward procedure calls. This way of doing things leads, as we have seen, to many unnecessary computational errors. The only reasonable and sensible way out of this difficulty consists of extending these languages so that the usual mathematical operator symbols for the corresponding predefined data types are provided for all operations in real and complex vector spaces and the corresponding interval spaces. If, say, A, B, C are variables which represent real matrices, then matrix multiplication is denoted simply by the assignment C := A * B. Here, the operation sign * calls at the machine level an elementary operation possibly realized in hardware which computes each component of the product matrix with maximal accuracy. This so-called “operator notation” for practically all ordinary arithmetic operations simplifies programming considerably. Many problems at once become really tractable and reliably programmable by means of operator notation. Programs become considerably easier to read. They become easier to test and are thus essentially more reliable. However, above all, the operations are maximally accurate.

V. Ku&h,

L.B. Rail/Numerics

with automatic

result verification

443

7. Availability of suitable language extensions

Already in the middle sixties, an ALGOL extension was created and implemented [7] which had in addition to types integer and real a type for real intervals including provision of the corresponding arithmetic and relational operators. This language served as an ancestor for the creation of ALGOL-68, and influenced the provision of an operator concept in it. Except for the necessary arithmetic, ALGOL-68 had all the necessary language concepts for implementation of numerical algorithms with automatic result verification. Unfortunately, ALGOL-68 did not catch on, perhaps because computers were not yet powerful enough. Subsequent languages which provide similar language capabilities are Ada, C+ + and FORTRAN-90. Ada has not yet caught on so well, perhaps because it is not upward compatible with any of the programming languages used in the meantime, so that all programs have to be written anew. Fast compilers for FORTRAN-90 are just now under development. At present, there are three language developments and corresponding compilers which contain all programming and arithmetic concepts for programming numerical routines with automatic result verification. These are the programming languages ACRITH-XSC, Pascal-XSC and C-XSC. They were all developed at the Institute for Applied Mathematics at the University of Karlsruhe, or in collaboration with this institute. ACRITH-XSC is an extension of the widely-used programming language FORTRAN-77. Together with a large number of problem solving routines with automatic result verification, ACRITH-XSC was developed and implemented at the Institute for Applied Mathematics of the University of Karlsruhe with the sponsorship and collaboration of IBM. It is now distributed worldwide as an IBM Program Product. Unfortunately, the use of ACRITH-XSC is limited to machines with IBM /370-architecture and here again only to those which operate under the VM CMS operating system. Pascal-XSC and C-XSC are extensions of the useful and widely-distributed programming languages Pascal and C, respectively. Compilers for Pascal-XSC and C-XSC are provided for a large number of machines of various manufacturers. The language description for Pascal-XSC has been published by Springer in German and English. A Russian edition is in preparation. A language description of C-XSC in English is likewise published by Springer. The language description also contains advice pertaining to application of the compiler 3, including some problem solving routines with automatic result verification. In addition to a large number of predefined arithmetic data types and the corresponding operations with maximal accuracy, the XSC-languages include among others the following language concepts, which were not present in the relatively basic version of the language. A module concept, an operator concept, functions and operators with general result types, overloading of functions, procedures, and operators, overloading of the assignment operator as well as read and write, dynamic arrays, access to subarrays, control of rounding by the user, and exact evaluation of expressions. For example, with these language concepts, modules and operators for complex arithmetic, matrix and vector arithmetic, rational arithmetic, double- or multiple-precision arithmetic or interval arithmetic can be developed in the language itself. Arithmetic expressions and algorithms with these data types can be written down in ordinary mathematical notation. This simplifies programming very greatly. Many programs can be read like a technical report. See [ 5 ] with regard to this. 3Compilers for Pascal-XSC and C-XSC are obtainable from Numerik Software GmbH, Postfach 2232, D7570 Baden-Baden, Germany, fax (+49) 721-694418, phone (+49) 721-694217, 44666, Madison, WI 53744-4666, United States, phone (608) 273-3702.

or FBSoftware,

P.O. Box

444

U. Kulisch, L.B. Rall /Numerics with automatic result verification

8. Automatic differentiation Automatic differentiation denotes a method which permits comput&ion of the value of the derivative of a function at a given location as well as the value of the function. This location can also be an interval. In this case, bounds for the derivative in this interval are computed. Further developments of automatic differentiation include the automatic generation of Taylor coefficients and gradient vectors and Hessian matrices of functions of several variables. Explanation of basic techniques of automatic differentiation and an overview of the state of the art are given in the proceedings of a recent conference on the subject [2]. Automatic differentiation and automatic generation of Taylor coefficients are essential tools for the development of problem solving routines with automatic result verification for problems in numerical analysis, such as for example numerical quadrature or the numerical integration of differential and integral equations. Automatic differentiation makes essential use of the programming language extensions mentioned in the previous section, especially the operator concept and operator overloading. Automatic differentiation is completely different from numerical differentiation, say by difference quotients. It also does not make use of symbolic differentiation of the expression. For the automatic differentiation of an arithmetic expression for a function, a number pair is substituted for each of the operands of the expression (which are constants, variables and standard functions), where the first component gives the value of this operand and the second the value of its derivative at the point at which they are evaluated. These number pairs are then combined with each other according to the operations given in the expressions by the known rules for differentiation of sums, differences, products, quotients, and the chain rule, as well as the rules for differentiation of standard functions. A number pair is again obtained as the result, the first component of which is the value of the expression and the second is the value of its derivative at the given location. For the automatic generation of Taylor coefficients, one computes with corresponding n-tuples, where the first component gives the value of the function at the location considered, while the further components give the values of its first, second,. . . Taylor coefficients. It is important to note that computations are done with numbers in each case, and not with symbolic expressions. Automatic differentiation and automatic generation of Taylor coefficients for somewhat complicated expressions can hardly be carried out by hand, since the chain rule or the quotient rule, say, have relatively complicated formulations, so that one regularly makes mistakes very quickly. However, a machine can carry out these complicated formulations very stably, quickly and correctly. Programming of automatic differentiation or automatic generation of Taylor coefficients turns out to be very simple with the programming language extensions described in the previous section, particularly overloading of operators, functions and procedures. The expression is simply written in ordinary notation, and the corresponding operations are carried out according to the rules for automatic differentiation of pairs or n-tuples. The machine furnishes the values of the function and its derivative or higher Taylor coefficients. In the case of interval data, it furnishes inclusions of these values. Automatic differentiation is carried out very quickly and efficiently on the computer. One distinguishes between forward and reverse methods. With the reverse method, for example, the expense of computation of a gradient is of the same order of magnitude as the calculation of the function.

U, Kulisch, L.B. Rall /Numerics

with automatic

result verification

445

9. Numerics with automatic result verification for problems in analysis Automatic differentiation is an essential tool in the development of methods with automatic result verification for problems in the area of numerical analysis, for example, numerical quadrature and the integration of differential or integral equations. The essential distinction between methods with automatic result verification and ordinary methods for numerical integration is that in the latter the remainder term is commonly neglected, while in the former it is completely included in the evaluation. Ordinarily, higher derivatives to be evaluated at an unknown intermediate point enter into remainder terms. Such an unknown intermediate point is simply replaced by an interval which contains the intermediate point. By means of interval computation and automatic differentiation, an inclusion can be obtained for the higher derivative which enters into the remainder term. By means of a multiplier (such as P/n!) likewise contained in the remainder term, the truncation error expressed by the remainder term can usually be reduced below the desired accuracy of the result. By expression of the remainder term and use of interval computation in this way, one can bound the truncation error and the roundoff error as well, and thus be led to mathematically completely valid assertions. Very efficient methods have been developed for numerical quadrature and cubature with automatic result verification on the basis of extrapolation methods. For each element, a representation for the remainder term is derived from the extrapolation tableau calculated in the ordinary way on the basis of the Euler-Maclaurin summation formula. The integration begins with evaluation of this remainder term for a given element of the tableau. If the resulting value lies inside prescribed error bounds, then follows the actual computation of the value of the integral as a linear combination (dot product) of the function values at the nodes and the weights which were previously calculated once and for all and stored in the computer. Very efficient methods have also been developed for integral equations, which are distinguished by their universality. For example, for Fredholm integral equations of the second kind, the kernel can be developed into a two-dimensional Taylor series with remainder term by automatic differentiation. It will then consist of the sum of a finite-rank kernel and a contractive kernel. Both parts can then be dealt with on the basis of known methods, and thus the solution can be enclosed. This method has been applied successfully to very large systems of nonlinear integral equations. By transformation into integral equations, problems in partial differential equations have also been dealt with successfully with automatic result verification. Conversely, problems with integral equations have been transformed into differential equations, and then solved with automatic result verification. In the case of ordinary differential equations, there are methods for initial-value problems as well as for boundary-value problems and eigenvalue problems. For initial-value problems, say in the case of a first-order system, the right side of the differential equation is developed into a Taylor series by automatic differentiation. This is followed by evaluating an inclusion of the remainder term for given step size. In each subinterval, continuous lower and upper bounds for the solution are obtained. Thus, in each individual step, one has information about the accuracy, which can be improved by marching backward and/or the choice of a smaller step size. Boundary-value and eigenvalue problems can be converted to initial-value problems by shooting methods, for example. Among others, great successes have been achieved in computing periodic solutions of differential equations. The wrapping effect for initial-value problems known from the literature can be controlled by introduction of local coordinates. In the case of differential and integral equations, methods

446

U. Kulisch, L.B. Rail /Numerics with automatic result verification

with automatic result verification furnish piecewise continuous upper and lower bounds for the solution. In the case of boundary-value and eigenvalue problems, the existence and uniqueness of the solution is verified by the algorithm, that is, with the computer. Experimentation with verified solutions of differential equations proves to be very interesting. For example, a “computed solution” can move on a completely different path if only two places of accuracy are lacking at a critical location, which naturally cannot be recognized if an ordinary approximate method is used. In the case of a verified method, the width of the inclusion will explode at such a location, which indicates that the accuracy must be entirely increased. In recent years, reliable solutions have been computed for a number of engineering and scientific problems with the problem solving routines with automatic result verification for basic problems of numerical analysis provided in the programming languages Pascal-XSC, ACRITH-XSC and C-XSC. Frequently, problems for which ordinary solvers and solution methods fail are handled. GAMM (Gesellschaft fur Angewandte Mathematik und Mechanik) and IMACS (International Association for Mathematics and Computers in Simulation) have held symposia on the theme “Computer Arithmetic and Scientific Computation with Automatic Result Verification” regularly since 1980. Proceedings of many of these meetings have appeared. Applications belong to the widest variety of scientific fields. By illustration, some of the problems handled successfully are listed here: highspeed turbines, filter calculations, plasma flows in a rectangular channel, gear-drive vibrations, high-temperature superconductivity, infiltration of pollution into ground water, inductive drives for buses, periodic solutions of the oregonator (chemical kinetics), geometric modeling, dimensioning of sewers, verified quadrature in chemical transport engineering, nonlinear pile driving, refutation of certain “chaotic” solutions in the three-body problem of celestial mechanics, development of solutions of the Schrodinger equation in wave functions, calculation of magnetohydrodynamic flows for large Hartmann numbers, optical properties of liquid crystal cells, simulation of semiconductor diodes, optimization of VLSI circuits, rotor oscillations, geometric methods for CAD systems, centered beam fracture, robotics. An extensive bibliography is found in [ 11.

10. Present availability of the necessary basic arithmetic in hardware In conclusion, we will make some additional comments on the state of developments in the field of the necessary arithmetic hardware. It is an interesting fact that the better mechanical calculators of 100 years ago already had all the arithmetic capabilities that are needed today. Above all, they were always able to calculate dot products exactly, that is, sums of products. The techniques used then can be considered to be direct predecessors of the electronic circuits available today. One such calculator had a keyboard of say nine decimal digits. The result register was essentially larger, and had perhaps 25 decimal digits. Above the input register was a movable carriage. It was a fixed-point register, that is, the decimal point always stayed in the same place during a calculation. This assigned a unique exponent to each digit. Thus, it was always possible to add a large number of summands positioned correctly in the keyboard. Products were accumulated by adding together multiples of the multiplicand corresponding to the digits of the multiplier. One could thus add a thousand products without further ado. As long as overflow or underflow of the result register did not occur, all digits of the result were correct. This operation was called simply a “running total”. It was used as often as possible, because it was the fastest way to carry out calculations. Intermediate results did not

U. Kulisch, L.B. Rall /Numerics with automatic result verification

447

need to be read out and then read in again for subsequent operations. Rounding was done by hand as needed downward, upward or to the closest number. Thus, we would like today to have the computer arithmetically as it already was 100 years ago. Apparently, the extremely useful, fast and accurate running total was lost in the transition to the first electronic computers. In addition, a completely wrong way was taken in the case of introduction of the vector processor. In order to gain speed, accumulation is not done simply in a fixed-point register, but instead in a pipeline and again in floating-point arithmetic. In this way, the sequence of summands is altered, which gives rise to additional computational errors. One really has to wonder why this has not given rise to an outcry from all mathematicians. Customers should be urgently advised not to overlook such defects when purchasing such installations, which generally cost several million dollars, and instead send in corresponding requirements for improvement to the manufacturers. For example, the organizations which fund civilian and military research in various countries finance purchases of large numbers of computers, and could well lend the necessary weight to such requirements. In the meantime, corresponding requirements have been set up by international organizations such as GAMM and IMACS. However, it is evident that these standards will be adopted only if backed up by large orders for conforming equipment and thus produce substantial income for manufacturers. The entirely correct accumulation of floating-point numbers as well as products of such numbers in a fixed-point register required here will likewise speed up the computer considerably. This kind of accumulation is really considerably simpler and thus essentially faster than an accumulation in floating-point arithmetic. For all kinds of computers such as personal computers, workstations, mainframes and supercomputers, the electronic circuitry has been developed so that arithmetic essentially contributes nothing to computation time. In the calculation of a dot product in a pipeline (in serial computation), a multiplication, a possible shift and the corresponding addition can be carried out without further ado in the time required for reading the two operands of the product. In contrast to this, in the case of accumulation in floating-point arithmetic after each individual addition and multiplication there is normalization and rounding, the latter giving rise to computational error, which are completely unnecessary and better omitted to speed up the course of the computation. Actually, the definition of computer arithmetic on the basis of mathematical criteria should be an inherent feature of every programming language. Users must know exactly what happens when they call an arithmetic operation in their program. Only in this way will the computer be a mathematical instrument. At the creation of the first programming languages in the fifties, no simple possibility for the definition of computer arithmetic was known. Thus, this complication was simply ignored then, and implementation of arithmetic was left to the manufacturer. Thus, a massive confusion resulted, which created specialists who detected all kinds of pitfalls in the arithmetic of individual computers. Basically, such a definition is very simple. One only needs to require that for a given data format and data type, the arithmetic operations given by the hardware or callable from the programming language satisfy the conditions of a semimorphism [ 61. This can be written down mathematically in four lines. All operations defined in this way, including those in product spaces, are then maximally accurate. That is, the computed result differs from the correct result by at most one rounding. The techniques of implementing these mathematically strict conditions differ only marginally from what is done now, that is, they are not essentially more complicated. If the vector and matrix operations are built into the design of the hardware, then the computer will consequently be a good bit faster.

448

W. Kulisch, LB. Rail/Numerics

with automatic result verijkation

The proposals of the IEEE, ANSI and IS0 microprocessor arithmetic standards 754 and 854 promulgated around 1985 are basically a step in the right direction. The operations given in them realize a semimorphic floating-point and interval arithmetic. Unfortunately, except for the XSC languages mentioned in Section 6, there are still no programming languages and compilers available which support simple use of interval arithmetic. This is all the more lamentable since prototypes for them were already available 25 years ago [ 7 1. Nevertheless, the initiators of the arithmetic standards mentioned cannot be spared a strong reproach. By 1985, technology had finally advanced so far and the necessary mathematical insight was available, so that one should not have stopped with a standardization of ordinary floating-point arithmetic. The always correct accumulation in a fixed-point register should have been included unconditionally. Now, many manufacturers believe with respect to arithmetic that realization of the microprocessor standards is all that is necessary to do. In this way, the existing standards also prove to be a great hindrance to further progress. An already outmoded concept has been fixed in place for many years by the standards. The simple examples mentioned in the introduction for the refutation of numerical algorithms should give thought to those mathematicians who continually assert that every solvable problem can be solved satisfactorily by ordinary floating-point arithmetic. The tragedy of this argument is that it is even correct; if one has determined that a computation has broken down, there is practically always a solution which leads to a correct result. In the end, one could indeed even use a universal computer as a Turing machine! The unsatisfactory and basically treacherous aspect of the present situation is that ordinary floating-point arithmetic gives no indication of whether it has broken down or not. However, in many applications it is important to obtain a correct solution at the first attempt and not have to use all possible skill and conjurer’s tricks after the computation has broken down, when the error that occurred could have been avoided. One should not first make wrong what can be made right from the very beginning! In the field of arithmetic, accumulation is the prime example of this. The programming languages Pascal-XSC, ACRITH-XSC and C-XSC cited in Section 6 bring in an essential advance here. They provide a universal computer arithmetic and permit simple handling of the semimorphic operations by means of ordinary mathematical operator symbols also in product spaces such as of complex numbers, as well as for vectors and matrices of types real, complex, interval and complex interval. Unfortunately, the latter operations are not supported by the IEEE Arithmetic Standards 754 and 854 mentioned above, so they have to be simulated in software. Thus, at a place where an increase in speed is really to be expected, a loss of speed results into the bargain. An additional difficulty comes in because processors which offer the IEEE Standard 754 do not furnish the double-length product, so this also has to be simulated, which means an additional considerable loss of speed. Double-length products are indispensable and essential for the computation of semimorphic vector and matrix multiplication. Multiple-precision arithmetic for various basic types can be realized very easily with the optimal dot product. Many programs for computer algebra can be speeded up very considerably by hardware support of the dot product. With this, the paper [5] also illustrates how very complex calculations can be programmed, read and executed very clearly and simply. In conclusion, it will be mentioned how far the required arithmetic properties have be realized in hardware. In 1982-l 983, IBM furnished all required operations on machines with /370 architecture (with the exception of the Vector Facility), on the larger machines by assembly software, and by

U. Kulisch, L.B. Rail/Numerics

with automatic result verification

449

microcode and later in VLSI on the smaller machines. These operations are addressable in the programming language ACRITH-XSC, a FORTRAN-77 extension, which is offered world-wide as an IBM program product. Unfortunately, ACRITH-XSC can be used only on IBM /370 machines which run under the operating system VM CMS. On more modern IBM architecture, the RS 6000 as well as the IBM PC, no hardware support for the optimal dot product is available. Following the development by IBM, Siemens apparently commenced offering the optimal dot product in microcode in about 1985 on all the computers with /370 architecture developed and built by Siemens themselves. However, support of these operations by means of a programming language is lacking up to now, so they are again of limited usefulness. Also in the middle of the eighties, Hitachi and NAS followed IBM and supported all the required operations on some of their machines in hardware (probably by microcode). Here the arithmetic can sometimes be addressed by the programming language ACRITH-XSC, because some of these machines operate under the IBM operating system VM CMS. Unfortunately, these advances achieved in the eighties have nearly dissolved into nothingness, since for some time scientific computation has already shifted to workstations and supercomputers. Here, there is no hardware support available at present, so that again one has to resort to software simulation. The situation with regard to supercomputers has already been mentioned above. Likewise, it is apparently up to customers to require advances in workstations. We would like this paper to help to strengthen the insight that an essential step forward has to be made in the field of computer arithmetic, as well as in computer hardware and also in standards. Basically, it is not comprehensible why guaranteed maximal accuracy is required of the four fundamental operations for real variables, but not for complex variables, vectors and matrices. The additional hardware cost is small in the case of a unified design of the arithmetic unit. However, it makes a big difference if one of these operations always gives a correct result, or only frequently. In the latter case, the user has to make an account of the desired accuracy in the case of each application, and certainly under conditions of many million times per second. In the first case, this is never necessary. In conclusion, to recall the theme of this paper in retrospect once more, it does not follow in any way that if an algorithm which has been verified to be correct in the sense of computer science, then a numerical result computed with it will be correct. In this connection, again recall the examples furnished. Here, numerics with automatic result veriIication is an important additional yet simple tool. Naturally, with it alone, one is not absolutely certain, because the verification step could still lead to a positive result because of a programming error. Verification of the algorithm, compiler, the operating system and also the computer hardware can never be made superfluous by it. In case of success, numerical result verification at least implies with high probability that all these components have furnished a correct result. Additionally, this result is certainly generally independent of whether a program verification has been carried out or not. On the other hand, a negative outcome of result verification, that is, when the computed result cannot be verified to be correct, also has positive value as a rule. That is to say, this result is established quickly and automatically by the computer, without anything having to be done externally. The user can provide program alternatives in this case, for example, choice of different algorithms or methods, or a change to higher precision. Thus, something can often succeed. In this way, correct and verified solutions have been computed for extremely unstable cases as well. For more details and references to this wide-ranging subject, see 11 I.

450

U. Kulisch, L.B. RaN /Numerics with automatic result verification

References [ 1 ] E. Adams and U. Kulisch, eds., Scientijc Computing with Automatic Result Verification (Academic Press, New York, 1993). [2] A. Griewank and G.F. Corliss, eds, Automatic Differentiation of Algorithms, Theory, Implementation, and Application (SIAM, Philadelphia, PA, 1992). [ 31 R. Klatte, U. Kulisch, Ch. Lawo, M. Rauch and A. Wiethoff, C-XSC (Springer, Berlin, 1992). [4] R. Klatte, U. Kulisch, M. Neaga, D. Ratz and Ch. Ullrich, PASCAL-XX’, Language Description with Examples (Springer, Berlin, 199 1) (in German); (Springer, New York,, 1992) (in English). [5] W. Kramer, Numerische Berechnung von Schranken fiir Pi, Jahrb. Uberblicke Math. 1993 (Vieweg, Braunschweig, 1993 ). [ 61 U. Kulisch and W.L. Miranker, Computer Arithmetic in Theory and Practice (Academic Press, New York, 1981). [ 71 H.-W. Wippermann, Realisierung einer Intervallarithmetik in einem ALGOL-60 System, Elektron. Rechenanlagen 9 (1967) 224-233. [ 81 ACRITH-XSC, IBM High Accuracy Arithmetic-Extended Scientific Computation. Version 1, Release 1, IBM, 1990; ( 1) General Information Manual (demonstration disc available), GC33-6461-01; (2) Reference, SC33-6462-00; (3) Sample Programs, SC33-6463-00; (4) How to Use, SC33-6464-00; (5) Syntax Diagrams, SC33-6466-00. [9] Patriot Missile Defense, United States General Accounting Office, GAO, IMTEC-92-26, 1992.