Reproducibility in scientific computing
David H. Bailey
http://www.davidhbailey.com
Lawrence Berkeley National Laboratory, USA (retired)Computer Science Department, University of California, Davis, USA
Collaborators:
Jonathan M. Borwein (CARMA, University of Newcastle, Australia)
Victoria Stodden (University of Illinois Urbana-Champaign, USA)

September 29, 2014

Reproducibility in physics
In March 2014, a team of researchers from Harvard
University made the dramatic announcement thatthey had discovered an interesting “twisting”
pattern in cosmic microwave background data,measured using their BICEP2 experimental system.
This pattern fit very well with the hypothesizedpattern of the most commonly assumed model of“inflation” in the first tiny fraction of a second after
the big bang, and thus has been trumpeted as thefirst experimental evidence of the inflationary
cosmology.
But other researchers had difficulty reconstructing
the claimed results. Finally, two teams challenged
the BICEP2 findings, saying that the results could
more readily be explained by dust in the Milky Way.

Ron Cowen, “Doubt grows about gravitational
waves detection,” Scientific American, 2 Jun 2014.

Reproducibility in biomedicine
The biomedical field has been stung by numerous cases where pharma products lookgood based on clinical trials, but later disappoint in real-world usage, or the resultscannot be reproduced in separate studies. Examples:
I In 2004, GlaxoSmithKline acknowledged that while some trials of Paxil found it

effective for depression in children, other unpublished studies showed no benefit.
I In 2011, Bayer researchers reported that they were able to reproduce the results of

only 17 of 67 published studies they examined.
I In 2012, Amgen researchers reported that they were able to reproduce the resultsof only 6 of 53 published cancer studies.
I In 2014, a review of Tamiflu found that while it made flu symptoms disappear a

bit sooner, it did not stop serious complications or keep people out of the hospital.
These experiences have exposed a fundamental flaw in methodology:
Only publicizing the results of successful trials introduces a bias into the results.
The AllTrials movement would require all results to be public: http://www.alltrials.net

High-performance scientific computing: An amazing success story
Performance of theworld’s top 500
supercomputers(1994 – present)
Red = #1
Orange = #500
Blue = Sum #1 thru #500

I “Performance development,” Top500.org, available at http://top500.org/statistics/perfdevel.

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.

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

Reproducibility in scientific computing: Four key issues

From the ICERM workshop and other recent discussions, four key issues have emerged:
1. The need to institute a culture of reproducibility.
2. The danger of statistical overfitting and other errors in data analysis.
3. The need for greater rigor and forthrightness in performance reporting.
4. The growing concern over numerical reproducibility.
We will deal with each of these items in the remainder of the talk.

Reproducibility in scientific computing:
1. A culture of reproducibility
The ICERM report noted:
Unfortunately the scientific culture surrounding computational work has
evolved in ways that often make it difficult to verify findings, efficiently build
on past research, or even to apply the basic tenets of the scientific method tocomputational procedures. Bench scientists are taught to keep careful lab
notebooks documenting all aspects of the materials and methods they use
including their negative as well as positive results, but computational work isoften done in a much less careful, transparent, or well-documented manner.
Often there is no record of the workflow process or the code actually used toobtain the published results, let alone a record of the false starts. Thisultimately has a detrimental effect on researchers’ own productivity, theirability to build on past results or participate in community efforts, and thecredibility of the research among other scientists and the public.

Some commonly heard opinions and practices

I “We can’t provide completely full details of our computational techniques —these are our competitive edge over other research teams...”I “Source code? A student wrote that. It’s not very important...”I “What was the system environment? We just used what the campus computer

center had available at the time...”I Which version of the source code did we use? I don’t know. After running thecode for our conference submission, we immediately started making some
additional enhancements, so that is not available...”Specific cultural issues mentioned in the ICERM report
I The need to carefully document the full context of computational

experiments—system environment, input data, code used, computed results, etc.
I The need to grant open access, whenever possible, to computer code and input

