Reproducible scientific computing: Progress and challenges
David H. Bailey
Lawrence Berkeley Natlional Laboratory (retired)
University of California, Davis, Department of Computer Science
This talk is available at:
http://www.davidhbailey.com/dhbtalks/dhb-sc2020-correct.pdf

1 / 36

1991 paper: “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.

2 / 36

Reproducibility crises in biomedicine, psychology, finance and economics
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 the past few years, “backtest overfitting” has

Reproducibility Project staff
Credit: NY Times

emerged as a major issue in computational finance.
http://www.nytimes.com/2015/08/28/science/
many-social-science-findings-not-as-strong-as-claimed-study-says.html
3 / 36

Backtest overfitting in computational financeBacktest overfitting: Statistical overfitting of historical market data. Examples:
I Proposing a theoretical model that inherently possesses a higher level ofcomplexity than the backtest data.
I Using a computer to try millions or billions of variations of a strategy on historicalmarket data, and then only presenting results from the one that works best.
I Constructing an investment fund by using a computer to explore millions orbillions of weighting factors, then only marketing the one that works best.
When a computer program tries many variations of a fund or strategy on a fixed
dataset, it is almost certain that the “optimal” selection will be statistically overfit.
Backtest overfitting is now thought to be the principal reason why financial funds andstrategies that look great on paper often disappoint when actually fielded.
I D. H. Bailey, J. M. Borwein, M. Lopez de Prado and J. Zhu, “Pseudo-mathematics and financial
charlatanism: The effects of backtest overfitting on out-of-sample performance,” Notices of the American
Mathematical Society, May 2014, 458–471.
4 / 36

ICERM report on reproducibility in mathematical and scientific computing
Unfortunately, the scientific culture surrounding computational work has evolved inways that often make it difficult to verify findings, efficiently build on past research, oreven to apply the basic tenets of the scientific method to computational procedures.
Bench scientists [in other fields] 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 is often 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 to obtain
the published results, let alone a record of the false starts. This ultimately has adetrimental effect on researchers’ own productivity, their ability to build on past results
or participate in community efforts, and the credibility of the research among otherscientists and the public.
I V. Stodden, D. H. Bailey, J. M. Borwein, R. J. LeVeque, W. Rider and W. Stein, “Setting the default toreproducible: Reproducibility in computational and experimental mathematics,” manuscript, 2 Feb 2013,https://www.davidhbailey.com/dhbpapers/icerm-report.pdf.

5 / 36

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.
6 / 36

Some all-too-common practices in the HPC community
I Did you only run the code once?

“We ran the code ten times, but only published the best results in our paper.”
I Where are the full details of your computational techniques?

“We do not share these details — they are our competitive edge over otherresearch teams...”I What was the system environment used for your runs?
“We used whatever the campus computer center had installed at the time...”I Where is the source code?
“A student wrote the code, but he/she graduated.”
I Which version of the source code did you use?
“That is no longer available — after running the code for our conference
submission, we immediately started making some additional modifications...”I How can you be sure that you used sufficient numeric precision?
“We just used standard IEEE floating-point arithmetic. Isn’t that good enough?”
I Have you compared your run results with those of other research teams?
“Our grant did not provide funding for reproducibility analyses.”
7 / 36

Numerical reproducibility
The ICERM report specifically addressed numerical reproducibility:
Numerical reproducibility has emerged as a particularly important issue, since the scaleof computations has greatly increased in recent years, particularly with computations
performed on many thousands of processors and involving similarly large datasets.
Large computations often greatly magnify the level of numeric error, so that numericaldifficulties that were once of little import now are large enough to alter the course of
the computation or to draw into question the overall validity of the results.

8 / 36

Reproducibility problems in a Large Hadron Collider code
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: five million line of C++ and Python code, developed by roughly 2000
physicists and engineers over 15 years.
In an attempt to speed up the calculation, researchers found that merely changing theunderlying math library (which should only affect at most the last bit) resulted in somecollisions being missed or misidentified.
Questions:
I How extensive are these numerical difficulties?
I How can numerically sensitive code be tracked down?
I How can a large code library such as this be maintained, producing numericallyreliable and reproducible results?
9 / 36

