Adaptive Algorithms: Flaws and Fixes Fred J. Hickernell Department of Applied Mathematics Center for Interdisciplinary Scientific Computation Illinois Institute of Technology [email protected] mypages.iit.edu/~hickernell Thank you for hosting me Thanks to the GAIL team Supported by NSF-DMS-1522687 Visit to NAG, April 10, 2018
Introduction Flaws Fixes Examples GAIL Development Closing References Personal Introduction Who I am professionally PhD in 1981 from MIT, professor at University of Southern California, Hong Kong Baptist University, and Illinois Institute of Technology Interested in computational math and statistics: adaptive algorithms, (quasi-)Monte Carlo, kriging (Gaussian process modeling, meshfree methods), design of experiments Contributed to NAG’s nag_quasi_init_scrambled (g05ync) Developer of the Guaranteed Automatic Integration Library (GAIL) , designed to make research code accessible and production code reliable Hong, H. S. & H., F. J. Algorithm 823: Implementing Scrambled Digital Nets. ACM Trans. Math. Software 29, 95–109 (2003). Choi, S.-C. T. et al. GAIL: Guaranteed Automatic Integration Library (Versions 1.0–2.2). MATLAB software. 2013–2017. http://gailgithub.github.io/GAIL_Dev/. 2/17
Introduction Flaws Fixes Examples GAIL Development Closing References Personal Introduction Who I am professionally PhD in 1981 from MIT, professor at University of Southern California, Hong Kong Baptist University, and Illinois Institute of Technology Interested in computational math and statistics: adaptive algorithms, (quasi-)Monte Carlo, kriging (Gaussian process modeling, meshfree methods), design of experiments Contributed to NAG’s nag_quasi_init_scrambled (g05ync) Developer of the Guaranteed Automatic Integration Library (GAIL) , designed to make research code accessible and production code reliable What I wish to accomplish Learn from your perspective on producing robust, accurate numerical software Share my perspective on adaptive numerical software with accuracy guarantees Explore how, as a developer and center director, I can better partner with NAG Hong, H. S. & H., F. J. Algorithm 823: Implementing Scrambled Digital Nets. ACM Trans. Math. Software 29, 95–109 (2003). Choi, S.-C. T. et al. GAIL: Guaranteed Automatic Integration Library (Versions 1.0–2.2). MATLAB software. 2013–2017. http://gailgithub.github.io/GAIL_Dev/. 2/17
Introduction Flaws Fixes Examples GAIL Development Closing References Good Theory Theorem Template If f ∈ C, where C is a cone of reasonable input functions, then the algorithm app(·, ·) satisfies sol(f) − app(f, ε) ε ∀ε > 0, f ∈ C. The computational cost of app(f, ε) is O(( f F /ε)1/p), and this cost is asymptotically optimal. 4/17
Introduction Flaws Fixes Examples GAIL Development Closing References Good Theory Theorem Template If f ∈ C, where C is a cone of reasonable input functions, then the algorithm app(·, ·) satisfies sol(f) − app(f, ε) ε ∀ε > 0, f ∈ C. The computational cost of app(f, ε) is O(( f F /ε)1/p), and this cost is asymptotically optimal. Some say: Check that your problem satisfies the hypothesis, then you can be confident in the conclusion 4/17
Introduction Flaws Fixes Examples GAIL Development Closing References Good Theory Theorem Template If f ∈ C, where C is a cone of reasonable input functions, then the algorithm app(·, ·) satisfies sol(f) − app(f, ε) ε ∀ε > 0, f ∈ C. The computational cost of app(f, ε) is O(( f F /ε)1/p), and this cost is asymptotically optimal. Some say: Check that your problem satisfies the hypothesis, then you can be confident in the conclusion I say: Practically, you cannot check that your problem satisfies the conclusion 4/17
Introduction Flaws Fixes Examples GAIL Development Closing References Good Theory Theorem Template If f ∈ C, where C is a cone of reasonable input functions, then the algorithm app(·, ·) satisfies sol(f) − app(f, ε) ε ∀ε > 0, f ∈ C. The computational cost of app(f, ε) is O(( f F /ε)1/p), and this cost is asymptotically optimal. Some say: Check that your problem satisfies the hypothesis, then you can be confident in the conclusion I say: Practically, you cannot check that your problem satisfies the conclusion Theory tells you what to expect for reasonable problems 4/17
Introduction Flaws Fixes Examples GAIL Development Closing References Good Theory Theorem Template If f ∈ C, where C is a cone of reasonable input functions, then the algorithm app(·, ·) satisfies sol(f) − app(f, ε) ε ∀ε > 0, f ∈ C. The computational cost of app(f, ε) is O(( f F /ε)1/p), and this cost is asymptotically optimal. Some say: Check that your problem satisfies the hypothesis, then you can be confident in the conclusion I say: Practically, you cannot check that your problem satisfies the conclusion Theory tells you what to expect for reasonable problems If you observe better than expected, your theory needs to be strengthened 4/17
Introduction Flaws Fixes Examples GAIL Development Closing References Good Theory Theorem Template If f ∈ C, where C is a cone of reasonable input functions, then the algorithm app(·, ·) satisfies sol(f) − app(f, ε) ε ∀ε > 0, f ∈ C. The computational cost of app(f, ε) is O(( f F /ε)1/p), and this cost is asymptotically optimal. Some say: Check that your problem satisfies the hypothesis, then you can be confident in the conclusion I say: Practically, you cannot check that your problem satisfies the conclusion Theory tells you what to expect for reasonable problems If you observe better than expected, your theory needs to be strengthened If you observe worse than expected, theory tells you what went wrong 4/17
Introduction Flaws Fixes Examples GAIL Development Closing References Good Theory Theorem Template If f ∈ C, where C is a cone of reasonable input functions, then the algorithm app(·, ·) satisfies sol(f) − app(f, ε) ε ∀ε > 0, f ∈ C. The computational cost of app(f, ε) is O(( f F /ε)1/p), and this cost is asymptotically optimal. Some say: Check that your problem satisfies the hypothesis, then you can be confident in the conclusion I say: Practically, you cannot check that your problem satisfies the conclusion Theory tells you what to expect for reasonable problems If you observe better than expected, your theory needs to be strengthened If you observe worse than expected, theory tells you what went wrong Don’t cheat by defining C as the cone for which the conclusion holds 4/17
Introduction Flaws Fixes Examples GAIL Development Closing References Non-Adaptive Trapezoidal Rule Theory Meets Practice sol(f) = 1 0 f(x) dx, Tn (f) = 1 n f(0) 2 + f 1 n + · · · + f n − 1 n + f(1) 2 |sol(f) − Tn (f)| f ∞ 12n2 As expected 5/17
Introduction Flaws Fixes Examples GAIL Development Closing References Non-Adaptive Trapezoidal Rule Theory Meets Practice sol(f) = 1 0 f(x) dx, Tn (f) = 1 n f(0) 2 + f 1 n + · · · + f n − 1 n + f(1) 2 |sol(f) − Tn (f)| f ∞ 12n2 Better than expected Integrand not smooth enough! 5/17
Introduction Flaws Fixes Examples GAIL Development Closing References Spiky Functions Chebfun misses spike because it samples too sparsely by default fminbnd misses spike because it assumes a convex function Can we set defaults to catch spikes? Yes, up to a point, in GAIL 6/17
Introduction Flaws Fixes Examples GAIL Development Closing References Fluky Functions 1 0 f(x) dx = 0.27882695 . . . integral(f,0,1,'AbsTol',1e−13, ... 'RelTol',1e−13) = 0.27879907 f is not spiky relative to the number of initial samples by integral integral can fail for integrands that are not even spiky because the “error bound” is the difference between two quadrature rules Can we get fluky functions right? Yes in GAIL Lyness, J. N. When Not to Use an Automatic Quadrature Routine. SIAM Rev. 25, 63–87 (1983). 7/17
Introduction Flaws Fixes Examples GAIL Development Closing References The GAIL Approach (We do not use interval arithmetic ) Want app(·, ·) that guarantees sol(f) − app(f, ε) ε for all ε > 0, f ∈ set of reasonable functions C Have sol(f) − soln (f) ∀f∈F C(n) f F Construct err {f(xi )}n i=1 ∀f∈C C(n) f F , where C is a cone ←− challenge Increase n until err {f(xi )}n i=1 ε and let app(f, ε) = soln (f) Rump, S. M. INTLAB - INTerval LABoratory. in Developments in Reliable Computing (ed Csendes, T.) URL: http://www.ti3.tuhh.de/rump/ (Kluwer Academic Publishers, Dordrecht, 1999), 77–104. 8/17
Introduction Flaws Fixes Examples GAIL Development Closing References The GAIL Approach Construct err {f(xi )}n i=1 ∀f∈C C(n) f F , where C is a cone ←− challenge Don’t cop out by defining C as the cone for which the inequality is satisfied 8/17
Introduction Flaws Fixes Examples GAIL Development Closing References The GAIL Approach Construct err {f(xi )}n i=1 ∀f∈C C(n) f F , where C is a cone ←− challenge Don’t cop out by defining C as the cone for which the inequality is satisfied For univariate integration integral_g, approximation funappx_g, and optimization funmin_g, bound f F by inflated finite difference approximations , e.g., Var(f ) C 2 n n−1 i=1 f(ti+1 ) − 2f(ti ) + f(ti−1 ) n , ti = i n , i = 0, . . . , n, C(h) = C(0)h h − h h says how wide ignored spikes can be, C(0) says how much Var(f ) can deviate from the finite difference approximation as n → ∞ Clancy, N. et al. The Cost of Deterministic, Adaptive, Automatic Algorithms: Cones, Not Balls. J. Complexity 30, 21–45 (2014), Choi, S.-C. T. et al. Local Adaption for Approximation and Minimization of Univariate Functions. J. Complexity 40, 17–33 (2017). 8/17
Introduction Flaws Fixes Examples GAIL Development Closing References The GAIL Approach Construct err {f(xi )}n i=1 ∀f∈C C(n) f F , where C is a cone ←− challenge Don’t cop out by defining C as the cone for which the inequality is satisfied For Monte Carlo cubature meanMC_g approximate var(f(X)) by an inflated sample variance assuming a bound on the kurtosis H., F. J. et al. Guaranteed Conservative Fixed Width Confidence Intervals Via Monte Carlo Sampling. in Monte Carlo and Quasi-Monte Carlo Methods 2012 (eds Dick, J. et al.) 65 (Springer-Verlag, Berlin, 2013), 105–128. 8/17
Introduction Flaws Fixes Examples GAIL Development Closing References The GAIL Approach Construct err {f(xi )}n i=1 ∀f∈C C(n) f F , where C is a cone ←− challenge Don’t cop out by defining C as the cone for which the inequality is satisfied For quasi-Monte Carlo cubature cubSobol_g and cubLattice_g track the decay of the discrete Fourier coefficients of f H., F. J. & Jiménez Rugama, Ll. A. Reliable Adaptive Cubature Using Digital Sequences. in Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014 (eds Cools, R. & Nuyens, D.) 163. arXiv:1410.8615 [math.NA] (Springer-Verlag, Berlin, 2016), 367–383, Jiménez Rugama, Ll. A. & H., F. J. Adaptive Multidimensional Integration Based on Rank-1 Lattices. in Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014 (eds Cools, R. & Nuyens, D.) 163. arXiv:1411.1966 (Springer-Verlag, Berlin, 2016), 407–422. 8/17
Introduction Flaws Fixes Examples GAIL Development Closing References The GAIL Approach Construct err {f(xi )}n i=1 ∀f∈C C(n) f F , where C is a cone ←− challenge Don’t cop out by defining C as the cone for which the inequality is satisfied For Bayesian cubature cubBayesianLattice construct the credible interval for the integral assuming that f comes from a Gaussian process fitted by maximum likelihood estimation Rathinavel, J. & H., F. J. Automatic Bayesian Cubature. in preparation. 2018+. 8/17
Introduction Flaws Fixes Examples GAIL Development Closing References Spiky Functions Chebfun misses spike because it samples too sparsely by default fminbnd misses spike because it assumes a convex function funappx_g and funmin_g catch the spike because they sample initially at enough sites and they detect the high curvature 9/17
Introduction Flaws Fixes Examples GAIL Development Closing References Fluky Functions 1 0 f(x) dx = 0.27882695 . . . integral(f,0,1,'AbsTol',1e−13, ... 'RelTol',1e−13) = 0.27879907 integral_g(f,0,1,'AbsTol',1e−8) ... = 0.27882695 integral can fail for integrands that are not even spiky because the “error bound” is the difference between two quadrature rules integral_g succeeds by sensing the curvature and sampling accordingly Lyness, J. N. When Not to Use an Automatic Quadrature Routine. SIAM Rev. 25, 63–87 (1983). 10/17
Introduction Flaws Fixes Examples GAIL Development Closing References Error of Monte Carlo Methods error = Rd f(x) (x) dx − 1 n n i=1 f(xi ) = CNF(f, {xi }n i=1 ) DSC({xi }n i=1 ) VAR(f) VAR(f) 0 and measures the variation of f: Can be decreased using importance sampling or control variates DSC({xi }n i=1 ) → 0 as n → ∞ and measures the discrepancy sampling distribution deviates from the distribution defining the integral: O(n−1/2) for IID xi O(n−1+ ) or better for low discrepancy xi O(n−r) for smooth f and unequal weights for the f(xi ) CNF(f, {xi }n i=1 ) = O(1) and measures how confounded f is with the difference between the sampling and true distributions H., F. J. The Trio Identity for Quasi-Monte Carlo Error Analysis. in Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Stanford, USA, August 2016 (eds Glynn, P. & Owen, A.) to appear, arXiv:1702.01487 (Springer-Verlag, Berlin, 2018), 13–37. 11/17
Introduction Flaws Fixes Examples GAIL Development Closing References Low Discrepancy Sampling Low discrepancy sampling places the xi more evenly than IID sampling IID points Sobol’ points Integration lattice points ··· Dick, J. et al. High dimensional integration — the Quasi-Monte Carlo way. Acta Numer. 22, 133–288 (2013), H., F. J. et al. SAMSI Program on Quasi-Monte Carlo and High-Dimensional Sampling Methods for Applied Mathematics. https://www.samsi.info/programs-and-activities/year-long-research-programs/2017-18-program-quasi- monte-carlo-high-dimensional-sampling-methods-applied-mathematics-qmc/. 12/17
Introduction Flaws Fixes Examples GAIL Development Closing References Option Pricing µ = fair price = Rd e−rT max 1 d d j=1 Sj − K, 0 e−zTz/2 (2π)d/2 dz ≈ $13.12, T = 1/4, d = 13 here Sj = S0 e(r−σ2/2)jT/d+σxj = stock price at time jT/d x = Az, AAT = Σ = min(i, j)T/d d i,j=1 Abs. Error Median Worst 10% Worst 10% Tolerance Method A Error Accuracy n Time (s) 1E−2 IID diff 2E−3 100% 6.1E7 33.000 1E−2 Scr. Sobol’ PCA 1E−3 100% 1.6E4 0.040 1E−2 Scr. Sob. cont. var. PCA 2E−3 100% 4.1E3 0.019 1E−2 Bayes. Latt. PCA 2E−3 99% 1.6E4 0.051 H., F. J. et al. Guaranteed Conservative Fixed Width Confidence Intervals Via Monte Carlo Sampling. in Monte Carlo and Quasi-Monte Carlo Methods 2012 (eds Dick, J. et al.) 65 (Springer-Verlag, Berlin, 2013), 105–128, Rathinavel, J. & H., F. J. Automatic Bayesian Cubature. in preparation. 2018+, H., F. J. et al. in Contemporary Computational Mathematics — a celebration of the 80th birthday of Ian Sloan (eds Dick, J. et al.) to appear, arXiv:1702.01491 [math.NA] (Springer-Verlag, 2018+). 13/17
Introduction Flaws Fixes Examples GAIL Development Closing References GAIL Goals and How NAG Might Be Interested Automatic, adaptive algorithms with theoretical underpinnings for the stopping criteria Published theorems that form the foundation 15/17
Introduction Flaws Fixes Examples GAIL Development Closing References GAIL Goals and How NAG Might Be Interested Automatic, adaptive algorithms with theoretical underpinnings for the stopping criteria Numerical efficiency Automatic so minimal wasted computer effort Our Bayesian cubature is O(n log n) instead of O(n3) because our sampling sites match our covariance kernels Univariate integration and function approximation need to be sped up with higher order algorithms Would like to run (quasi-)Monte Carlo algorithms in parallel, but no success yet; NSF & DOE are pushing advanced architectures; what can be done in NAG? 15/17
Introduction Flaws Fixes Examples GAIL Development Closing References GAIL Goals and How NAG Might Be Interested Automatic, adaptive algorithms with theoretical underpinnings for the stopping criteria Numerical efficiency Start with simple algorithms and increase sophistication Hence, need to develop higher order univariate integration and function approximation Want to implement multi-level Sobol’ cubature , but not sure if NAG’s Sobol’ generator is flexible enough Would like higher dimension for Sobol’ generator web.maths.unsw.edu.au/~fkuo/sobol/ Would like higher order nets through interlacing , which achieve higher order rates of convergence Giles, M. Multilevel Monte Carlo Methods. Acta Numer. 24, 259–328 (2015). Joe, S. & Kuo, F. Y. Constructing Sobol sequences with better two-dimensional projections. SIAM J. Sci. Comput. 30, 2635–2654 (2008). Dick, J. Higher order scrambled digital nets achieve the optimal rate of the root mean square error for smooth integrands. Ann. Statist. 39, 1372–1398 (2011). 15/17
Introduction Flaws Fixes Examples GAIL Development Closing References GAIL Goals and How NAG Might Be Interested Automatic, adaptive algorithms with theoretical underpinnings for the stopping criteria Numerical efficiency Start with simple algorithms and increase sophistication Minimal user intervention required Hence, automatic Absolute, relative, and hybrid error tolerances accommodated Algorithm parameters have pre-set defaults Use classes to perform input parsing and validation; what do you do? 15/17
Introduction Flaws Fixes Examples GAIL Development Closing References GAIL Goals and How NAG Might Be Interested Automatic, adaptive algorithms with theoretical underpinnings for the stopping criteria Numerical efficiency Start with simple algorithms and increase sophistication Minimal user intervention required Software tested and documented Unit tests and doctests Hard to keep documentation attractive and comprehensive 15/17
Introduction Flaws Fixes Examples GAIL Development Closing References GAIL Goals and How NAG Might Be Interested Automatic, adaptive algorithms with theoretical underpinnings for the stopping criteria Numerical efficiency Start with simple algorithms and increase sophistication Minimal user intervention required Software tested and documented Software widely available My academic community promotes theorems, software is rare ; this needs to change Want GAIL and related software to be adopted by the community (company, volunteers, or both) in the sense of use and development Available on Github at gailgithub.github.io/GAIL_Dev/ How can NAG better support the developer community to create a robust collection of software with broader capability? Giles, M. https://people.maths.ox.ac.uk/gilesm/mlmc/, L’Ecuyer, P. http://simul.iro.umontreal.ca, Nuyens, D. https://people.cs.kuleuven.be/~dirk.nuyens/. 15/17
Introduction Flaws Fixes Examples GAIL Development Closing References Center for Interdisciplinary Scientific Computation (CISC) Goals Director of newly-founded CISC, want to leverage our strengths across departmental silos and outside of Illinios Tech Want to promote advanced scientific computation in research, development, and education Students see MATLAB, Mathematica, R, etc. in the classroom, but need to see it at the cutting edge of scientific research and software development What mutually beneficial partnerships are there with NAG? 16/17
Introduction Flaws Fixes Examples GAIL Development Closing References Hong, H. S. & H., F. J. Algorithm 823: Implementing Scrambled Digital Nets. ACM Trans. Math. Software 29, 95–109 (2003). Choi, S.-C. T. et al. GAIL: Guaranteed Automatic Integration Library (Versions 1.0–2.2). MATLAB software. 2013–2017. http://gailgithub.github.io/GAIL_Dev/. Lyness, J. N. When Not to Use an Automatic Quadrature Routine. SIAM Rev. 25, 63–87 (1983). Rump, S. M. INTLAB - INTerval LABoratory. in Developments in Reliable Computing (ed Csendes, T.) URL: http://www.ti3.tuhh.de/rump/ (Kluwer Academic Publishers, Dordrecht, 1999), 77–104. Clancy, N., Ding, Y., Hamilton, C., H., F. J. & Zhang, Y. The Cost of Deterministic, Adaptive, Automatic Algorithms: Cones, Not Balls. J. Complexity 30, 21–45 (2014). Choi, S.-C. T., Ding, Y., H., F. J. & Tong, X. Local Adaption for Approximation and Minimization of Univariate Functions. J. Complexity 40, 17–33 (2017). 17/17
Introduction Flaws Fixes Examples GAIL Development Closing References H., F. J., Jiang, L., Liu, Y. & Owen, A. B. Guaranteed Conservative Fixed Width Confidence Intervals Via Monte Carlo Sampling. in Monte Carlo and Quasi-Monte Carlo Methods 2012 (eds Dick, J., Kuo, F. Y., Peters, G. W. & Sloan, I. H.) 65 (Springer-Verlag, Berlin, 2013), 105–128. H., F. J. & Jiménez Rugama, Ll. A. Reliable Adaptive Cubature Using Digital Sequences. in Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014 (eds Cools, R. & Nuyens, D.) 163. arXiv:1410.8615 [math.NA] (Springer-Verlag, Berlin, 2016), 367–383. Jiménez Rugama, Ll. A. & H., F. J. Adaptive Multidimensional Integration Based on Rank-1 Lattices. in Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014 (eds Cools, R. & Nuyens, D.) 163. arXiv:1411.1966 (Springer-Verlag, Berlin, 2016), 407–422. Rathinavel, J. & H., F. J. Automatic Bayesian Cubature. in preparation. 2018+. H., F. J. The Trio Identity for Quasi-Monte Carlo Error Analysis. in Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Stanford, USA, August 2016 (eds Glynn, P. & Owen, A.) to appear, arXiv:1702.01487 (Springer-Verlag, Berlin, 2018), 13–37. Dick, J., Kuo, F. & Sloan, I. H. High dimensional integration — the Quasi-Monte Carlo way. Acta Numer. 22, 133–288 (2013). 17/17
Introduction Flaws Fixes Examples GAIL Development Closing References H., F. J., Kuo, F. Y., L’Ecuyer, P. & Owen, A. B. SAMSI Program on Quasi-Monte Carlo and High-Dimensional Sampling Methods for Applied Mathematics. https://www.samsi.info/programs-and-activities/year-long-research- programs/2017-18-program-quasi-monte-carlo-high-dimensional-sampling- methods-applied-mathematics-qmc/. H., F. J., Jiménez Rugama, Ll. A. & Li, D. in Contemporary Computational Mathematics — a celebration of the 80th birthday of Ian Sloan (eds Dick, J., Kuo, F. Y. & Woźniakowski, H.) to appear, arXiv:1702.01491 [math.NA] (Springer-Verlag, 2018+). Giles, M. Multilevel Monte Carlo Methods. Acta Numer. 24, 259–328 (2015). Joe, S. & Kuo, F. Y. Constructing Sobol sequences with better two-dimensional projections. SIAM J. Sci. Comput. 30, 2635–2654 (2008). Dick, J. Higher order scrambled digital nets achieve the optimal rate of the root mean square error for smooth integrands. Ann. Statist. 39, 1372–1398 (2011). Giles, M. https://people.maths.ox.ac.uk/gilesm/mlmc/. 17/17
Introduction Flaws Fixes Examples GAIL Development Closing References L’Ecuyer, P. http://simul.iro.umontreal.ca. Nuyens, D. https://people.cs.kuleuven.be/~dirk.nuyens/. (eds Cools, R. & Nuyens, D.) Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014. 163 (Springer-Verlag, Berlin, 2016). 17/17
fjhickernell
aya_watanabe
stktu
tatsuropianooo
speed1313
ryou0634
isayan
hpprc
kosuken
shotasasaki
eumesy
rina_glutton
ddemaree
geoffreycrofte
michaelherold
bkeepers
jennybc
jrom
colly
afnizarnur
morganepeng
rmw
danielanewman
dotmariusz