datasets, permitting other researchers to rerun the calculation and reproduce thepublished results.
I The need to save the code and relevant data in a permanent repository.
I The need for journal editors to acknowledge importance of computational details,
and to encourage full documentation of the computational techniques, perhaps on
an auxiliary website if not in the journal itself.
I The need for reviewers, research institutions and funding agencies to recognize theimportance of computing and computing professionals, and to allocate funding forafter-the-grant support and repositories.
I The need to encourage publication of negative results—other researchers canoften learn from them.

Best practices for publishing research in scientific computing
1. A precise statement of assertions to be made in the paper.
2. The computational approach, and why it constitutes a rigorous test.
3. Complete statements of, or references to, every algorithm employed.
4. Auxiliary software (both research and commercial software).
5. Test environment, including hardware, software and number of processors.
6. Data reduction and statistical analysis methods.
7. Adequacy of precision level and grid resolution.
8. Full statement or summary of experimental results.
9. Verification and validation tests performed.
10. Availability of computer code, input data and output data.
11. Curation: where are code and data available?
12. Instructions for repeating computational experiments.
13. Terms of use and licensing. Ideally code and data “default to open.”
14. Avenues explored and negative findings.
15. Proper citation of all code and data used.

Consequences of high performance
Very high performance computer systems, either massively parallel networks or
stand-alone workstations, provide unprecedented capability to process huge datasets
and simulate complex natural phenomena.
But there are dangers ahead:
1. Since computer programs can easily explore many variations of a model, given an
input dataset, the resulting “optimal” models are often statistically overfit.
2. Performance analyses of these complex systems are fraught with potential error
and distortion.
3. Large computations greatly magnify sensitivity to numerical roundoff error,
yielding results that are in some cases unreliable and nonreproducible.
We will explore each of these dangers in this talk.

Reproducibility in scientific computing:
2. Statistical overfitting and data analysis errors
I Statistical overfitting means either the development of a model that isfundamentally more complex than the data that it is based on, or else presenting
a single “optimal” model, based on an input dataset, not mentioning many othervariations of the model that were explored and discarded.
I Overfitting is increasingly common in scientific computing, because the rapidly

increasing performance of modern scientific computing systems makes it very easy
to explore thousands, millions or even billions of variations of a model.
I Overfitting can be seen as an instance of selection bias, namely the well-known

effect where a researcher only discloses the results of successful experiments, andignores or hides results of unsuccessful experiments.

I David J. Hand, The Improbability Principle, Macmillan, 2014.

Selection bias
Consider the following absurd financial scheme:
I A financial advisor sends letters to 10, 240 = 10 × 210 potential clients, with 5120

letters predicting a certain security will go up, and the other half predicting it will
go down.
I One month later, the advisor sends letters only to the 5120 investors who werepreviously sent the correct prediction, with 2560 letters predicting a certain
security will go up, and the other half predicting it will go down.
I The advisor continues this process for 10 months.
I The remaining ten investors, having seen ten consecutive spot-on predictions,doubtless would be impressed at the advisor’s remarkable prescience!
This string of ten predictions is a clear instance of selection bias, since the final ten
investors are never told of the many other failed predictions.

Statistical overfitting in mathematical finance
Investment strategies derived from computer simulation often look great on paper, butfall flat in practice, due to backtest overfitting:
I Since present-day high-performance computer systems (or even a workstation) can

easily analyze millions of variations of a given investment strategy, each evaluatedagainst a single historical market dataset (a “backtest”), it is quite likely that the“optimal” strategy will be overfit.
I Using this “optimal” strategy on other datasets, or in actual market trading, will

not reproduce the strategy’s claimed performance — it was “optimized” against
idiosyncrasies of the input dataset, and is otherwise ineffective.
I Backtest overfitting is widespread in the field, since published schemes are

typically not accompanied by any information regarding testing methodology or
how many variations of the strategy were explored when developing the strategy.

How easy is it to overfit a strategy derived from backtests? Very!
I If only 2 years of daily market data are

