Numerical reproducibility in high-performance computing
David H. Bailey
http://www.davidhbailey.com
Lawrence Berkeley National Laboratory (retired)Computer Science Department, University of California, Davis

November 19, 2015

1 / 24

Big data science: DANGER AHEAD

Supercomputers operating on big data cangenerate nonsense faster than ever before!

Key concerns:
I Are the results statistically sound?
I Are the results numerically reliable?
I Have the results been validated using independent rigorous tests?
I Are the algorithms, data sources and processing methods well documented?

2 / 24

Reproducibility crises in physics and cosmology
I In 2011, an international team of researchers at the Gran Sasso Laboratory in Italy

announced that neutrinos had exceeded the speed of light, thus directly
challenging Einstein’s relativity. However, after months of careful checking, asubtle flaw was found in the measurement apparatus (a major embarrassment).
I In 2013, CERN researchers confirmed the discovery of the long-sought Higgs

boson. But more recently, scientists have raised questions as to whether theparticle discovered is really the Higgs – it might be some other particle or particles
masquerading as the Higgs; additional research studies are required.
I In March 2014, researchers announced with considerable fanfare that they had

detected the fingerprint of the long-hypothesized inflationary epoch, a tiny
fraction after the big bang. Sadly, within a few weeks critics pointed out that
their experimental results might well be due to dust in the Milky Way, pending
better data.
3 / 24

Reproducibility crises in biomedicine, psychology, economics, finance, autos
I In 2011, Bayer researchers reported that they wereable to reproduce only 17 of 67 pharma studies.
I In 2012, Amgen researchers reported that they wereable to reproduce only 6 of 53 cancer studies.
I In August 2015, the Reproduciblity Project inVirginia reported that they were able to reproduce
only 39 of 100 psychology studies.
I In September 2015, the U.S. Federal Reserve wasable to reproduce only 29 of 67 economics studies.
I In 2014-2015, “backtest overfitting” emerged as a

Reproducibility Project staff
Credit: NY Times

major problem in computational finance.
I In March 2014, West Virginia researchers reported that they were unable toreproduce Volkswagen’s claimed emission figures. This has now exploded into amajor international scandal.
4 / 24

ICERM workshop report on numerical reproducibility
Numerical round-off error and numerical differences are greatly magnified ascomputational simulations are scaled up to run on highly parallel systems. As a result,
it is increasingly difficult to determine whether a code has been correctly ported to a
new system, because computational results quickly diverge from standard benchmark
cases. And it is doubly difficult for other researchers, using independently writtencodes and distinct computer systems, to reproduce published results.

I V. Stodden, D. H. Bailey, J. M. Borwein, R. J. LeVeque, W. Rider and W. Stein, “Setting thedefault to reproducible: Reproducibility in computational and experimental mathematics,” 2012,
http://www.davidhbailey.com/dhbpapers/icerm-report.pdf.

5 / 24

Numerical reproducibility
Many applications routinely use either 32-bit or 64-bit IEEE arithmetic, and employ
fairly simple algorithms, assuming that all is well. But problems can arise:
I Highly ill-conditioned linear systems.
I Large summations, especially those involving +/- terms and cancellations.
I Long-time, iterative simulations (such as molecular dynamics or climate models).
I Computations to resolve small-scale phenomena.
I Studies in computational physics or experimental mathematics often require hugeprecision levels to produce reliable results.
Large-scale, highly parallel computations greatly magnify numerical sensitivities.

6 / 24

Analysis of collisions at the Large Hadron Collider
The 2012 discovery of the Higgs boson at the ATLAS experiment in the LHC reliedcrucially on the ability to track charged particles with exquisite precision (10 microns
over a 10m length) and high reliability (over 99% of roughly 1000 charged particles per
collision correctly identified).
Recently, in an attempt to speed up the calculation, researchers found that merely
changing the underlying math library resulted in some collisions being missed ormisidentified.
I How serious are these numerical difficulties?
I How can they be tracked down?
I How can the library be maintained, producing numerically reliable results?

7 / 24

Approaches to solve accuracy and numerical reproducibility problems
I Employ an expert numerical analyst to re-examine every algorithm employed in thecomputation to ensure that only the most stable known schemes are being used.
I Carefully analyze every step of the computation to ascertain the level of numericalsensitivity at each step.
I Employ interval arithmetic for large portions of the application (which greatlyincreases run time and code complexity).
I Employ some “smart” tools to help identify the numerically sensitive parts of anapplication, then use higher-precision arithmetic arithmetic if needed to fix them.
The last item is the only practical, time-critical solution.

8 / 24

Numerical analysis expertise among U.C. Berkeley graduates
Of the 2010 U.C. Berkeley graduating class, 870 were in disciplines likely to require
technical computing:
I Division of Mathematical and Physical Sciences (Math, Physics, Statistics).
I College of Chemistry.
I College of Engineering (including Computer Science).

Other fields (not counted) that will likely involve significant computing:
I Biology, geology, medicine, economics, psychology, sociology.

Enrollment in numerical analysis courses:
I Math 128A (Introductory numerical analysis required of math majors): 219.
I Math 128B (A more advanced course, required to do serious work): 24.

