Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 88
Текст из файла (страница 88)
The definition of p isgeneralized and a method of computing it is developed that involves solvinga generalized eigenvalue problem. The book by Miller and Wrathall [762,1980] gives a thorough development of the work of [760, 1978], includinga description of the graph theory techniques used to compute the partialderivatives, and it provides further examples of the use of the software.
Thepotential of Miller and Spooner’s software for exposing numerical instability isclearly demonstrated by the case studies in these references, yet the softwarehas apparently not been widely used. This is probably largely due to theinability of the software to analyse algorithms expressed in Fortran, or anyother standard language.A different approach to algorithm analysis is taken by Larson and Sameh[688, 1978], [689, 1980], and implemented in software by Larson, Pasternak,and Wisniewski [687, 1983].
Here, errors are measured in a relative ratherthan an absolute sense, and the stability is analysed at fixed data instead ofattempting to maximize instability over all data; however, the analysis is stilllinearized.The idea of applying automatic differentiation to a computational graphto obtain estimates for the forward error in an algorithm is found not only inthe references cited above, but also in the automatic differentiation literature;see Rail [857, 1981], Iri [606, 1991], and Kubota [675, 1991], for example.Hull [589, 1979] discusses the possibility of applying program verificationtechniques to computations that are subject to rounding errors, with the aimof proving a program “correct”. Difficulties include deciding what “correct”means and formulating appropriate assertions.
Although he reports someprogress, Hull concludes that “there is not a great deal to be gained by tryingto apply the techniques of program verification to programs for numericalcalculations”.Bliss, Brunet, and Gallopoulos [126, 1992] develop Fortran preprocessortools for implementing the local relative error approach of Larson and Sameh.Their tools also implement a statistical technique of Brunet and Chatelin[152, 1989], [202, 1990] in which the result of every floating point operationis randomly perturbed and statistical methods are used to measure the effecton the output of the algorithm.Rowan [882, 1990] develops another way to search for numerical instability.
For an algorithm with data d he maximizes S(d) = e(d)/ k (d) using adirect search maximizer he has developed called the subplex method (which24.6 N OTESANDR EFERENCES489is based on the Nelder–Mead simplex method). Here, e(d) = yacc - is anapproximation to the forward error in the computed solution where yacc is amore accurate estimate of the true solution than and the condition numberk(d) is estimated using finite difference approximations. The quantity S(d) isa lower bound on the backward error of the algorithm at d.
Fortran softwaregiven in [882, 1990] implements this “functional stability analysis”. The software takes as input two user-supplied Fortran subprograms; one implementsthe algorithm to be tested in single precision, and the other provides a moreaccurate solution, typically by executing the same algorithm in double precision. The examples in [882, 1990] show that Rowan’s software is capable ofdetecting numerical instability in a wide variety of numerical algorithms.A technique called “significance arithmetic” was studied by Ashenhurst,Metropolis, and others in the 1960s. It involves performing arithmetic onunnormalized numbers in such a way that the number of significant digitsin the answers provides estimates of the accuracy of those answers.
Significance arithmetic is therefore a form of automatic error analysis. For moredetails see, for example, Ashenhurst and Metropolis [31, 1965], [750, 1977]and Sterbenz [938, 1974, §7.2]; Sterbenz explains several drawbacks of thetechnique.Stoutemyer [958, 1977] describes the use of the symbolic manipulationpackage REDUCE to analyse the propagation of data and rounding errors innumerical algorithms.Finally, we note that running error analysis (see §3.3) is a form of automatic error analysis and is an attractive alternative to interval arithmeticas a means of computing a posteriori error bounds for almost any numericalalgorithm.24.6.
Notes and ReferencesSections 24.1 – 24.3.2 and §24.5 are based on Higham [555, 1993].The Release Notes for MATLAB 4.1 state that “This release of MATLABfixes a bug in the rcond function. Previously, rcond returned a larger thanexpected estimate for some matrices . . . rcond now returns an estimate thatmatches the value returned by the Fortran LINPACK library.” In the directsearch experiments in [555, 1993] we used MATLAB 3.5, and we found it mucheasier to generate counterexamples to rcond than we do now with MATLAB4.2. It seems that the maximizations in [555, 1993] were not only defeatingthe algorithm underlying rcond, but also, unbeknown to us, exploiting a bugin the implemental ion of the function.
The conclusions of [555, 1993] areunaffected, however.Another way to solve a cubic is to use Newton’s method to find a realzero, and then to find the other two zeros by solving a quadratic equation.A UTOMATIC E RROR A NALYSIS490In a detailed investigation, Kahan [631, 1986] finds this iterative method tobe preferable to (sophisticated) use of the explicit formulae.
Other usefulreferences on the numerical computation of roots of a cubic are Lanczos [686,1956, Chap. 1], Press, Teukolsky, Vetterling, and Flannery [842, 1992, §5.6],and Uspensky [1038, 1948].Problems24.1. Let Aand letandbe the computed orthogonalQR factors produced by the classical and modified Gram-Schmidt methods,respectively. Use direct search to maximize the functionsIn order to keep k2(A) small, try maximizing fi (A) – θ max( k2(A) – µ, 0),where θ is a large positive constant and µ is an upper bound on the acceptablecondition numbers.24.2.
It is well known that ifthenis nonsingular and vTA–1 u–1This is known as the Sherman-Morrison formula. For a history and generalizations see Henderson and Searle [512, 1981]. A natural question is whetherthis formula provides a stable way to solve a rank-1 perturbed linear system.That is, is the following algorithm stable?% Solve Cx: = (A + uvT)x = b.Solve Ay = b for y.Solve Az = u for z.x = y - (v T y)(1 + v T z) - 1 z(a) Investigate the stability using direct search.
Let both linear systemswith coefficient matrix A be solved by GEPP. Take A, u, and v as the dataand let the function to be maximized be the normwise backward error η C,b inthenorm.(b) (RESEARCH PROBLEM) Obtain backward and forward error bounds forthe method (for some existing analysis see Yip [1118, 1986]).24.3.
(RESEARCH PROBLEM) Investigate the stability of the formulae of $24.3.3for computing the roots of a cubic.Previous Home NextChapter 25Software Issues in Floating PointArithmeticThe first American Venus probe was lost due to aprogram fault caused by theinadvertent substitution of a statement of the formDO 3 I = 1.3 for one of the form DO 3 I = 1,3.— JIM HORNING, Note on Program Reliability 2 2 ( 1 9 7 9 )Numerical subroutines should deliver results that satisfy simple,useful mathematical laws whenever possible.— DONALD E.
KNUTH, The Art of ComputerProgramming, Volume 2, Seminumerical Algorithms (1981)No method of solving a computational problem isreally available to a user until it iscompletely described in an algebraic computing languageand made completely reliable.Before that, there are indeterminate aspects in any algorithm.— GEORGE E. FORSYTHE,Today’s Computational Methods of Linear Algebra (1967)The extended precision calculation of pi has substantial application as atest of the “global integrity” of a supercomputer .
. .Such calculations . . . are apparently now used routinelyto check supercomputers before they leave the factoty.A large-scale calculation of pi is entirely unforgiving;it soaks into all parts of the machine and asingle bit awry leaves detectable consequences.— J. M. BORWEIN, P. B. BORWEIN, and D. H. BAILEY,Ramanujan, Modular Equations, and Approximations to Pior How to Compute One Billion Digits of Pi (1989)22Quoted, with further details, in Tropp [1021, 1984].491S OFTWARE I SSUES492FLOATING P OINT A RITHMETICINIn this chapter we discuss some miscellaneous aspects of floating point arithmetic that have an impact on software development.25.1.
Exploiting IEEE ArithmeticIEEE standard 754 and 854 arithmetic can be exploited in several waysin numerical software, provided that proper access to the arithmetic’s features is available in the programming environment. Unfortunately, althoughmost commercially significant floating point processors at least nearly conformto the IEEE standards, language standards and compilers generally providepoor support (exceptions include the Standard Apple Numerics Environment(SANE) [233, 1988], Apple’s PowerPC numerics environment [234, 1994], andSun’s SPARCstation compilers [969, 1992], [970, 1992]).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