available, then no more than 7 strategy
variations should be tried.
I If only 5 years of daily market data are

available, then no more than 45 strategy
variations should be tried.

I D. H. Bailey, J. M. Borwein, M. Lopez de Prado and Q. J. Zhu, “Pseudo-mathematics and financial charlatanism: The
effects of backtest over fitting on out-of-sample performance,” Notices of the American Mathematical Society, May
2014, pg. 458–471.

I D. H. Bailey, J. M. Borwein, M. Lopez de Prado and Q. J. Zhu, “The probability of backtest overfitting,” manuscript, 12
Feb 2014, http://ssrn.com/abstract=2326253.

Demonstration of backtest overfitting
The following three viewgraphs present the results of different steps in an attempt tofind an “optimal” investment strategy using a computational search process, based ona simulated (pseudorandom) backtest market dataset.
As you can see, by tweaking some very basic parameters (entry date, sell date, lengthof time in market, stop-loss, etc. — there are many thousands of such combinations),
the computer search process can fit and “predict” the underlying dataset quite well,and yet not possess any fundamental “skill” for any other dataset.
You can try this simulator yourself at: http://datagrid.lbl.gov/backtest

I D. H. Bailey, S. Ger, M. Lopez del Prado, A. Sim and K. Wu, “Statistical overfitting and backtest performance,”
manuscript, 23 Sep 2014, http://www.davidhbailey.com/dhbpapers/overfitting.pdf.

Initial strategy applied to input (In Sample) dataset

Final optimized strategy applied to input (In Sample) dataset

Final optimized strategy applied to new (Out of Sample) dataset

Analysis of the optimization results
I Note that the Sharpe ratio (“SR,” a standard measure of investment

performance) increases from 0.06 for the initial strategy to 1.32 (fairly good) for
the final optimized strategy.
I But when the final “optimized” strategy is applied to the new (Out of Sample)

dataset, the Sharpe ratio is -0.03, indicating a completely ineffective strategy.
I Any attempt to reproduce the promising-looking results obtained from the firstdataset on any other dataset (or in real-life finance) will be disappointed.

Backtest overfitting is now thought to be a chief reason that many new financial
products and strategies look great on paper (based on backtest evaluations) but thenfall flat when actually fielded in the market.

Reproducibility in scientific computing:
3. Performance reporting practices
High-performance computing, 1990-1994:
I Many new parallel systems had been introduced; each claimed theirs was best.
I Many researchers were excited about the potential of highly parallel systems.
I Few standard benchmarks and testing methodologies had been established.
I It was hard to reproduce published performance results; much confusion reigned.
I Overall, the level of rigor and peer review in the field was rather low.

In response, in 1991 DHB published a humorous essay ”Twelve Ways to Fool the
Masses,” poking fun at some of the abuses. Since abuses continued, DHB presented atalk at Supercomputing 1992 and published a paper with specific examples.
I D. H. Bailey, “Misleading performance reporting in the supercomputing field,” Scientific Programming, vol. 1., no. 2(Winter 1992), pg. 141–151.

“Twelve ways to fool the masses in highly parallel computing”
1. Quote 32-bit performance results, not 64-bit results, but don’t mention this in paper.
2. Present performance figures for an inner kernel, then represent these figures as theperformance of the entire application.
3. Quietly employ assembly code and other low-level language constructs.
4. Scale up the problem size with the number of processors, but omit any mention of this.
5. Quote performance results projected to a full system.
6. Compare your results against scalar, unoptimized code on conventional systems.
7. When run times are compared, compare with an old code on an obsolete system.
8. Base Mflop/s rates on the operation count of the parallel implementation, instead of thebest practical serial algorithm.
9. Quote performance as processor utilization, parallel speedups or Mflop/s per dollar.
10. Mutilate the algorithm used in the parallel implementation to match the architecture.
11. Measure parallel run times on a dedicated system, but measure conventional run times ina busy environment.
12. If all else fails, show pretty pictures and animated videos, and don’t discuss performance.

