Numerical Reproducibility in High-Performance Computing
David H. Bailey
http://www.davidhbailey.com
Computational Research Department, Lawrence Berkeley National Laboratory
Computer Science Department, University of California, Davis

17 Jun 2013

Reproducibility in scientific computing
A December 2012 workshop on reproducibility in computing, held at Brown University
in Rhode Island, USA, noted that
Science is built upon the foundations of theory and experiment validated andimproved through open, transparent communication. With the increasingly
central role of computation in scientific discovery this means communicating
all details of the computations needed for others to replicate the experiment.
... The “reproducible research” movement recognizes that traditional
scientific research and publication practices now fall short of this ideal, andencourages all those involved in the production of computational science ...
to facilitate and practice really reproducible research.
V. Stodden, D.H. Bailey, J. Borwein, R.J. LeVeque, W. Rider and W. Stein, “Setting the default to reproducible: Reproducibility in computationaland experimental mathematics,” http://www.davidhbailey.com/dhbpapers/icerm-report.pdf.

Numerical reproducibility

The reproducibility report further noted:
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 beencorrectly ported to a new system, because computational results quickly
diverge from standard benchmark cases. And it is doubly difficult for otherresearchers, using independently written codes and distinct computer
systems, to reproduce published results.
V. Stodden, D.H. Bailey, J. Borwein, R. J. LeVeque, W. Rider and W. Stein, “Setting the default to reproducible: Reproducibility in computationaland experimental mathematics,” http://www.davidhbailey.com/dhbpapers/icerm-report.pdf.

The growing problem of numerical reliability
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.
Particularly vulnerable are:
1. Large-scale, highly parallel simulations, running on systems with hundreds ofthousands or millions of processors — numerical sensitivities are greatly magnified.
2. Certain applications with highly ill-conditioned linear systems.
3. Large summations, especially those involving cancellations.
4. Long-time, iterative simulations (such as molecular dynamics or climate models).
5. Computations to resolve small-scale phenomena.
6. Studies in computational physics or experimental mathematics often require hugeprecision levels.
D.H. Bailey, R. Barrio, and J.M. Borwein, “High precision computation: Mathematical physics and dynamics,” Applied Mathematics andComputation, vol. 218 (2012), pg. 10106–10121.

Analysis of collisions at the Large Hadron Collider
I The 2012 discovery of the Higgs boson at the ATLAS experiment in the LHC

relied crucially on the ability to track charged particles with exquisite precision (10microns over a 10m length) and high reliability (over 99% of roughly 1000
charged particles per collision correctly identified).
I Software: 5 millions line of C++ and python code, developed by roughly 2000

physicists and engineers over 15 years.
I Recently, in an attempt to speed up the calculation, researchers found that merely

changing the underlying math library resulted in some collisions being missed or
misidentified.
Questions:
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?

Ways to enhance numerical reproducibility and solve accuracy problems
1. 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.
2. Carefully analyze every step of the computation to ascertain the level of numericalsensitivity at each step.
3. Employ interval arithmetic for large portions of the application (which greatly
increases run time and code complexity).
4. Employ higher-precision arithmetic arithmetic, assisted with some “smart” tools
to help determine where extra precision is needed and where it is not.
Item #4 is the only practical solution.

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 had rigorous training in numerical analysis!

Enhancing reproducibility with high-precision arithmetic
Problem: Findof the irregular function
P the arc length−kg (x) = x + 0≤k≤10 2 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.0

0.5

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

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

0.5

1.0

1.5

2.0

2.5

3.0

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

Aren’t 64 bits enough?
There has been considerable resistance in the scientific computing community to thenotion that more than 64-bit arithmetic is not only useful, but may even be essential insome scientific computations. Why?
I Many are persuaded that physical reality fundamentally does not require highprecision, beyond, say, the limits of 64-bit IEEE arithmetic.
I Many believe that numerical sensitivity problems are always due to the usage ofinferior algorithms.
I Some researchers regard the usage of high-precision arithmetic as “cheating” or“sinful” — because it is too easy?
I Until a few years ago, easy-to-use high-precision arithmetic software facilities were

not widely available to experiment with, or the computational cost was too great.
D.H. Bailey, R. Barrio, and J.M. Borwein, “High precision computation: Mathematical physics and dynamics,” Applied Mathematics andComputation, vol. 218 (2012), pg. 10106–10121.

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..
P .

nn
k=1 xkPn
Pnk=1 xk2
k=1 xk..
P . n+1
k=1 xk

···
···
..
.
···

Pn

  Pnna0k
Pn k=1 xn+1
Pn k=1 yk  a1  
k=1 xkk=1 xk yk    ..  = ..
..
 .  

.
.
P
Pn
n
n
2n
a
x
y
x
n
k=1 k kk=1 k

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 8

Algorithm changes versus double-double?
Double-double is the pragmatic choice
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
In contrast, a straightforward Linpack scheme, implemented with double-double
arithmetic, works fine for this and a wide range of similar problems.

