Experimental Mathematics and High-Performance Computing David H - - PowerPoint PPT Presentation

Experimental Mathematics and High-Performance Computing

David H Bailey Lawrence Berkeley National Lab

“All truths are easy to understand once they are discovered; the point is to discover them.” – Galileo Galilei

The NERSC Computer Center at the Berkeley Laboratory

Seaborg: 6656-CPU IBM P3 system, 10 Tflop/s peak, 7.8 Tbye memory. Bassi: 976-CPU IBM P5 system, 6.7 Tflop/s peak, 3.5 Tbyte memory. Franklin: 9672 dual-core Opteron CPUs, 100 Tflop/s peak is now being installed.

Example of NERC Computations: Astrophysics

• Multi-physics and multi-scale phenomena.
• Large dynamic range in time and length.
• Requires adaptive mesh refinement.
• Dense linear algebra.
• FFTs and spherical harmonic transforms.

Supernova simulation:

• Future 3-D model calculations will require 1,000,000

CPU-hours per run, on 100 Tflop/s peak system.

Analysis of cosmic microwave background data:

• WMAP (now)

3x1021 flops, 16 Tbyte mem

• PLANCK (2007)

2x1024 flops, 1.6 Pbyte mem

• CMBpol (2015)

1x1027 flops, 1 Ebyte mem

Graphic: T. Mezzacappa, J. Blondin, K.-L. Ma, et al (ORNL)

Characteristics of Modern High- Performance Scientific Computing

The ultimate objective is to advance the applied discipline:

• Physics, chemistry, astronomy, biology, climate, engineering, biotech.

Advanced numerical algorithms and computing techniques:

• FFTs, dense linear algebra, sparse linear algebra, iterative solvers,

multigrid, highly parallel processing, dynamic data structures, etc.

State-of-the-art calculations require highly parallel computers:

• Enormous computational requirements are common.
• 1000+ CPUs are used in many calculations.

A pragmatic attitude prevails: “If it works, use it.”

• Some combinatorial optimization algorithms are observed to work

significantly better in practice than theory might suggest.

• Gaussian elimination with partial pivoting is not guaranteed to work in

all cases, yet it works fine in real applications.

• The QR algorithm was used for many years before it was found to

cycle in a simple 4x4 case. A proof of convergence is still elusive.

What Is Experimental Mathematics?

“Experimental mathematics” is a term for the emerging discipline where state-of-the-art computing technology is aggressively applied to problems in mathematical research: Actively exploring mathematical questions. Computing explicit numerical examples. Performing large symbolic manipulations. Testing (and often falsifying) conjectures. Investigating possible paths for formal proof. Hamming: “The purpose of computing is insight, not numbers.”

Books on Experimental Mathematics

Mathematics by Experiment: Plausible Reasoning in the 21st Century Experiments in Mathematics: Computational Paths to Discovery Authors: Jonathan Borwein, DHB and (for vol. 2) Roland Girgensohn. Coming soon (Mar 2007): Experimental Mathematics in Action. Authors: David Bailey, Jon Borwein, Neil Calder, Roland Girgensohn, Russell Luke and Victor Moll. Both books are now available on CD-ROM in a hyperlinked, searchable PDF format. Also, a FREE condensed version is available at:

http://www.experimentalmath.info

Experimental Mathematics as High-Performance Computing

The ultimate objective is to advance the applied discipline:

• Here the “applied discipline” is pure mathematics!

Advanced numerical algorithms and computing techniques:

• PSLQ, high-precision arithmetic, symbolic computing, FFTs, numerical

analysis, evaluation of integrals and series, etc.

State-of-the-art calculations require highly parallel computers:

• High-precision arithmetic greatly magnifies run times.
• 1000+ CPUs have been used in several calculations.

A pragmatic attitude prevails: “If it works, use it.”

• We do not know ahead of time what terms to use in an integer relation

search – guessing which terms to try is still a black art.

• Whereas the standard PSLQ algorithm is guaranteed to find relations,

no proof is known for multi-pair PSLQ.

• We do not fully understand why tanh-sinh quadrature works so well,

especially in 2-D, 3-D, etc.

Examples of Large Experimental Math Computations

Identification of the 4th bifurcation point of the logistic iteration:

Integer relation of size 121; 10,000 digit arithmetic. Required 67 min on 48 CPUs = 54 CPU-hours.

Finding a relation derived from roots of Lehmer’s polynomial:

Integer relation of size 125; 50,000 digit arithmetic. Required 16 hours on 64 CPUs = 1024 CPU-hours.

Numerical verification of a mathematical physics integral:

1-D quadrature calculation; 20,000-digit arithmetic. Required 45 min on 1024 CPUs = 768 CPU-hours.

Numerical evaluation of an Ising theory integral:

3-D quadrature of a very complicated function; 500-digit arithmetic. Required 18.2 hours on 256 CPUs = 4659 CPU-hours.

Authors: David Broadhurst, Jonathan Borwein, Richard Crandall, Roland Girgensohn and DHB

Computational Methods Used in Experimental Math

High-precision computation. PSLQ (integer relation detection). Symbolic computing tools. Function evaluations: sin, exp, log, erf, gamma, zeta, polylog. Fast Fourier transforms (FFTs). Dense and sparse linear algebra. Evaluation of definite integrals. Evaluation of infinite series sums. Error bounds on computed results. Highly parallel computing. Computer graphics. Note that except for the first three, these are all staples of modern applied mathematics and numerical analysis.

LBNL’s High-Precision Software (ARPREC and QD)

Low-level routines written in C++. C++ and F-90 translation modules permit use with existing programs with only minor code changes. Double-double (32 digits), quad-double, (64 digits) and arbitrary precision (>64 digits) available. Special routines for extra-high precision (>1000 dig). High-precision integer, real and complex datatypes. Includes many common functions: sqrt, cos, exp, gamma, etc. PSLQ, root finding, numerical integration. An interactive “Experimental Mathematician’s Toolkit” is also available. Available at: http://www.experimentalmath.info This software is being used by physicists, climate modelers, chemists and engineers, in addition to mathematicians.

Authors: Xiaoye Li, Yozo Hida, Brandon Thompson and DHB

ARPREC vs GMP

ARPREC advantages:

Comparatively simple install procedure. Simple arrays facilitate parallel implementations. High-level Fortran-90/95 interface (not available for GMP). High-level C++ interface (ARPREC’s is nicer than GMP’s). FFT-based arithmetic for very high precision (> 1000 digits).

GMP/MPFR advantages:

Better performance, especially for over 1000 digit precision and for transcendental functions. Support of a large community.

What is needed: Combine the high-level ARPREC Fortran and C++ interfaces with the GMP low-level routines. Issue: How can this be done and still facilitate parallel applications?

The PSLQ Integer Relation Algorithm

Let (xn) be a vector of real numbers. An integer relation algorithm finds integers (an) such that At the present time, the PSLQ algorithm of mathematician-sculptor Helaman Ferguson is the best-known integer relation algorithm. PSLQ was named one of ten “algorithms of the century” by Computing in Science and Engineering. High-precision arithmetic software is required: at least d x n digits, where d is the size (in digits) of the largest of the integers ak.

Authors: Helaman Ferguson, Stephen Arno and DHB

The BBP Formula for Pi

In 1996, a computer program running the PSLQ algorithm discovered this formula for pi: This formula permits one to directly calculate binary or hexadecimal (base-16) digits of pi beginning at an arbitrary starting position n, without needing to calculate any of the first n-1 digits, by means of a very simple algorithm that requires almost no memory. This formula is now used in the G95 compiler.

Authors: Peter Borwein, Simon Plouffe and DHB

Some Other Similar BBP-Type Identities

Authors: Peter Borwein, Simon Plouffe, David Broadhurst, Richard Crandall and DHB

Is There a Base-10 Formula for Pi?

Note that there is both a base-2 and a base-3 BBP-type formula for π2. Base-2 and base-3 formulas are also known for a handful of other constants. Question: Is there any nonbinary (base-n, where n ≠ 2b) BBP-type formula for π? Answer: No. This is ruled out in a 2004 paper. This does not rule out some completely different scheme for finding individual non-binary digits of π.

Authors: Jon Borwein, David Borwein and Will Galway

Normality

A real number x is said to be b-normal (or normal base b) if every m-long string of base-b digits appears, in the limit, with frequency b-m. Whereas it can be shown that almost all real numbers are b-normal (for any b), there are only a handful of explicit examples. It is not known whether any of the following are b-normal (for any b):

A Connection Between BBP Formulas and Normality