1992 paper: Scaling performance results to full-sized system
In some published papers and conference presentations, performance results on
small-sized parallel systems were linearly scaled to full-sized systems, without evenclearly disclosing this fact.
Example: 8,192-CPU performance results were linearly scaled to 65,536-CPU results,
simply by multiplying by eight.
Excuse: “We can’t afford a full-sized system.”This and the other examples mentioned in the next few viewgraphs are presented in:
I D. H. Bailey, “Misleading performance reporting in the supercomputing field,” Scientific Programming, vol. 1., no. 2(Winter 1992), pg. 141–151.

1992 paper: Using inefficient algorithms on highly parallel systems
In many cases, inefficient algorithms were employed for the highly parallel
implementation, requiring many more operations, thus producing artificially highMflop/s rates:
I Numerous researchers cited parallel PDE performance based explicit schemes,

where implicit schemes were known to be much better.
Excuse: Explicit schemes “run better” on the researchers’ parallel system.
I One paper cited performance for computing a 3D discrete Fourier transform bydirect evaluation of the defining formula (8n2 operations), rather than by using afast Fourier transform (5n log2 n).
Excuse: Direct computation was “more appropriate” for the architecture being
analyzed.
Both examples violate a rule of professional performance reporting, namely to base theoperation count (when computing Mflop/s or Gflop/s rates) on the best practical serial
algorithm, no matter what scheme was actually used on the parallel system.

1992 paper: Not actually performing a claimed computation
Abstract of published paper: “The current Connection Machine implementation runs
at 300-800 Mflop/s on a full [64K] CM-2, or at the speed of a single processor of a
Cray-2 on 1/4 of a CM-2.”
I Excerpt from text: “This computation requires 568 iterations (taking 272

seconds) on a 16K Connection Machine.”
In other words, the computation was run on a 16K system, not on a 64K system;
the figures cited in the Abstract were merely multiplied by four.
I Excerpt from text: “In contrast, a Convex C210 requires 909 seconds to compute

this example. Experience indicates that for a wide range of problems, a C210 isabout 1/4 the speed of a single processor Cray-2.”
In other words, the computation mentioned in the Abstract was not actually run
on a Cray-2; instead, it was run on a Convex system, and a questionable
rule-of-thumb scaling factor was used to produce the Cray-2 rate.

1992 paper: Performance plot [parallel (lower) vs vector (upper)]3.5

3Time (Hours)

2.5

2

1.5

1

0.5

0

0

1000

2000

3000

4000

5000

6000

Number of Objects

7000

8000

9000

10000

1992 paper: Data for performance plot
Problem size(x axis)
20
40
80
160
990
9600

Parallel system
run time
8:18
9:11
11:59
15:07
21:32
31:36

Vector system
run time
0:16
0:26
0:57
2:11
19:00
3:11:50*

Details in text of paper:
I In last entry, the 3:11:50 figure is an “estimate.”
I The vector system code is “not optimized.”

Note that the parallel system is actually slower than the vector system for all cases,
except for the last (estimated) entry. Also, except for the last entry, all real data in thegraph is in the lower left corner.

Fast forward to 2014: New ways to fool the masses
I Citing performance rates for a run with only one processor core active in a shared-

memory multi-core node, producing artificially inflated performance (since there is
no shared memory interference) and wasting resources (since most cores are idle).
I

Example: Citing performance on “1024 cores,” even though the code was run on1024 multicore nodes, using only one core per node, with 15 out of 16 cores idle.

I Claiming that since one is using a graphics processing unit (GPU) system, that

efficient parallel algorithms must be discarded in favor of more basic algorithms.
I Citing performance rates only for a core algorithm (such as FFT or linear systemsolution), even though full-scale applications have been run on the system.
I Listing only the best performance figure in the paper, even though the run wasmade numerous times (recall the experience of pharmaceutical tests).
I Employing special hardware, operating system or compiler settings that are notappropriate for real-world production usage.