Free software for high-precision computation
1. ARPREC. Arbitrary precision, with numerous algebraic and transcendental
functions. High-level interfaces for C++ and Fortran-90.
http://crd.lbl.gov/~dhbailey/mpdist.
2. GMP. Produced by a volunteer effort and distributed under the GNU license.
http://gmplib.org.
3. MPFR. C library for multiple-precision floating-point computations with exact
rounding, based on GMP. http://www.mpfr.org.
4. MPFR++. High-level C++ interface to MPFR.
http://perso.ens-lyon.fr/nathalie.revol/software.html.
5. MPFUN90. Similar to ARPREC, but is written entirely in Fortran-90 and provides
only a Fortran-90 interface. http://crd.lbl.gov/~dhbailey/mpdist.
6. QD. Performs “double-double” (31 digits) and “quad-double” (62 digits)
arithmetic. High-level interfaces for C++ and Fortran-90.
http://crd.lbl.gov/~dhbailey/mpdist.

U.C. Berkeley’s CORVETTE project and the “Precimonious” toolObjective: Develop software facilities to find and ameliorate numerical anomalies inlarge-scale computations:
I Facilities to test the level of numerical accuracy required for an application.
I Facilities to delimit the portions of code that are inaccurate.
I Facilities to search the space of possible code modifications.
I Facilities to repair numerical difficulties, including usage of high-precision

arithmetic.
I Facilities to navigate through a hierarchy of precision levels (32-bit, 64-bit, 80-bit

or higher as needed).
The current version of this tool is known as “Precimonious.”
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,” manuscript, May 2013.

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.

Some applications where high precision is useful or essential
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).
D.H. Bailey, R. Barrio, and J.M. Borwein, “High precision computation: Mathematical physics and dynamics,” Applied Mathematics andComputation, vol. 218 (2012), pg. 10106–10121.

Supernova simulations
I Researchers at LBNL have used quad-double arithmetic

to solve for non-local thermodynamic equilibrium
populations of iron and other atoms in the atmospheresof supernovas.
I Iron may exist in several species, so it is necessary tosolve for all species simultaneously.
I Since the relative population of any state from thedominant state is proportional to the exponential of theionization energy, the dynamic range of these values canbe very large, and cancellations may occur.
I Quad-double arithmetic (62 digits) was employed toreduce these errors.

P.H. Hauschildt and E. Baron, “Thenumerical solution of the expanding stellar
atmosphere problem,” Journal
Computational and Applied Mathematics,
vol. 109 (1999), pg. 41–63.

Climate modeling: High-precision for reproducibility

I Climate and weather simulations are fundamentallychaotic: if microscopic changes are made to the current
state, soon the future state is quite different.
I In practice, computational results are altered even if

minor changes are made to the code or the system.
I This numerical variation is a major nuisance for code

maintenance.
I He and Ding found that by using double-double

arithmetic in two key inner loops, most of thisnumerical variation disappeared.

Y. He and C. Ding, “Using accurate
arithmetics to improve numericalreproducibility and stability in parallel
applications,” Journal of Supercomputing,vol. 18, no. 3 (Mar 2001), pg. 259–277.

Experimental mathematics:
Discovering new mathematical results by computer
Methodology:
1. Compute various mathematical entities (limits, infinite series sums, definite
integrals, etc.) to high precision, typically 100–10,000 digits.
2. Use algorithms such as PSLQ to recognize these numerical values in terms ofwell-known mathematical constants.
3. When results are found experimentally, seek formal mathematical proofs of the
discovered relations.
Many results have recently been found using this methodology, both in puremathematics and in mathematical physics.
“If mathematics describes an objective world just like physics, there is no reason why
inductive methods should not be applied in mathematics just the same as inphysics.” – Kurt GodelIsing integrals from mathematical physics
We recently applied our methods to study three classes of integrals (one of which wasreferred to us by Craig Tracy of U.C. Davis) that arise in the Ising theory ofmathematical physics:
Z
Z ∞
1
du1
dun
4 ∞
···
···
Cn :=
P
2
n! 0
u1
un
n
0
j=1 (uj + 1/uj )Q  ui −uj 2
Z ∞Z ∞i<j ui +uj
4
dun
du1
Dn :=
···
···
P
2
n! 0
u1
un
n
0
j=1 (uj + 1/uj )

2
Z 1
Z 1
Y uk − uj dt2 dt3 · · · dtn
En = 2···
uk + uj
0
0
1≤j<k≤nwhere in the last line uk = t1 t2 · · · tk .
D.H. Bailey, J.M. Borwein and R.E. Crandall, “Integrals of the Ising class,” Journal of Physics A: Mathematical and General, vol. 39 (2006), pg.
12271–12302.