Conclusion: Fewer than 2% of Berkeley graduates who will do technical computing
have rigorous training in numerical analysis.
9 / 24

Selective usage of high-precision arithmetic
High-precision arithmetic is an attractive solution for many present-day reproducibility
problems, because it can be implemented on an application with only minor changes tosource code.
Free software packages for high-precision computation
1. QD. Performs “double-double” (31 digits) and “quad-double” (62 digits)
arithmetic. High-level interfaces for C++ and Fortran-90.
http://crd-legacy.lbl.gov/˜dhbailey/mpdist.
2. MPFR. C library for arbitrary-precision floating-point computations with exact
rounding, based on GMP. http://www.mpfr.org.
3. MPFR++. High-level C++ interface to MPFR.
http://www.holoborodko.com/pavel/mpfr/.
4. MPFUN2015. High-level Fortran interface to MPFR.
http://crd-legacy.lbl.gov/˜dhbailey/mpdist.
High-precision does not guarantee bit-for-bit reproducibility, but it can ensure
reproducibility (to a given number of digits) with high probability.
10 / 24

Enhancing reproducibility by selective use of double-double arithmetic
Problem: Find the arc length of the irregular function
Pg(x ) = x + 0≤k≤10 2−k sin(2k x ), over the interval (0, π)(using 107 abscissa points).
3.0

I If this computation is done with ordinary double

precision arithmetic, the calculation takes 2.59 seconds
and yields the result 7.073157029008510.
I If it is done using all double-double arithmetic (31-digit

2.5

2.0

1.5

1.00.5accuracy), it takes 47.39 seconds seconds and yields theresult 7.073157029007832 (correct to 15 digits).
I But if only the summation is changed to double-double,

the result is identical to the double-double result (to 15digits), yet the computation only takes 3.47 seconds.

0.5

1.0

1.5

2.0

2.5

3.0

Graph of g(x ) =Px + 0≤k≤10 2−k sin(2k x ),over (0, π).

11 / 24

U.C. Berkeley’s “Precimonious” toolObjective: Develop software facilities to find and ameliorate numerical anomalies inlarge-scale computations:
I Test the level of numerical accuracy required for an application.
I Delimit the portions of code that are inaccurate.
I Search the space of possible code modifications.
I Repair numerical difficulties, including usage of high-precision arithmetic.
I Navigate through a hierarchy of precision levels (32-bit, 64-bit, 80-bit or higher asneeded).
The current version of this tool is known as “Precimonious.”

I C. Rubio-Gonzalez, C. Nguyen, H. D. Nguyen, J. Demmel, W. Kahan, K. Sen, D. H. Bailey and

C. Iancu, “Precimonious: Tuning assistant for floating-point precision,” Proceedings of SC13,
Nov 2013.
12 / 24

Precimonious in action for the arc length computation
Original code:

Automatically tuned code:

int main() {int i, n = 1000000;
long double h, t1, t2, dppi;
long double s1;...
for (i = 1; i <= n; i++) {t2 = fun (i * h);s1 = s1 + sqrt (h*h + (t1-t2)*(t1-t2));
t1 = t2;}
return 0;}

int main() {int i, n = 1000000;
double h, t1, t2, dppi;
long double s1;...
for (i = 1; i <= n; i++) {t2 = fun (i * h);s1 = s1 + sqrt (h*h + (t1-t2)*(t1-t2));
t1 = t2;}
return 0;}

The tuned code produces the same answers, to 15-digit accuracy, yet runs faster.

13 / 24

The ReproBLAS project
Researchers who have supported the basic linear algebra subroutine (BLAS) library aredeveloping versions that guarantee reproducibility (to the last bit) for results. Some ofthese routines are being proposed to be incorporated into the message passing
interface (MPI) standard.
Example: A summation routine returns the sum of a vector of 64-bit summands, withthe same result regardless of the order of the summands.
For details, see http://bebop.cs.berkeley.edu/reproblas

14 / 24

ReproBLAS summation on 1024-processor Cray XC30 (1M summands)

15 / 24

Bit-for-bit reproducibility: pluses and minuses
Full-application bit-for-bit reproducibilit is nearly feasible at the present time.
Advantages:
I Some applications have legal requirements for bit-for-bit reproducibility.
I Aids parallel debugging and porting to other systems.

Disadvantages:
I Runs several times longer.
I Potentially conflicts with advances in compiler optimization.
I Masks a host of potential numerical accuracy problems, since one never sees theeffect of slight perturbations.
Would interval arithmetic be a better solution for those with legal requirements?
We need to discuss this more in the high-performance computing community.
16 / 24

John Gustafson’s new “unum” system for floating-point
Gustafson proposes a radical new definition of floating-point arithmetic, based onunums (short for “universal numbers”), a flexible format that contains a “ubit” todesignate that the true number lies between two adjacent binary floats.
Example: Avogadro’s constant = 6.022 × 1023 is (29 bits total):
fracubit exp size frac size
sign exp
11001101 111111100001 11111011
0Advantages:
I Obeys commutative law, transitive law and other laws of algebra.
I No rounding error, overflow, underflow, negative zero.
I Consistent treatment of positive and negative infinity and NaN.
I Safe to parallelize; portable between systems.
I Often requires fewer bits than conventional IEEE 64-bit arithmetic.

