1 / 54

Experimental Mathematics and Optimization

Experimental Mathematics and Optimization. David H Bailey Lawrence Berkeley National Laboratory http://crd.lbl.gov/~dhbailey. Outline. Introduction Application: the uncertainty principle. Integer relation detection The PSLQ algorithm. Applications: bifurcation constants, sculpture.

odina
Download Presentation

Experimental Mathematics and Optimization

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Experimental Mathematics and Optimization David H Bailey Lawrence Berkeley National Laboratory http://crd.lbl.gov/~dhbailey

  2. Outline • Introduction • Application: the uncertainty principle. • Integer relation detection • The PSLQ algorithm. • Applications: bifurcation constants, sculpture. • Sequence limit extrapolation • Application: the Quinn-Rand-Strogatz constant. • Infinite series summation • The Euler-Maclaurin formula. • Euler’s transformation. • Applications: Euler sums, Apery-like series, Ramanujan-like series. • Numerical integration • Gaussian quadrature. • Tanh-sinh quadrature. • Applications: tetrahedral integral, box integrals, Ising integrals, spin integrals. • Summary

  3. Methods Used in Experimental and Computer-Assisted Mathematics • Computer-based symbolic computation (e.g., Mathematica, Maple, Magma). • High-precision (hundreds or thousands of digits) numerical computation. • Integer relation detection algorithms (notably the PSLQ algorithm). • Computer-based sequence and constant recognition facilities. • Minimization and maximization schemes. • Prime verification algorithms and other computational number theory methods. • The Wilf-Zeilberger algorithm for proving certain infinite series identities. • Formal verification methods – for instance, the effort to verify Thomas Hales’ recent proof of the Kepler conjecture using symbolic logic.

  4. The Uncertainty Principle Uncertainty principle from quantum mechanics: The uncertainty in position times the uncertainty in momentum must be at least some minimum value. Uncertainty principle from signal processing: A signal cannot be both time-limited and band-limited -- the dispersion of the signal times the dispersion of the Fourier transform must be at least some minimum value. The mathematical formulations of these two principles are identical: where

  5. A Numerical Minimization Approach to the Uncertainty Principle A formal proof of the uncertainty principle (UP) is not very difficult, but it is hardly intuitive at first reading. Here is a more computationally and intuitively appealing approach: 1. Start with a simple “tent” function. 2. Evaluate the UP product for the all perturbations of the current function on a grid with the current resolution. 3. Select the function from step 2 that has the smallest UP product. 4. Refine the grid and go to 2. 5. End when grid is sufficiently fine. Minimizing function: Gaussian function Minimum UP product: 0.25 Jonathan M. Borwein and David H. Bailey, Mathematics by Experiment, AK Peters, 2004, pg. 183-188.

  6. The Integer Relation Problem Let (xn) be a given vector of real numbers. An integer relation algorithm finds integers (an)such that (or at least within “epsilon” of zero). This can be viewed as an integer programming problem: find a set of integers (ai) that minimizes subject to the constraint for some M. The input (xi) values must be specified to at least d*n-digit precision, where d is the size (in digits) of the largest of the (ai), and the integer relation detection algorithm must be performed using at least d*n-digit precision.

  7. The PSLQ Integer Relation Algorithm • 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. • A “multi-pair” variant of PSLQ is known that is well-suited for parallel computation. • PSLQ constructs a sequence of integer-valued matrices Bn that reduces the vector y = x * Bn, until either the relation is found (as one of the columns of Bn), or else precision is exhausted. • When a relation is found, the size of smallest entry of the vector y suddenly drops to roughly “epsilon” (i.e. 10-p, where p is the number of digits of precision). • The size of this drop can be viewed as a “confidence level” that the relation is real and not merely a numerical artifact -- a drop of 20+ orders of magnitude almost always indicates a real relation. 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.

  8. Statement of Basic PSLQ Algorithm Initialize: 1. Set the n x n matrices A and B to the identity. 2. For k := 1 to n do: set s(k) := sqrt[sum_(k <j<=n) x(j)^2]; enddo; set t = 1 / s(1); for k := 1 to n do: set y(k) := t x(k); s(k) := t s(k); enddo. 3. Compute the initial n x (n-1) matrix H as H(I,j) = 0 if i < j, H(j,j) := s(j+1) / s(j), and H(I,j) := - y(i) y(j) / (s(j) s(j+1)) if i > j. 4. Reduce H: for i := 2 to n do: for j := i - 1 to 1 step -1 do: set t := nint (H(i,j) / H(j,j)); and y(j) := y(j) + t y(i); for k := 1 to j do: set H(i,k) := H(i,k) - t H(j,k); enddo; for k := 1 to n do: set A(i,k) := A(i,k) - t A(j,k) and B(k,j) := B(k,j) + t B(k,i); enddo; enddo; enddo. Iterate the following until an entry of y is within a reasonable tolerance of zero, or precision has been exhausted: 1. Select m such that gamma^i |H(i,i)| is maximal when i = m. 2. Exchange the entries of y indexed m and m + 1, the corresponding rows of A and H, and the corresponding columns of B. 3. Remove the corner on H diagonal: If m <= n - 2 then set t(0) := sqrt[H(m,m)^2 + H(m,m+1)^2)], t(1) := H(m,m) / t(0) and t(2) := H(m,m+1) / t(0); for i := m to n do: set t(3) := H(i,m), t(4) := H(i,m+1), H(i,m) := t(1) t(3) + t(2) t(4) and H(i,m+1) := - t(2) t(3) + t(1) t(4); enddo; endif. 4. Reduce H: for i := m + 1 to n do: for j := min(i - 1, m + 1) to 1 step -1 do: set t := nint (H(i,j) / H(j,j)) and y(j) := y(j) + t y(i); for k := 1 to j do: set H(i,k) := H(i,k) - t H(j,k); enddo; for k := 1 to n do: set A(i,k) := A(i,k) - t A(j,k) and B(k,j) := B(k,j) + t B(k,i); enddo; enddo; enddo. 5. Norm bound: Compute M := 1 / \max_j |H(j,j)|. Then there can exist no relation vector whose Euclidean norm is less than M. Upon completion, the desired relation is found in the column of B corresponding to the zero entry of y.

  9. Decrease of Smallest Entry of the Vector y = x * Bn in a Typical PSLQ Run x axis: n (iteration number) y axis: log10 (mink |y(k)|)

  10. Typical Scenario for Using PSLQ in Experimental Mathematics • Compute the numerical value of some constant alpha to high precision (typically several hundred digits). • Examples of alpha: limit of sequence, sum of infinite series, value of definite integral, result of minimization or maximization, etc. • Conjecture, based on experience with related problems, the terms ti of an analytic formula that this constant might satisfy (with unknown integer or rational coefficients). • Compute the value of each of these terms to the same precision as alpha. • Form the vector X = (alpha, t1, t2, …, tn) as input to an integer relation finding algorithm. • If the algorithm finds an integer relation for X, then a formula for alpha can be found by simply solving the relation for alpha. • Even if no relation is found, PSLQ produces useful bounds on the size of any possible relation involving the given terms.

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

  12. Bifurcation Points in Chaos Theory Let B3 be the smallest r such that the “logistic iteration” exhibits 8-way periodicity instead of 4-way periodicity. By means of a sequential approximation scheme, one can obtain the numerical value B3 = 3.54409035955… to any desired precision.

  13. PSLQ and Logistic Iteration Bifurcation Points By computing the vector (1, t, t2, t3, …, t12), where t = B3, to 200-digit precision, then applying PSLQ, one can find that B3 is a root of Similarly define B4 as the smallest r such that the logistic iteration exhibits 16-way periodicity instead of 8-way periodicity. One can compute B4 = 3.564407266095… to arbitrarily high precision as before. Recently B4(2-B4) was identified as the root of a 128-degree polynomial by a much more challenging computation. Thus B4 is a root of a 256-degree polynomial. These results have subsequently been proven formally. 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.

  14. PSLQ and Sculpture The complement of the figure-eight knot, when viewed in hyperbolic space, has finite volume that can be computed to be V = 2.029883212819307250042… Recently David Broadhurst found that This was done by computing V/sqrt(3) and the six terms Tk (k = 0 to 5) to high precision, then applying PSLQ. Jonathan M. Borwein and David H. Bailey, Mathematics by Experiment, AK Peters, 2004, pg. 53-56.

  15. Some Supercomputer-Class PSLQ Computations • Identification of B4, the fourth bifurcation point of the logistic iteration: • Integer relation of size 121; 10,000-digit arithmetic. • Identification of Apery sums: • 15 integer relation problems, with size up to 118, requiring up to 5,000-digit arithmetic. • Discovery of recursions in the rows of Ising Integrals: • 2660 individual PSLQ runs, each of size 36; 1200-digit arithmetic. • Identification of Euler-zeta sums: • Hundreds of integer relation problems, each of size 145 and requiring 5,000-digit arithmetic. • Run on IBM SP parallel system. • Finding a relation involving root of a Lehmer polynomial: • Integer relation of size 125; 50,000 digit arithmetic. • Utilizes 3-level, multi-pair parallel PSLQ program. • Run on IBM SP using ARPEC; 16 hours on 64 CPUs. Papers by DHB, Jonathan Borwein, David Bradley, David Broadhurst, Richard Crandall and Roland Girgensohn.

  16. BBP-Type Formulas Until recently, it was widely believed that it is impossible to compute digits of various well-known mathematical constants beginning at an arbitrary position n, except by computing all digits up through position n. In 1996, Peter Borwein and Simon Plouffe observed that a well-known formula permitted this to be done for the binary digits of log(2): They immediately began to seek other constants with formulas of the type where p(n) and q(n) are integer polynomials, with q nonzero at nonnegative integer arguments. Any rational-linear combination of such constants also has this property.

  17. PSLQ and the BBP Formula for Pi Simon Plouffe used DHB’s PSLQ program, with pi, log(2) and about 25 other BBP-type constants as input (all computed to 200-digit precision), in an attempt to see if pi satisfies a BBP-type formula. After several attempts, the program discovered this formula: Indeed, this formula permits one to compute binary (or hexadecimal) digits of pi beginning at an arbitrary starting position. The largest such computation to date (by Colin Percival of Simon Fraser University) found approximately 100 binary digits of pi beginning at the quadrillionth (1015) position. David H. Bailey, Peter B. Borwein and Simon Plouffe, “On the Rapid Computation of Various Polylogarithmic Constants,” Mathematics of Computation, vol. 66, no. 218 (Apr 1997), pg. 903-913.

  18. Some Other BBP-Type Identities Papers by DHB, Peter B. Borwein, Simon Plouffe, David Broadhurst and Richard Crandall.

  19. Is There a Base-10 Formula for Pi? For some constants, both a base-2 and a base-3 formula are known. Question: Is there any base-n BBP-type formula for pi for nonbinary n? Answer: A 2004 paper by Borwein, Galway and Borwein ruled out the existence of any degree-1 BBP formula for pi in a non-binary base. This does not rule out some completely different scheme for computing non-binary digits of pi beginning at an arbitrary starting point. 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.

  20. A Connection Between BBP Formulas and 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. Let {} denote fractional part. Consider the sequence defined by x0 = 0, 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:pi is 16-normal if and only if this sequence is equidistributed in the unit interval. A similar result holds for any constant with a BBP-type formula. David H. Bailey and Richard E. Crandall, "On the Random Character of Fundamental Constant Expansions," Experimental Mathematics, vol. 10, no. 2 (Jun 2001), pg. 175-190.

  21. Richardson’s Algorithm for Sequence Limit Extrapolation Richardson’s algorithm is one of the most widely used techniques for extrapolating the limit of a sequence to high precision. Algorithm: Pick r (we typically set r = 2 or r = 4). Given a sequence (xm), for each m > 0 set Am,1 = xm, and then for k = 2 to k = m, successively set This recursive scheme generates a triangular matrix A. The best estimates for the limit of xm are the diagonal values Am,m. Avram Sidi, Practical Extrapolation Methods: Theory and Applications, Cambridge University Press, 2002, pg. 21-41.

  22. Application: The Quinn-Rand-Strogatz Constant of Nonlinear Physics Quinn, Rand, and Strogatz recently described a nonlinear Winfree-oscillator mean-field system, by means of the formula For large N, s = 1 – c / N, where c = 0.605443657... What is this constant? By means of a Richardson extrapolation scheme, implemented on 64-CPUs of a highly parallel computer system, we computed c = 0.6054436571967327494789228424472074752208996… This led to the conclusion that c is the root of: where the Hurwicz zeta function is defined as D. H. Bailey, J. M. Borwein and R. E. Crandall, “Resolution of the Quinn-Rand-Strogatz Constant of Nonlinear Physics,” manuscript, Jun 2007, http://crd.lbl.gov/~dhbailey/dhbpapers/QRS.pdf.

  23. The Euler-Maclaurin Formula and Infinite Series Summation The Euler-Maclaurin summation formula approximates a finite sum as an integral with high-order corrections: [Here h = (b - a)/n and xj = a + j h. Dm f(x) means m-th derivative of f.] One can use the E-M to compute a high-precision sum for an infinite series: Explicitly compute, to high precision, the sum of the first N terms of the series, where N = 10p (we typically set p = 8, so that N = 100,000,000). Then use the E-M formula to calculate a high-precision value for the “tail.” Each term of the E-M formula adds roughly p more correct digits. A related technique, “Boole summation,” can be used for alternating series.

  24. Euler’s Transformation for Summing Alternating Infinite Series For example, Catalan’s constant can be computed to 500-digit precision by setting n = 1000, then evaluating 400 terms of the second series (a total of 1400 function evaluations). William H. Press, et al, Numerical Recipes, Cambridge University Press, 1966, pg. 133-134.

  25. Converting All-Positive Series to Alternating Series Given an all-positive series (xn), one can construct an alternating series (yn) with the same sum as follows: Set y0 = x0, then for n > 0 Each of these individual summations converges quite rapidly, so only a modest number of terms typically need to be computed. Euler’s transformation can then be applied to find the sum This method works fairly well, but is many times more costly than the alternating series case. Is there an efficient, general-purpose, numerically robust scheme for finding high-precision values for infinite series sums?

  26. Application: Euler Sum Identities and many other specific and general formulas. These were found by computing high-precision values of the left-hand side (via the E-M formula), speculating what terms might be involved in the right-hand side, then applying PSLQ to find the coefficients. David H. Bailey, Jonathan M. Borwein and Roland Girgensohn, "Experimental Evaluation of Euler Sums," Experimental Mathematics, vol. 3, no. 1 (1994), pg. 17-30.

  27. Apery-Like Sum Identities These were found using PSLQ and straightforward numerical summation: D. H. Bailey, J. M. Borwein and D. M. Bradley, “Experimental Determination of Apery-Like Identities for Zeta(2n+2),” Experimental Mathematics, vol. 15 (2006), pg. 281-289 .

  28. General Apery-Like Sum Identities: Discovered and Proven by Computer Based on specific results obtained via a “bootstrap” process, we found the following general formulas and then proved them using the Wilf-Zeilberger algorithm: PSLQ: Great for discovering new identities, but not helpful for proof. Wilf-Zeilberger: Great for proving identities, but not helpful for discovery. Together they are a great combination!

  29. Ramanujan-Like Identities Guillera recently found some Ramanujan-like identities, including where Guillera proved the first two of these using the Wilf-Zeilberger algorithm. He ascribed the third to Gourevich, who found it using integer relation methods. Are there any higher-order analogues?

  30. Relations Found by PSLQ(in addition to Guillera’s three relations) David H. Bailey and Jonathan M. Borwein, "Computer-Assisted Discovery and Proof," Feb 2007, http://crd.lbl.gov/~dhbailey/dhbpapers/comp-disc-proof.pdf.

  31. Gaussian Quadrature Gaussian quadrature approximates the integral of a function on [-1, 1] as Here the abscissas xj are the roots of the Legendre polynomial Pn(x) and the weights wj are given by Values of Pn(x) can be computed using the recurrence The abscissas and weights can be pre-computed. Gaussian quadrature is the best integration scheme for continuous, well-behaved functions at moderate levels of precision (< 1000 digits).

  32. Numerical Integration and the Euler-Maclaurin Formula [Here h = (b - a)/n and xj = a + j h. Dm f(x) means m-th derivative of f.] Note when f(t) and all of its derivatives are zero at a and b (as in a bell-shaped curve), the error E(h) of a simple trapezoidal approximation to the integral goes to zero more rapidly than any power of h. Kendall Atkinson, An Introduction to Numerical Analysis, John Wiley, 1989, pg. 289.

  33. Trapezoidal Approximation to a Bell-Shaped Function

  34. Tanh-Sinh Quadrature Given f(x) defined on (-1,1), define g(t) = tanh (pi/2 sinh t). Then setting x = g(t) yields where xj = g(hj) and wj = g’(hj). Since g’(t) goes to zero very rapidly for large t, then even if f(x) has a vertical derivative or blow-up singularity at an endpoint, the product f(g(t)) g’(t) typically is a nice bell-shaped function for which the E-M formula applies. Reducing h by half typically doubles the number of correct digits. Tanh-sinh quadrature is the best integration scheme for functions with vertical derivatives or blow-up singularities at endpoints, or for any function at very high precision (> 1000 digits). 1. 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. 2. H. Takahasi and M. Mori, “Double Exponential Formulas for Numerical Integration,” Publications of RIMS, Kyoto University, vol. 9 (1974), pg. 721–741.

  35. Original and Transformed Integrand Functions Original function on [-1,1] with singularities at endpoints: Transformed function on real line using tanh-sinh substitution:

  36. A Tetrahedral Integral Identity This conjectured identity arises from analysis of volumes of ideal tetrahedra in hyperbolic space. We have verified this numerically to 20,000 digits (using highly parallel tanh-sinh quadrature), but no formal proof is known. D.H. Bailey, J.M. Borwein, V. Kapoor and E. Weisstein, “Ten Problems in Experimental Mathematics,” American Mathematical Monthly, vol. 113, no. 6 (Jun 2006), pg. 481-409 .

  37. Parallel Run Times for the 20,000-Digit Tetrahedral Integral Computation 1-CPU timings are sums of timings from a 64-CPU run, where barrier waits and communication were not timed. We believe this is the largest single numerical integration ever done, parallel or serial.

  38. Box Integrals In 2006, Luis Goddyn of SFU suggested we examine integrals of the form which can be thought of as the average distance from the origin to the sides of an n-dimensional hypercube. The following evaluations are now known: where D. H. Bailey, J. M. Borwein and R. E. Crandall, “Box Integrals,” Journal of Computational and Applied Mathematics, vol. 206 (2007), pg. 196-208.

  39. Ising Integrals We recently applied our methods to study three classes of integrals that arise in the Ising theory of mathematical physics: David H. Bailey, Jonathan M. Borwein and Richard E. Crandall, “Integrals of the Ising Class,” Journal of Physics A: Mathematical and General, vol. 39 (2006), pg. 12271-12302.

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

  41. Limiting Value of Cn The Cn numerical values approach a limit: What is this limit? We copied the first 50 digits of this numerical value into the online Inverse Symbolic Calculator tool, available at http://oldweb.cecm.sfu.ca/projects/ISC/ISCmain.html The result was: where gamma denotes Euler’s constant. Note: An improved and updated version of the ISC (2.0) is now available at: http://ddrive.cs.dal.ca/~isc

  42. Other Ising Integral Evaluations

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

  44. Recursions in Ising Integrals Consider the 2-parameter class of Ising integrals After computing 1000-digit numerical values for all n up to 36 and all k up to 75 (a total of 2660 individual quadrature calculations, performed on a highly parallel computer system), we discovered (using PSLQ) linear relations in the rows of this array. For example, when n = 3: These recursions have been proven for n = 1, 2, 3, 4. Similar, but more complicated, recursions have been found for larger n (see next page). David H. Bailey, David Borwein, Jonathan M. Borwein and Richard Crandall, “Hypergeometric Forms for Ising-Class Integrals,” Experimental Mathematics, to appear, http://crd.lbl.gov/~dhbailey/dhbpapers/meijer/pdf.

  45. Experimental Recursion for n = 24

  46. General Recursion Formulas We were then able to find general recursion formulas for each n up to 36: New Result: Jonathan Borwein and Bruno Salvy have given an explicit form for these recursions, together with code to compute any desired case. Jonathan M. Borwein and Bruno Salvy, “A Proof of a Recursion for Bessel Moments,” manuscript, 2007, http://users.cs.dal.ca/~jborwein/recursion.pdf.

  47. Heisenberg Spin Integrals In another application of experimental high-precision integration to mathematical physics, we recently investigated some integrals first studied by Boos and Korepin: Here C denotes the contour {x - i/2, x on real line), and D. H. Bailey, J. M. Borwein, R. E. Crandall, D. Manna, “New Representations for Spin Integrals,” manuscript, 2007 (work in progress).

  48. Transformation of Spin Integrals We were able to transform this expression to the following more manageable form over a finite n-dimensional interval: By evaluating this n-dimensional integral numerically, we have verified some analytic evaluations given by Boos and Korepin, and hope to extend their results.

  49. Evaluations of P(n)Derived Analytically, Confirmed Numerically

  50. Evaluation of P(6) where

More Related