How to solve accuracy problems and ensure reproducibility
1. Employ an expert numerical analyst to examine every algorithm employed in thecode, to ensure that only the most stable and efficient schemes are being used.
2. Employ an expert numerical analyst to analyze every section of code for numericalsensitivity.
3. Employ an expert numerical analyst to convert large portions of the code to useinterval arithmetic (greatly increasing run time and code complexity).
4. Employ significantly higher precision than is really necessary (greatly increasing
run time).
5. Employ variable-precision arithmetic, assisted with some “smart” tools to help
determine where extra precision is needed and where it is not.
Item 5 is the only practical solution for real-world scientific computing.

10 / 36

Numerical analysis expertise among U.C. Berkeley graduates
Of the 2010 U.C. Berkeley graduating class, 870 were indisciplines requiring technical computing (count by DHB):
I Division of Mathematical and Physical Sciences (Math,

Physics, Statistics).
I College of Chemistry.
I College of Engineering (including Computer Science).

Other fields whose graduates will likely do significant computing:
I Finance, biology, geology, medicine, economics, psychology and sociology.
The total count is very likely over 1000; probably closer to 2000.
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: At most, only about 2% of U.C. Berkeley graduates who will do technicalcomputing in their careers have had rigorous, expert-level training in numerical analysis. 11 / 36

Variable precision computing
I Mixed 16-bit/32-bit: Artificial intelligence, machine learning, graphics.
I Mixed 32-bit/64-bit: A broad range of current scientific applications.
I Mixed 64-bit/128-bit: Large applications with numerically sensitive portions.
I Mixed 64-bit/multiple: Applications in computational mathematics and physics.

12 / 36

Commonly used formats for floating-point computing
Formalname
IEEE 16-bit(none)(none)IEEE 32-bit
IEEE 64-bit
IEEE 80-bit
IEEE 128-bit(none)(none)(none)(none)

Nickname
IEEE half
ARM half
bfloat16
IEEE single
EEE double
IEEE extended
IEEE quad
double-double
quad-double
double-quad
multiple

Sign
1
1
1
1
1
1
1
1
1
1
variesNumber of bits
Exponent Mantissa
5
10
5
10
8
7
7
24
11
52
15
64
15
112
11
104
11
208
15
224
variesvaries

Hidden
1
1
1
1
1
0
1
2
4
1
variesDigits
3
3
2
7
15
19
34
31
62
68
varies13 / 36

Advantages of variable precision
Compared with using a fixed level of precision for the entire application, employing avariable precision framework features:
I Faster processing.
I Better cache utilization.
I Lower run-time memory usage.
I Lower offline data storage.
I Lower energy costs.

Challenge:
I Unless the precision level is carefully managed, there is significantly greater

potential for loss of numerical accuracy and reproducibility.

14 / 36Mixed half-single applications: Machine learning
NVIDIA researchers found that a combination of IEEE 32-bit and bfloat16 arithmetic
achieved nearly the same training loss reduction curve on a English language learningmodel as with 100% IEEE 32-bit:
Maintain a master copy of training system weights in IEEE 32-bit, and select a scaling
factor S. Then for each iteration:
I Convert the array of 32-bit weights to a bfloat16 array.
I Perform forward propagation with bfloat16 weights and activations.
I Multiply the resulting loss with a scaling factor S.
I Perform backward propagation with bfloat16 weights, activations and gradients.
I Multiply the weight gradient by 1/S.
I Update the 32-bit weights using the bfloat16 data.

I “Deep learning SDK documentation,” NVIDIA, 2019,https://docs.nvidia.com/deeplearning/sdk/index.html.
15 / 36

DeepMind’s AlphaGo Zero teaches itself to play Go
I In October 2017, DeepMind researchers programmed amachine learning system with the rules of Go, then had theprogram play itself — teaching itself with no human input.
I After just three days of training, the resulting program“AlphaGo Zero” defeated their earlier program (which had
defeated Ke Jie, the highest rated human) 100 games to 0.
I The Elo rating of Ke Jie, the world’s highest rated Go

