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 NAG, April 10, 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 (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
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
Theorem Template If f ∈ C, where C is a cone of reasonable input functions, then the algorithm app(·, ·) satisﬁes 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
Theorem Template If f ∈ C, where C is a cone of reasonable input functions, then the algorithm app(·, ·) satisﬁes 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 satisﬁes the hypothesis, then you can be conﬁdent in the conclusion 4/17
Theorem Template If f ∈ C, where C is a cone of reasonable input functions, then the algorithm app(·, ·) satisﬁes 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 satisﬁes the hypothesis, then you can be conﬁdent in the conclusion I say: Practically, you cannot check that your problem satisﬁes the conclusion 4/17
Theorem Template If f ∈ C, where C is a cone of reasonable input functions, then the algorithm app(·, ·) satisﬁes 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 satisﬁes the hypothesis, then you can be conﬁdent in the conclusion I say: Practically, you cannot check that your problem satisﬁes the conclusion Theory tells you what to expect for reasonable problems 4/17
Theorem Template If f ∈ C, where C is a cone of reasonable input functions, then the algorithm app(·, ·) satisﬁes 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 satisﬁes the hypothesis, then you can be conﬁdent in the conclusion I say: Practically, you cannot check that your problem satisﬁes 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
Theorem Template If f ∈ C, where C is a cone of reasonable input functions, then the algorithm app(·, ·) satisﬁes 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 satisﬁes the hypothesis, then you can be conﬁdent in the conclusion I say: Practically, you cannot check that your problem satisﬁes 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
Theorem Template If f ∈ C, where C is a cone of reasonable input functions, then the algorithm app(·, ·) satisﬁes 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 satisﬁes the hypothesis, then you can be conﬁdent in the conclusion I say: Practically, you cannot check that your problem satisﬁes 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 deﬁning C as the cone for which the conclusion holds 4/17
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
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)| Var(f ) 8n2 As expected Integrand ﬁts revised theory 5/17
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)| Var(f ) 8n2 Better than expected 5/17
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)| Var(f ) 8n2 and |sol(f) − Tn (f)| f (0) − f (1) 12n2 − f (0) − f (1) 720n4 + · · · Better than expected Integrand ﬁts revised theory 5/17
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)| Var(f ) 8n2 and |sol(f) − Tn (f)| f (0) − f (1) 12n2 − f (0) − f (1) 720n4 + · · · Worse than expected 5/17
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)| Var(f ) 8n2 and |sol(f) − Tn (f)| f (0) − f (1) 12n2 − f (0) − f (1) 720n4 + · · · Worse than expected Integrand not smooth enough! 5/17
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
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). 7/17
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
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 8/17
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 , 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 ﬁnite diﬀerence 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
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. 8/17
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. 8/17
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+. 8/17
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
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). 10/17
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. 11/17
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
µ = 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+). 13/17
and How NAG Might Be Interested Automatic, adaptive algorithms with theoretical underpinnings for the stopping criteria Published theorems that form the foundation 15/17
and How NAG Might Be Interested 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 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
and How NAG Might Be Interested 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 not sure if NAG’s Sobol’ generator is ﬂ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). 15/17
and How NAG Might Be Interested 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; what do you do? 15/17
and How NAG Might Be Interested 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 Unit tests and doctests Hard to keep documentation attractive and comprehensive 15/17
and How NAG Might Be Interested 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 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
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 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 NAG? 16/17
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
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). 17/17
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
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