Applied Mathematics Center for Interdisciplinary Scientiﬁc 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 MathWorks, March 14-15, 2018
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
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
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
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 diﬀerence between two quadrature rules Can we get ﬂuky functions right? Yes in GAIL Lyness, J. N. When Not to Use an Automatic Quadrature Routine. SIAM Rev. 25, 63–87 (1983). 5/15
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
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 deﬁning C as the cone for which the inequality is satisﬁed 6/15
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 deﬁning C as the cone for which the inequality is satisﬁed For univariate integration integral_g, approximation funappx_g, and optimization funmin_g, bound f F by inﬂated ﬁnite diﬀerence 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
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 deﬁning C as the cone for which the inequality is satisﬁed For Monte Carlo cubature meanMC_g approximate var(f(X)) by an inﬂated sample variance assuming a bound on the kurtosis H., F. J. et al. Guaranteed Conservative Fixed Width Conﬁdence 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
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 deﬁning C as the cone for which the inequality is satisﬁed For quasi-Monte Carlo cubature cubSobol_g and cubLattice_g track the decay of the discrete Fourier coeﬃcients 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
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 deﬁning C as the cone for which the inequality is satisﬁed For Bayesian cubature cubBayesianLattice construct the credible interval for the integral assuming that f comes from a Gaussian process ﬁtted by maximum likelihood estimation Rathinavel, J. & H., F. J. Automatic Bayesian Cubature. in preparation. 2018+. 6/15
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
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 diﬀerence 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
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 deﬁning 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 diﬀerence 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
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
µ = 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 diﬀ 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 Conﬁdence 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
and How MathWorks Might Help Automatic, adaptive algorithms with theoretical underpinnings for the stopping criteria Published theorems that form the foundation 13/15
and How MathWorks Might Help Automatic, adaptive algorithms with theoretical underpinnings for the stopping criteria Numerical eﬃciency Automatic so minimal wasted computer eﬀort 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 ﬁxed? 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
and How MathWorks Might Help Automatic, adaptive algorithms with theoretical underpinnings for the stopping criteria Numerical eﬃciency 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 ﬂexible 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
and How MathWorks Might Help Automatic, adaptive algorithms with theoretical underpinnings for the stopping criteria Numerical eﬃciency 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
and How MathWorks Might Help Automatic, adaptive algorithms with theoretical underpinnings for the stopping criteria Numerical eﬃciency 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
and How MathWorks Might Help Automatic, adaptive algorithms with theoretical underpinnings for the stopping criteria Numerical eﬃciency 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
Interdisciplinary Scientiﬁc 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 scientiﬁc 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 scientiﬁc research and software development What mutually beneﬁcial partnerships are there with MathWorks? Moler, C. Numerical Computing with MATLAB. (SIAM, 2004). 14/15
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
J., Jiang, L., Liu, Y. & Owen, A. B. Guaranteed Conservative Fixed Width Conﬁdence 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
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
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