17 / 24

Sample problem (due to Jean-Michel Muller)
Define:E (z) = (e z − 1)/z,Q(x ) = |x −p

E (0) = 1,x 2 + 1| − 1/(x +px 2 + 1),

H(x ) = E (Q 2 (x )).
Problem: Compute H(x ) for x = 15, 16, 17, 9999.
I IEEE 32-bit: (0, 0, 0, 0).
I IEEE 64-bit: (0, 0, 0, 0).
I IEEE 128-bit: (0, 0, 0, 0).
I Unum (correct): (1, 1, 1, 1).

18 / 24

Another sample problem (due to Siegfried Rump)
DefineF (x , y ) = (333 + 3/4)y 6 + x 2 (11x 2 y 2 − y 6 − 121y 4 − 2) + (5 + 1/2)y 8 + x /(2y )Problem: Compute F (77617, 33096).
I IEEE 32-bit: 1.172603.
I IEEE 64-bit: 1.1726039400531.
I IEEE 128-bit: 1.172603940053178.
I Unum (correct to 23 digits): −0.82739605994682136814116 . . ..

19 / 24

Unum: Current status and plans
I Gustafson’s book has been published by CRC Press (see below).
I A full implementation is nearly complete in C.
I Much testing and analysis will be done.
I Will require a flexible working precision level.
I Hardware implementation?

For additional details:
I John Gustafson, The End of Error: Unum Computing, CRC Press, 2015.

20 / 24

Aren’t 64 bits enough?
Some applications actually require more than 64-bit arithmetic:
1. Planetary orbit calculations (32 digits).
2. Supernova simulations (32–64 digits).
3. Climate modeling (32 digits).
4. Coulomb n-body atomic system simulations (32–120 digits).
5. Schrodinger solutions for lithium and helium atoms (32 digits).
6. Electromagnetic scattering theory (32–100 digits).
7. Scattering amplitudes of fundamental particles (32 digits).
8. Discrete dynamical systems (32 digits).
9. Theory of nonlinear oscillators (64 digits).
10. The Taylor algorithm for ODEs (100–600 digits).
11. Ising integrals from mathematical physics (100–1000 digits).
12. Problems in experimental mathematics (100–50,000 digits).
I D. H. Bailey, R. Barrio, and J. M. Borwein, “High precision computation: Mathematical physics

and dynamics,” Applied Mathematics and Computation, vol. 218 (2012), pg. 10106–10121.
21 / 24

Innocuous example where standard 64-bit precision is inadequate
Problem: Find a polynomial to fit the data (1, 1048579, 16777489,84941299, 268501249, 655751251, 1360635409, 2523398179, 4311748609) forarguments 0, 1, · · · , 8. The usual approach is to solve the linear system:





Pn n

k=1 xk..
Pn .

nk=1 xk

Pn
Pnk=1 xk2
k=1 xk..
P . n+1
k=1 xk······..
.
···Pn  Pn
na0kPn k=1 xn+1Pn k=1 yk  a1  
k=1 xkk=1 xk yk  =

..
.
..
  ..  .
.
PnPnn2nan
k=1 xk yk
k=1 xk

A 64-bit computation (e.g., using Matlab, Linpack or LAPACK) fails to find thecorrect polynomial in this instance, even if one rounds results to nearest integer.
However, if Linpack routines are converted to use double-double arithmetic (31-digit
accuracy), the above computation quickly produces the correct polynomial:f (x ) = 1 + 1048577x 4 + x 8 = 1 + (220 + 1)x 4 + x 822 / 24

Innocuous example, cont.
The result on the previous page can be obtained with double precision using Lagrange
interpolation or the Demmel-Koev algorithm. But few scientists, outside of expert
numerical analysts, are aware of these schemes.
Besides, even these schemes fail for higher-degree problems.
For example: (1, 134217731, 8589938753, 97845255883, 549772595201,
2097396156251, 6264239146561, 15804422886323, 35253091827713,
71611233653971, 135217729000001, 240913322581691, 409688091758593)
is generated by:f (x ) = 1 + 134217729x 6 + x 12 = 1 + (227 + 1)x 6 + x 12
Lagrange interpolation and Demmel-Koev fail for this problem, but a straightforward
Linpack scheme, implemented with double-double arithmetic, works fine.
23 / 24

Summary
I Numerical reproducibility is looming as a major challenge for large-scale scientificcomputations.
I One practical solution is to apply intellgient tools that judiciously employ

higher-precision arithmetic for numerically sensitive portions of the code.
I Numerous other tools are being developed to enhance reproducibility or to ensure

bit-for-bit reproducibility.
I Do we want bit-for-bit reproducibility? Will this only mask serious numericalproblems? This issue needs to be further discussed in the community.

This talk is available at
http://www.davidhbailey.com/dhbtalks/dhb-num-repro.pdf.
24 / 24