player, is 3661. After 40 days of training, AlphaGo Zero’s
Elo rating was over 5000.
I AlphaGo Zero was as far ahead of the world champion as

the world champion is ahead of a good amateur.
I “The latest AI can work things out without being taught,” The Economist, 21 Oct 2017,https://www.economist.com/news/science-and-technology/
21730391-learning-play-go-only-start-latest-ai-can-work-things-out-without.
16 / 36Challenges for mixed half-single computing
I Software availability: So far there are relatively few specialized packages andwell-supported high-level language interfaces for 16-bit computing.
I Hardware compatibility: The proliferation of hardware formats (IEEE half,bfloat16, arm16) complicates the development of software.
I Reproducibility: Only eight mantissa bits severely limits any measure of accuracy

and reliability. At the least, software must provide a simple means to re-run anapplication using 32-bit or 64-bit precision to verify some level of reliability andreproducibility.

17 / 36Mixed single-double application: LU decomposition
I The usual scheme is to factor a pivoted coefficient matrix PA = LU, where L islower triangular and U is upper triangular, then solve Lu = Pb and Ux = y , allusing double precision.
I In a mixed single-double scheme, the factorization and solutions are done in single
precision; then the residual and full solutions are done in double.
I Note that the factorization of A, which is the only step that is O(n3 ), is done
entirely in single.
I A similar scheme can be used for sparse systems.
I A broad range of iterative refinement methods can be performed in this way.
I Baboulin et al. (below) achieved speedups of 1.86 on a direct sparse system on asuite of test problems.

I M. Baboulin, A. Buttari, J. Dongarra, J. Kurzak, J. Langou, P. Luszczek and S. Tomov, “Accelerating
scientific computations with mixed precision algorithms,” Computational Physics Communications, vol.
180 (2009), 2256–2253.
18 / 36Mixed double-quad application: Computational biology
I Certain models of biochemical reaction networks, e.g., “metabolic expression

models,” involve solving large multiscale linear programming problems.
I Since double precision implementations of these models often produce numericallyunreliable results, researchers have resorted to exact rational arithmetic methods,
which are enormously expensive — run times are typically weeks or months.
I Ding Ma and Michael Sanders of Stanford have achieved excellent results using
double codes enhanced with quad precision for numerically sensitive operations.
I Their code was a straightforward modification of a double precision Fortran code,using the REAL(16) datatype, which is supported in the GNU gfortran compiler.

I D. Ma and M. Saunders, “Solving multiscale linear programs using the simplex method in quadruple
precision,” in M. Al-Baali, L. Grandinetti and A. Purnama, ed., Recent Developments in Numerical
Analysis and Optimization, Springer, NY, 2017,
http://web.stanford.edu/group/SOL/reports/quadLP3.pdf.
19 / 36

Innocuous example where IEEE 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.
Common approach: Use a least-squares scheme, solving the linear system:






P n

n
k=1 xk..
Pn .

n
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    ..  = .
..
..
 .  Pn . nPn . 2nan
k=1 xk yk
k=1 xkAn implementation using IEEE 64-bit arithmetic fails to find the correct polynomial,
even if one rounds results to nearest integer.
But an implementation using IEEE quad or double-double arithmetic quickly produces
the correct polynomial:f (x ) = 1 + 1048577x 4 + x 8 = 1 + (220 + 1)x 4 + x 8 .
20 / 36

Innocuous example, cont.
The result on the previous page can be obtained with 64-bit arithmetic using Lagrange
interpolation or the Demmel-Koev algorithm. But few outside of expert numericalanalysts are aware of these schemes.
Besides, even these two schemes fail for larger problems, such as:(1, 134217731, 8589938753, 97845255883, 549772595201,
2097396156251, 6264239146561, 15804422886323, 35253091827713,
71611233653971, 135217729000001, 240913322581691, 409688091758593).
These values are generated by the simple polynomialf (x ) = 1 + 134217729x 6 + x 12 = 1 + (227 + 1)x 6 + x 12 .
This is easily found by an IEEE quad or double-double implementation of any of thethree algorithms (least-squares, Lagrange or Demmel-Koev).