Consider the sequence defined by x0 = 0, and where { } denotes fractional part as before. Result: log(2) is 2-normal if and only if this sequence is equidistributed in the unit interval. In a similar vein, consider the sequence x0 = 0, and Result: p is 16-normal if and only if this sequence is equidistributed in the unit interval.

Authors: Richard Crandall and DHB

A Class of Provably Normal Constants

We have also shown that an infinite class of mathematical constants is normal, including α2,3 was proven 2-normal by Stoneham in 1971, but we have extended this to the case where (2,3) are any pair (p,q) of relatively prime integers. We have also extended this result to an uncountably infinite class, as follows [here rk is the k-th bit of r in (0,1)]: This result has led to a practical and efficient pseudo-random number generator based on the binary digits of α2,3.

Authors: Richard Crandall and DHB

The “Hot Spot” Lemma for Proving Normality (2005)

Recently we were able to prove normality for these alpha constants very simply, by means of a new result that we call the “hot spot” lemma, proven using ergodic theory: Hot Spot Lemma: Let {} denote fractional part. Then x is b-normal if and only if there is no y in [0,1) such that Paraphrase: x is b-normal if and only if it has no base-b hot spots. Sample Corollary: If, for each m and n, no m-long string of digits appears in the first n digits of the base-2 expansion of x more often than 1,000 n 2-m times, then x is 2-normal.

Authors: Michal Misiurewicz and DHB

Another PSLQ Application: Multivariate Zeta Sums

Consider this example: Using the EZFACE+ tool on the CECM website, one obtains the value:

0.1561669333811769158810359096879881936857767098403038729 57529354497075037440295791455205653709358147578...

Using PSLQ, one can then find this evaluation: Dozens of general and specific results have now been established.

Authors: Jonathan Borwein, David Borwein, Roland Girgensohn and DHB

Apery-Like Identities

Authors: Jonathan Borwein and David Bradley

The following were recently found using extensive integer relation searches:

New Apery-Like Identities (Nov 2005)

Authors: Jonathan Borwein, David Bradley and DHB

Following an even more extensive “bootstrapping” experimental approach, similar results have now been found for found for zeta(2n):

The “Multi-Pair” PSLQ Algorithm

Recently a new variant of PSLQ was developed that is well suited for parallel computing (and even runs faster on a single processor). Here are some parallel timings for three benchmark integer relation problems:

Tanh-Sinh Integration

Given f(x) defined on (-1,1), let g(t) = tanh(sinh (t)). Then x = g(t) yields Here xj = g(hj) and wj = g’(hj). Note g’(t) goes to zero very rapidly for large t. Thus even if f(x) has a vertical derivative or blow-up singularity at an endpoint, the product f(g(t)) g’(t) usually is a nice bell-shaped function. For such functions, the Euler-Maclaurin formula implies that the error in the approximation above decreases faster than any power of h. The tanh-sinh scheme often achieves quadratic convergence – reducing h by half produces twice as many correct digits.

Authors: Xiaoye Li, Karthik Jeyabalan and DHB

Parallel Implementation

• f High-Precision Quadrature

The individual function evaluations required for tanh-sinh quadrature (or for other schemes such as Gaussian quadrature) are “embarrassingly parallel.” In most cases, it is NOT necessary to perform individual high- precision arithmetic operations in parallel – high-level, application- loop-level parallelism suffices. Thus highly parallel implementations are relatively straightforward.

Application of Tanh-Sinh Quadrature

This arises from analysis of volumes

• f ideal tetrahedra in hyperbolic
• space. This “identity” has now been

verified numerically to 20,000 digits, but no proof is known. Note that the integrand function has a nasty singularity.

Authors: Jonathan Borwein, David Broadhurst and DHB

Parallel Evaluation

• f the log-tan Integral

1-CPU timings are sums of timings from a 64-CPU run, where barrier waits and communication were not timed. The performance rate for 1024-CPU run is 690 Gflop/s.

Box Integrals

Spurred by a question posed in January 2006 by Luis Goddyn of SFU, we examined some related integrals of the form: The following evaluations are now known: where

Authors: Jonathan Borwein, Richard Crandall and DHB

Ising Integrals

We recently (April 2006) applied our methods to study three classes of integrals that arise in the Ising theory of mathematical physics:

Authors: Jonathan Borwein, Richard Crandall and DHB

Computing and Evaluating Cn

where K0 is the modified Bessel function. We used this formula to compute 1000-digit numerical values of various Cn, from which these results and others were found and proven: Richard Crandall showed that the multi-dimensional Cn integrals can be transformed to 1-D integrals:

Limiting Value of Cn

Cn appear to approach a limit: What is this limit?

Limiting Value of Cn

Cn appear to approach a limit: We pasted the first 50 digits of this numerical value into the Inverse Symbolic Calculator tool, available at

http://oldweb.cecm.sfu.ca/projects/ISC/ISCmain.html

The result was: where gamma denotes Euler’s constant. This experimental result is now proven.

Other Evaluations

The Ising Integral E5

We were able to reduce E5, which is a 5-D integral, to an extremely complicated 3-D integral (see below). We computed this integral to 250-digit precision, using a highly parallel high- precision 3-D quadrature program. Then we used a PSLQ program to discover the evaluation given on the previous page.

Recursions in Ising Integrals

Consider this 2-parameter class of Ising integrals: We computed 1000-digit numerical values for all n up to 36 and all k up to 75 -- a total of 2660 individual quadrature calculations, which were performed independently on a highly parallel computer system. Using PSLQ, we then discovered linear relations in the rows of this array. For example, when n = 3: Similar, but more complicated, recursions were found for larger n (next page).

Experimental Recursion for n = 24

General Recursion Formulas

We were able to find general recursion formulas for each n up to 36:

Compact Recursion Formulas

Let cn,k = n! k! 2-n Cn,k and let M be the largest integer in (n+1)/2. We found (using extensive high-precision polynomial regression) that all of these recursions can be written in the compact form for certain relatively simple polynomials pn,i(x). Here are the polynomials for n = 5 and n = 6:

Spin Integrals

Recently (Jan 2007) we investigated integrals such as: where Note that these integrals involve some complex-arithmetic calculations, even though the final results are real.

Numerical Results for Spin Integrals

A new manuscript by several mathematical physicists asserts that Our parallel multi-dimensional quadrature program affirms P(2) and P(3). We are currently attempting to compute P(4) to high precision.

Conclusions

Experimental mathematics is rapidly becoming a high- performance computing discipline. The huge computing requirements are mostly rooted in the pervasive usage of high-precision arithmetic. In many cases, state-of-the-art results can be obtained in reasonable time only by using a highly parallel computer. Effective parallel implementations have been developed for several key operations, such as PSLQ and numerical integration. BUT we need: Faster high-precision arithmetic, which works properly in a parallel

• environment. Note: in most cases we do NOT need to parallelize

an individual high-precision arithmetic operation. Parallel versions of other useful experimental math computations, such as high-precision infinite series summation.

Some References

Preprints are available at http://crd.lbl.gov/~dhbailey/dhbpapers:

• David H. Bailey, Xiaoye S. Li and Karthik Jeyabalan, "A Comparison of Three High-Precision

Quadrature Schemes," Experimental Mathematics, vol. 14 (2005), no. 3, pg 317-329.

• David H. Bailey, Jonathan M. Borwein, Vishal Kapoor and Eric Weisstein, "Ten Problems in

Experimental Mathematics," American Mathematical Monthly, vol. 113, no. 6 (Jun 2006), pg. 481-409.

• David H. Bailey and Michal Misiurewicz, "A Strong Hot Spot Theorem," Proceedings of the

American Mathematical Society, to appear, 2006.

• David H. Bailey and Jonathan M. Borwein, "Experimental Mathematics: Examples, Methods

and Implications," Notices of the American Mathematical Society, vol. 52, no. 5 (May 2005),

• pg. 502-514.
• David H. Bailey and Jonathan M. Borwein, "Highly Parallel, High-Precision Numerical

Integration," manuscript, Jun 2006.

• David H. Bailey, Jonathan M. Borwein and David M. Bradley, "Experimental Determination
• f Apery-Like Identities for Zeta(2n+2)," Experimental Mathematics, vol. 15 (2006), pg. 281-

289.

• David H. Bailey, Jonathan M. Borwein and Richard E. Crandall, "Box Integrals," Journal of

Computational and Applied Mathematics, to appear, Jun 2006.

• David H. Bailey, Jonathan M. Borwein and Richard E. Crandall, "Integrals of the Ising

Class," Journal of Physics A: Mathematical and General, to appear, Jun 2006.

• David H. Bailey, David Borwein, Jonathan M. Borwein and Richard Crandall,

"Hypergeometric Forms for Ising-Class Integrals," Experimental Mathematics, to appear, Jul 2006.