Limiting value of Cn : What is this number?
Key observation: The Cn integrals can be converted to one- dimensional integrals
involving the modified Bessel function K0 (t):Z
2n ∞ nCn =tK0 (t) dt
n! 0
High-precision numerical values, computed using this formula and tanh-sinhquadrature, approach a limit. For example:
C1024 = 0.6304735033743867961220401927108789043545870787 . . .
What is this number? We copied the first 50 digits into the online Inverse Symbolic
Calculator (ISC) at http://carma-lx1.newcastle.edu.au:8087. The result was:lim Cn = 2e −2γ .

n→∞where γ denotes Euler’s constant. This is now proven.

Other Ising integral evaluations found using PSLQ
D2 = 1/3
D3 = 8 + 4π 2 /3 − 27 Li−3 (2)D4 = 4π 2 /9 − 1/6 − 7ζ(3)/2E2 = 6 − 8 log 2E3 = 10 − 2π 2 − 8 log 2 + 32 log2 2
E4 = 22 − 82ζ(3) − 24 log 2 + 176 log2 2 − 256(log3 2)/3+16π 2 log 2 − 22π 2 /3?

E5 = 42 − 1984 Li4 (1/2) + 189π 4 /10 − 74ζ(3) − 1272ζ(3) log 2+40π 2 log2 2 − 62π 2 /3 + 40(π 2 log 2)/3 + 88 log4 2+464 log2 2 − 40 log 2
where ζ(x) is the Riemann zeta function and Lin (x) is the polylogarithm function.

Algebraic numbers in Poisson potential functions associated with latticesums
Lattice sums arising from the Poisson equation have been studied widely inmathematical physics and also in image processing. We numerically discovered, andthen proved, that for rational (x, y ), the two-dimensional Poisson potential function
satisfies
1 X cos(mπx) cos(nπy )1φ2 (x, y ) =
=
log α
π2
m2 + n2πm,n odd

where α is an algebraic number, i.e., the root of an integer polynomial
0 = a0 + a1 α + a2 α 2 + · · · + an α n
The minimal polynomials for these α were found by PSLQ calculations, with the(n + 1)-long vector (1, α, α2 , · · · , αn ) as input, where α = exp(8πφ2 (x, y )). PSLQ
returns the vector of integer coefficients (a0 , a1 , a2 , · · · , an ) as output.
1. D.H. Bailey, J.M. Borwein, R.E. Crandall and J. Zucker, “Lattice sums arising from the Poisson equation,” Journal of Physics A:
Mathematical and Theoretical, vol. 46 (2013), pg. 115201, http://www.davidhbailey.com/dhbpapers/PoissonLattice.pdf.
2. D.H. Bailey and J.M. Borwein, “Compressed lattice sums arising from the Poisson equation: Dedicated to Professor Hari Sirvastava,”
Boundary Value Problems, vol. 75 (2013), DOI: 10.1186/1687-2770-2013-75,
http://www.boundaryvalueproblems.com/content/2013/1/75.

Samples of minimal polynomials found by PSLQ
k
5
6
7
8
9

10

Minimal polynomial for exp(8πφ2 (1/k, 1/k))
1 + 52α − 26α2 − 12α3 + α4
1 − 28α + 6α2 − 28α3 + α4
−1 − 196α + 1302α2 − 14756α3 + 15673α4 + 42168α5 − 111916α6 + 82264α7−35231α8 + 19852α9 − 2954α10 − 308α11 + 7α12
1 − 88α + 92α2 − 872α3 + 1990α4 − 872α5 + 92α6 − 88α7 + α8
−1 − 534α + 10923α2 − 342864α3 + 2304684α4 − 7820712α5 + 13729068α6−22321584α7 + 39775986α8 − 44431044α9 + 19899882α10 + 3546576α11−8458020α12 + 4009176α13 − 273348α14 + 121392α15−11385α16 − 342α17 + 3α18
1 − 216α + 860α2 − 744α3 + 454α4 − 744α5 + 860α6 − 216α7 + α8

The minimal polynomial for exp(8πφ2 (1/32, 1/32)) has degree 128, with individual
coefficients ranging from 1 to over 1056 . This PSLQ computation required 10,000-digit
precision. See next page.
Other polynomials required up to 50,000-digit precision.

Summary
I Reproducibility in general, and numerical reproducibility in particular, is looming

as a major challenge for future exascale computations.
I The most practical solution to such difficulties is to judiciously employ

higher-precision arithmetic, combined with some “smart” tools to identify
sensitivities and make the requisite code modifications.
I In most cases, only a handful of particularly troublesome loops need to becomputed using higher precision — thus, the total run time is only slightly
increased.
I While 32-digit or 64-digit arithmetic suffices for almost all conventional scientificcomputations, a growing body of scientifically interesting computations require
hundreds or even thousands of digits.
I Fortunately, such computations are now feasible using modern high-performancecomputing technology, together with relatively easy-to-use software facilities.
This talk is available at http://www.davidhbailey.com/dhbtalks/dhb-reproducibility.pdf.