21 / 36Mixed double/multiprecision applications
1. Solutions of ordinary differential equations in planetary orbit models (32 digits).
2. Supernova simulations (32–67 digits).
3. Scattering amplitudes of fundamental particles (32 digits).
4. Solutions of certain discrete dynamical systems (32 digits).
5. Coulomb n-body atomic system simulations (32–120 digits).
6. Solving ordinary differential equations using the Taylor algorithm (100–300 digits).
7. Recognizing integrals from the Ising theory of mathematical physics in terms of
simple mathematical identities (200–1000 digits).
8. Finding minimal polynomials connected to the Poisson equation of mathematicalphysics (1000–64,000 digits).

I D. H. Bailey, R. Barrio, and J. M. Borwein, “High precision computation: Mathematical physics anddynamics,” Applied Mathematics and Computation, vol. 218 (2012), 10106–10121.
22 / 36Mixed double/multiprecision: Computational math and physics
Frequently-used “experimental math” methodology:
1. Compute a limit, infinite series sum, definite integral or other mathematicalquantity to high precision, typically 100–10,000 digits.
2. Use an integer relation algorithm such as “PSLQ” to recognize the numericalvalue in terms of a mathematical expression. In other words, if the computed
quantity is α, and if one suspects that α = x0 is given by a formula involving
known constants (x1 , x2 , · · · , xn−1 ), then use PSLQ to find integers(a0 , a1 , a2 , · · · , an−1 ) such that a0 x0 + a1 x1 + · · · + an−1 xn−1 = 0.
3. When a result is found experimentally, seek a formal mathematical proof of the
discovered relation.
4. Efficient implementations typically require three levels of precision, such as
double/500-digit/10,000-digit, varying up and down many times within the run,performing as much computation as possible using only double precision.

23 / 36Experimental mathematics and physics results
I New formulas and identities for various mathematical constants, including

π, log(2) and others, some with remarkable properties (250–1000 digits).
I Evaluation of Euler sums, used in quantum physics (250–1000 digits).
I Evaluation of Ising integrals, used in quantum physics and string theory(250–2000 digits).
I Evaluation of Poisson polynomials, used in mathematical physics and image

processing (1000–64,000 digits).

I D. H. Bailey and J. M. Borwein, “High-precision arithmetic in mathematical physics,” Mathematics, vol.
3 (2015), 337–367.
24 / 36Mixed double/multiprecision application: Finding Poisson polynomials
In 2012 Richard Crandall of Apple Computers found that each pixel in the image
produced by a certain sharpening algorithm is given by the 2-D Poisson potentialfunction of mathematical physics:
1 X cos(mπx ) cos(nπy )φ2 (x , y ) = 2
π m,n odd
m2 + n2
Crandall and others then found that for rational x and y ,1φ2 (x , y ) = log απwhere α is algebraic, i.e., the root of a some integer polynomial of degree m.
What is the connection between x and y and the degree of α?
I 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), 115201.
I D. H. Bailey, J. M. Borwein, J. Kimberley and W. Ladd, “Computer discovery and analysis of large
Poisson polynomials,” Experimental Mathematics, 27 Aug 2016, vol. 26, 349–363.
25 / 36

Finding Poisson polynomials, cont.
By computing high-precision numerical values of φ2 (x , y ) for various specific rational x
and y , using three precision levels (e.g., double/500-digit/10,000-digit), and applying avariant of the PSLQ algorithm, we found the following results among others:
s
5
6
7
8
9

10

Minimal polynomial of α corresponding to x = y = 1/s
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

26 / 36

Finding Poisson polynomials, cont.
Based on these preliminary results, Jason Kimberley conjectured that the degree m(s)of the minimal polynomial associated with the case x = y = 1/s is given by this rule:Set m(2) = 1/2. Otherwise for primes p congruent to 1 mod 4, setm(p) = int2 (p/2), where int denotes greatest integer, and for primes pcongruent to 3 mod 4, set m(p) = int (p/2)(int (p/2) + 1). Then for any
other positive integer s whose prime factorization is s = p1e1 p2e2 · · · prer ,r −1m(s) = 4rY
2(ei −1)pim(pi ).