Reproducibility in scientific computing:
4. Numerical reproducibility
The ICERM report 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.
I V. Stodden, D. H. Bailey, J. M. Borwein, R. J. LeVeque, W. Rider and W. Stein, “Setting the default to reproducible:
Reproducibility in computational and experimental mathematics,”
http://www.davidhbailey.com/dhbpapers/icerm-report.pdf.

Sensitivity to numerical round-off error
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.
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.

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?

Solutions to numerical sensitivity and reproducibility 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.
Enough excuses!
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.

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

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, May 2013,
http://www.davidhbailey.com/dhbpapers/precimonious.pdf.

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

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 ofSupercomputing, vol. 18, no. 3 (Mar2001), 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

The PSLQ integer relation algorithm
Let (xn ) be a given vector of real numbers. An integer relation algorithm either finds
integers (an ) such thata1 x1 + a2 x2 + · · · + an xn = 0(to within the “epsilon” of the arithmetic being used), or else finds bounds within
which no relation can exist.
The “PSLQ” algorithm of mathematician-sculptor Helaman Ferguson is the mostwidely used integer relation algorithm.
Integer relation detection requires very high precision (at least n × d digits, where d isthe size in digits of the largest ak ), both in the input data and in the operation of thealgorithm.
1. H.R.P. Ferguson, D.H. Bailey and S. Arno, “Analysis of PSLQ, An Integer Relation Finding Algorithm,” Mathematics of Computation, vol.
68, no. 225 (Jan 1999), pg. 351–369.
2. D.H. Bailey and D.J. Broadhurst, “Parallel Integer Relation Detection: Techniques and Applications,” Mathematics of Computation, vol. 70,
no. 236 (Oct 2000), pg. 1719–1736.

Decrease of log 10(min |yi |) in multipair PSLQ run

The BBP formula for π
In 1996, a computer program, utilizing the PSLQ algorithm and high-precision(200-digit) computation, discovered this new formula for π:


∞
X
4
2
1
1
1
−
−
−
π =
16k 8k + 1 8k + 4 8k + 5 8k + 6
n=0This formula permits one to compute binary (or hexadecimal) digits of π beginning at
an arbitrary starting position, using a very simple scheme that requires only standard
64-bit or 128-bit arithmetic.
In 2004, Borwein, Galway and Borwein proved that no base-n formulas of this type
exist for π, except when n = 2m .
BBP-type formulas (discovered with PSLQ) are now known for numerous othermathematical constants.
1. D.H. Bailey, P.B. Borwein and S. Plouffe, “On the rapid computation of various polylogarithmic constants,” Mathematics of Computation,
vol. 66, no. 218 (Apr 1997), pg. 903–913.
2. J.M. Borwein, W.F. Galway and D. Borwein, “Finding and excluding b-ary Machin-type BBP formulae,” Canadian Journal of Mathematics,
vol. 56 (2004), pg. 1339–1342.

Poisson potential functions associated with lattice sums
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 )1log αφ2 (x, y ) =
=
2
2
2
π
m +n
π
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, 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), 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

For k = 32, the minimal polynomial has degree 128, with individual coefficients
ranging up to 56 digits long. This PSLQ computation required 10,000-digit precision.
Other polynomials required up to 50,000-digit precision.

How much precision is enough?

I For any of these high-precision computations, there is no known way to know apriori what precision level should be used.
I Current ad-hoc approach: Increase precision until a consistent, credible-looking

result is achieved, then reduce until one finds roughly the minimum precisionneeded.
I Clearly research is needed to better understand this process.

Summary
Reproducibility looms as a major challenge for the future of high-performancecomputing:
1. The need to institute a culture of reproducibility.
2. The danger of statistical overfitting and other errors in data analysis.
3. The need for greater rigor and forthrightness in performance reporting.
4. The growing concern over numerical reproducibility.
Technical and organizational solutions are available, but will the community adopt
them?

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