MathWorks 2018 March

Presentation at MathWorks as part of exchange of ideas

Transcript

1. Adaptive Algorithms: Flaws and Fixes Fred J. Hickernell Department of Applied Mathematics Center for Interdisciplinary Scientific Computation Illinois Institute of Technology hickernell@iit.edu mypages.iit.edu/~hickernell Thank you for hosting me Thanks to the GAIL team Supported by NSF-DMS-1522687 Visit to MathWorks, March 14-15, 2018
2. 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, design of experiments MATLAB user for thirty years, contributed to MATLAB’s sobolset 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/15
3. 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, design of experiments MATLAB user for thirty years, contributed to MATLAB’s sobolset Developer of the Guaranteed Automatic Integration Library (GAIL) , designed to make research code accessible and production code reliable What I wish to accomplish Technical discussions regarding the components of MATLAB that I do or should use Learn from your perspective on producing robust, accurate numerical software Share my perspective on adaptive numerical software, like integral, fminbnd, and Chebfun, but with accuracy guarantees Explore how, as a developer and center director, I can better partner with MATLAB 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/15
4. Introduction Flaws Fixes Examples GAIL Development Closing References Adaptive Algorithms sol(f) =                Rd f(x) (x) dx Bayesian inference, financial risk, statistical physics, ... f surrogate modeling for computer experiments, ... min x∈X f(x) process improvement, ... . . . 3/15
5. Introduction Flaws Fixes Examples GAIL Development Closing References Adaptive Algorithms sol(f) =                Rd f(x) (x) dx Bayesian inference, financial risk, statistical physics, ... f surrogate modeling for computer experiments, ... min x∈X f(x) process improvement, ... . . . soln (f) = n i=1 wi,n f(xi ) =            integral Chebfun fminbnd . . . Theory: sol(f) − soln (f) C(n) f 3/15
6. Introduction Flaws Fixes Examples GAIL Development Closing References Adaptive Algorithms sol(f) =                Rd f(x) (x) dx Bayesian inference, financial risk, statistical physics, ... f surrogate modeling for computer experiments, ... min x∈X f(x) process improvement, ... . . . soln (f) = n i=1 wi,n f(xi ) =            integral Chebfun fminbnd . . . Theory: sol(f) − soln (f) C(n) f Want app(·, ·) based on some sequence of soln , with sol(f) − app(f, ε) ε, ∀ε > 0, reasonable f 3/15
7. Introduction Flaws Fixes Examples GAIL Development Closing References Adaptive Algorithms sol(f) =                Rd f(x) (x) dx Bayesian inference, financial risk, statistical physics, ... f surrogate modeling for computer experiments, ... min x∈X f(x) process improvement, ... . . . soln (f) = n i=1 wi,n f(xi ) =            integral Chebfun fminbnd . . . Theory: sol(f) − soln (f) C(n) f Need: sol(f) − soln (f) err {f(xi )}n i=1 Want app(·, ·) based on some sequence of soln , with sol(f) − app(f, ε) ε, ∀ε > 0, reasonable f When n is large enough to make err {f(xi )}n i=1 ε, then set app(f, ε) = soln (f) 3/15
8. 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 in GAIL 4/15
9. 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). 5/15
10. 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. 6/15
11. 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 6/15
12. 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 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). 6/15
13. 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. 6/15
14. 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. 6/15
15. 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+. 6/15
16. 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 7/15
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). 8/15
18. 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. 9/15
19. 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/. 10/15
20. 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+). 11/15
21. Introduction Flaws Fixes Examples GAIL Development Closing References GAIL Algorithms vs. MathWorks/Chebfun Adaptive Algorithms MathWorks/Chebfun GAIL integral/sum no guarantee integral_g guarantee high order low order (trapezoidal rule) locally/globally adaptive globally adaptive ??/chebfun no guarantee funappx_g guarantee ??/high order low order (linear spline) ??/globally adaptive locally adaptive fminbnd/min no guarantee funmin_g guarantee local/global minimum global minimum low/high order low order (linear spline) locally/globally adaptive locally adaptive integral2/sum no guarantee meanMC_g guarantee integral3/sum dimension 3 cubSobol_g arbitrary dimension locally adaptive cubLattice_g globally adaptive cubBayesLattice 12/15
22. Introduction Flaws Fixes Examples GAIL Development Closing References GAIL Goals and How MathWorks Might Help Automatic, adaptive algorithms with theoretical underpinnings for the stopping criteria Published theorems that form the foundation 13/15
23. Introduction Flaws Fixes Examples GAIL Development Closing References GAIL Goals and How MathWorks Might Help 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 We need custom fast Fourier and Walsh transforms for Sobol’, lattice and Bayesian cubature error because our data are in non-standard order; but our transforms written in MATLAB are much slower than the built in ones; can this be fixed? 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 MATLAB? 13/15
24. Introduction Flaws Fixes Examples GAIL Development Closing References GAIL Goals and How MathWorks Might Help 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 MATLAB’s Sobol’ generator is not 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). 13/15
25. Introduction Flaws Fixes Examples GAIL Development Closing References GAIL Goals and How MathWorks Might Help 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; are we doing it right? 13/15
26. Introduction Flaws Fixes Examples GAIL Development Closing References GAIL Goals and How MathWorks Might Help 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 Generating documentation is clumsy; hard to re-use help xxx information for doc xxx; Chebfun documentation is nice Unit tests and doctests Would appreciate more guidance from MathWorks 13/15
27. Introduction Flaws Fixes Examples GAIL Development Closing References GAIL Goals and How MathWorks Might Help 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 MathWorks better support the volunteer developer community to create a robust collection of MATLAB 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/. 13/15
28. 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 Is there a good update of or replacement for Moler’s book ? 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 MathWorks? Moler, C. Numerical Computing with MATLAB. (SIAM, 2004). 14/15
29. Thank you I am grateful for your hospitality and our interactions Slides available at speakerdeck.com/fjhickernell/mathworks-2018-march
30. 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). 15/15
31. 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). 15/15
32. 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/. 15/15
33. 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/. Moler, C. Numerical Computing with MATLAB. (SIAM, 2004). (eds Cools, R. & Nuyens, D.) Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014. 163 (Springer-Verlag, Berlin, 2016). 15/15