i=1

Additional three-level precision computations, (e.g., double/2000-digit/64,000-digit),
confirmed Kimberley’s formula for all s up to 40, and for even s up to 52.
By examining the computed results, connections were found to a sequence ofpolynomials defined in a 2010 paper by Savin and Quarfoot. These observations
ultimately led to a proof of Kimberley’s formula and several other facts.
27 / 36

Producing reliable results in large variable-precision calculations
Producing reliable results in this type of computation requires considerable effort andmachinery:
I The current accuracy and potential for exhaustion of accuracy must becontinually monitored all levels of precision.
I As a lower level of precision is nearing exhaustion, higher-precision arrays must be
updated, and then the lower-level arrays re-initialized.
I If, in spite of careful monitoring, the program detects that precision at one levelhas been exhausted, perhaps due to a statistically unlikely event, there must be afacility to “back up” the calculation to a previous, certified iteration.
I For large three-level calculations, such as the Poisson polynomial studies, all of
this machinery increases the length and complexity of the resulting code by afactor of four, compared with a straightforward implementation using only thehighest level of precision.

28 / 36

Software and language issues: Conversions between types
Existing programming languages have significant issues with varying precision levels,
particularly in conversions between types. For example, consider:program testprec
real(8) d1, d2, d3
d1 = sqrt (3. + 1./7.)
d2 = sqrt (3.d0 + 1./7.)
d3 = sqrt (3. + 1.d0/7.)
write (6, ’(f20.15)’) d1, d2, d3
stop
end

This code, when compiled with the gfortran compiler with default options, produces:
1.772810459136963
1.772810522656991
1.772810520855837

Only the last is a correct, full-precision value!
One could set a compiler flag to treat all single constants as double (default in somelanguages), but this defeats the goal of using varying precision levels.

29 / 36

Software and language issues: Conversions between types, cont.
Similar difficulties, greatly magnified, afflict mixed double-multiprecision applications:program testprec
use mpmodule
type (mp real) r1, r2, r3
r1 = sqrt (3.d0 + 1.d0/7.d0)
r2 = sqrt (mpreal (3.d0, 10) + 1.d0/7.d0)
r3 = sqrt (3.d0 + mpreal (1.d0, 10)/7.d0)
call mpwrite (6, 80, 60, r1, r2, r3)
stop
end

DHB’s MPFUN-MPFR package checks for conversion problems such as this atruntime, producing an error message like this:*** MP CHECKDP: DP value has more than 40 significant bits:
1.772810520855837D+00
and thus very likely represents an unintended loss of accuracy.

One solution: A compiler option to outlaw all mixed-mode operations and automatic
type conversions. Unpopular? Perhaps. But finding precision problems is even harder.
30 / 36Thread safety in high precision libraries
Most multiprecision libraries in use today are not thread-safe. Thus applications
written to use them do not work in a shared-memory parallel environment, such as ona single node of a large highly parallel calculation.
The solution is to include the following in each multiprecision datum:
I The full length of the storage allocated to this datum.
I The current working precision associated with this datum.
I The length of nonzero data associated with this datum.
Currently the MPFR package (see below) is the only low-level multiprecision librarythat is thread-safe (provided thread-safe option is set); it is also the fastest.
DHB’s MPFUN-MPFR high-level Fortran multiprecision package is built on MPFR.

I L. Fousse, G. Hanrot, V. Lefevre, P. Pelissier and P. Zimmermann, “MPFR: A multiple precision binary
floating-point library with correct rounding,” ACM Transactions on Mathematical Software, vol. 33(2007), https://doi.org/10.1145/1236463.1236468.

31 / 36

The U.C. Berkeley/U.C. Davis “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).
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, November 2013,https://www.davidhbailey.com/dhbpapers/precimonious.pdf.
I C. Nguyen, C. Rubio-Gonzalez, B. Mehne, K. Sen, C. Iancu, J. Demmel, W. Kahan, W. Lavrijsen,
D. H. Bailey and D. Hough, “Floating-point precision tuning using blame analysis,” Proceedings of the
38th International Conference on Software Engineering (ICSE 2016), 14–22 May 2016,
https://www.davidhbailey.com/dhbpapers/blame-analysis.pdf.
32 / 36

Bit-for-bit floating-point replicability between systems
I Although the IEEE floating-point standards guarantee bit-for-bit results (provided

the same rounding mode is used) for individual operations, different systems may
produce different higher-level results, because of non-associativity:(A ⊕ B) ⊕ C 6= A ⊕ (B ⊕ C ).
I Recently some researchers have proposed an addition to the IEEE-754 standard
that would permit guaranteed bit-for-bit replicability, not only on a single systembut also on a parallel system — see references below.
I ARM (a U.K. computer processor vendor) has experimented with a high-precision
accumulator, using an adaptable vector unit.
Dangers: If this standard becomes the default on computer systems:
I Will it lock into place computations that are seriously inaccurate?
I Will it make it more difficult for users to even be aware that they have a problem
with numerical sensitivity and/or inaccuracy?
I P. Ahrens, J. Demmel and H. D. Nguyen, “Algorithms for efficient reproducible floating pointsummation,” manuscript available from authors.
I J. Riedy and J. Demmel, “Augmented arithmetic operations proposed for IEEE-754 2018,”
http://www.ecs.umass.edu/arith-2018/pdf/arith25_34.pdf.

33 / 36

New approaches to floating-point arithmetic
I Although the IEEE-754 floating-point standard has served us very well for several

decades, some wonder if it time to completely rethink floating-point arithmetic.
I Gustafson and Yonemoto have proposed “Unum,” a floating-point standard that is

flexible in size and precision — see references below.
I Lindstrom summarizes Unums and some other proposals — see below.

I J. L. Gustafson, The End of Error: Unum Computing, Chapman and Hall, New York, 2017.
I J. Gustafson and I. Yonemoto, “Beating floating point at its own game: Posit arithmetic,” 2017,
http://www.johngustafson.net/pdfs/BeatingFloatingPoint.pdf.
I P. Lindstrom, “Fixed-rate compressed floating-point arrays,” IEEE Transactions on Visualization andComputer Graphics, vol. 20 (2014), 2674–2683, https://doi.org/10.1109/TVCG.2014.2346458.
34 / 36

Research questions to be explored
1. How are various scientific disciplines currently using variable precision? What newopportunities are on the horizon?
2. What does the explosive growth in artificial intelligence and machine learning mean for
hardware, software and tools to support variable precision? What are the fundamental
accuracy requirements of such computations?
3. Are there new mathematical approaches to efficiently and accurately compute various
mathematical operations and library functions to both very modest precision (2-3 digit
accuracy) and very high precision (hundreds or thousands of digits)?
4. What new formats and representation systems to facilitate variable level precisioncomputing are available? To what extent do these new proposals overcome strengths andweaknesses in existing systems such as IEEE-754 hardware and MPFR software?
5. What new support is needed from computer hardware and software vendors for variable
precision computing? Could some modest hardware changes be made to support low,
variable precision and high precision, rather than relying on all-software implementations?
6. Almost all of published research in numerical mathematics and numerical analysis to datehas focused on either 32-bit or 64-bit results. Are there other techniques that may be
35 / 36superior in a variable precision environment?

Research questions to be explored, cont.
6. Classical error bounds used in numerical analysis, such as those developed by Wilkinson,
typically assume worst-case scenarios and thus are not very useful in a limited precision orvariable precision environment. Can more realistic error estimates be derived?
7. How can iterative refinement best be performed in a variable precision environment?
What new opportunities are there for using iterative refinement with variable precision?
8. What algorithms and software techniques can effectively compress floating-point data,on-the-fly, to save data storage and transfer? Can improved schemes be devised?
9. What software tools are available to automatically or semi-automatically analyze an
application code, identify numerically sensitive spots and propose or implement remedies?
What other tools are needed?
10. How can we ensure that computations with variable precision levels are reproducible?
What does reproducibility even mean when computing, say, with only 16-bit precision?

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