Hickernell Department of Applied Mathematics Center for Interdisciplinary Scientiﬁc Computation Illinois Institute of Technology hickernell@iit.edu mypages.iit.edu/~hickernell Thanks to my friends and collaborators, including the GAIL team and also NSF-DMS-1522687 and NSF-DMS-1638521 (SAMSI) Computational Mathematics and Statistics Seminar, September 13, 2018
It Mean to Solve a Mathematics or Statistics Problem? Problem Solution x2 + 2x − 3 = 0 x = −3, 1 numbers d dx sin(x) cos(x) elementary function 1 √ 2π 1 −∞ e−x2/2 dx Φ(1) special function Φ is the standard Gaussian cumulative distribution function 2/13
It Mean to Solve a Mathematics or Statistics Problem? Problem Solution x2 + 2x − 3 = 0 x = −3, 1 numbers d dx sin(x) cos(x) elementary function 1 √ 2π 1 −∞ e−x2/2 dx Φ(1) special function Φ is the standard Gaussian cumulative distribution function 1 0 f(x) dx Need a numerical solution, i.e., A(f, ε) f known via an algorithm depending on integrand values f(x1), f(x2), . . . i.e., can get f(x) ∀x ∈ [0, 1] such that 1 0 f(x) dx − A(f, ε) ε 2/13
the Trapezoidal Rule? 1 0 f(x) dx ≈ Tn(f) = 1 n f(0) 2 + f 1 n + · · · + f n − 1 n + f(1) 2 Let A(f, ε) = Tn(f) where n = ρ 8ε Then 1 0 f(x) dx − A(f, ε) ε ∀f ∈ F = {f : f 1 ρ}, ε > 0 T50(f) = 0.3957 We say that integration is solvable for integrands in F 3/13
problem is called solvable for a set of inputs F if there exists a numerical algorithm that provides a solution with error no greater than ε. Speciﬁcally, integration is solvable for F if there exists an A for which 1 0 f(x) dx − A(f, ε) ε ∀f ∈ F, ε > 0 Kunsch, R. J., Novak, E. & Rudolf, D. Solvable Integration Problems and Optimal Sample Size Selection. submitted for publication. 2018+. 4/13
problem is called solvable for a set of inputs F if there exists a numerical algorithm that provides a solution with error no greater than ε. Speciﬁcally, integration is solvable for F if there exists an A for which 1 0 f(x) dx − A(f, ε) ε ∀f ∈ F, ε > 0 We demonstrated that integration is solvable for F = {f : f 1 ρ} Kunsch, R. J., Novak, E. & Rudolf, D. Solvable Integration Problems and Optimal Sample Size Selection. submitted for publication. 2018+. 4/13
problem is called solvable for a set of inputs F if there exists a numerical algorithm that provides a solution with error no greater than ε. Speciﬁcally, integration is solvable for F if there exists an A for which 1 0 f(x) dx − A(f, ε) ε ∀f ∈ F, ε > 0 We demonstrated that integration is solvable for F = {f : f 1 ρ} Is integration solvable for F = {f : f 1 < ∞} since 1 0 f(x) dx − Tn(f) f 1 8n2 ? Kunsch, R. J., Novak, E. & Rudolf, D. Solvable Integration Problems and Optimal Sample Size Selection. submitted for publication. 2018+. 4/13
Not Solvable for F = {f : f 1 < ∞} Suppose that there exists A such that 1 0 f(x) dx − A(f, ε) ε ∀f ∈ F = {f : f 1 < ∞}, ε > 0 Observe where A(0, ε) gets function values 5/13
Not Solvable for F = {f : f 1 < ∞} Suppose that there exists A such that 1 0 f(x) dx − A(f, ε) ε ∀f ∈ F = {f : f 1 < ∞}, ε > 0 Observe where A(0, ε) gets function values Construct fbump ∈ F, as tall as you like; note that A(0, ε) = A(fbump , ε) 5/13
Not Solvable for F = {f : f 1 < ∞} Suppose that there exists A such that 1 0 f(x) dx − A(f, ε) ε ∀f ∈ F = {f : f 1 < ∞}, ε > 0 Observe where A(0, ε) gets function values Construct fbump ∈ F, as tall as you like; note that A(0, ε) = A(fbump , ε) Since 1 0 0 dx and 1 0 fbump (x) dx are so diﬀerent, A cannot be successful 5/13
Not Solvable for F = {f : f 1 < ∞} Suppose that there exists A such that 1 0 f(x) dx − A(f, ε) ε ∀f ∈ F = {f : f 1 < ∞}, ε > 0 Observe where A(0, ε) gets function values Construct fbump ∈ F, as tall as you like; note that A(0, ε) = A(fbump , ε) Since 1 0 0 dx and 1 0 fbump (x) dx are so diﬀerent, A cannot be successful Why doesn’t this argument hold for F = {f : f 1 ρ} ? 5/13
Not Solvable for F = {f : f 1 < ∞} Suppose that there exists A such that 1 0 f(x) dx − A(f, ε) ε ∀f ∈ F = {f : f 1 < ∞}, ε > 0 Observe where A(0, ε) gets function values Construct fbump ∈ F, as tall as you like; note that A(0, ε) = A(fbump , ε) Since 1 0 0 dx and 1 0 fbump (x) dx are so diﬀerent, A cannot be successful Why doesn’t this argument hold for F = {f : f 1 ρ} ? Because fbump / ∈ 5/13
Rule the Best Algorithm? Let cost(A, f, ε) denote the number of function values, n required for A(f, ε) and cost(A, F, ε) = sup f∈F cost(A, f, ε) For our trapezoidal rule algorithm, A cost(A, F, ε) = cost(A, f, ε) = ρ 8ε = O ρ ε ∀f ∈ F = {f : f 1 ρ} This is the (up to a constant) the smallest possible cost 6/13
Rule the Best Algorithm? Let cost(A, f, ε) denote the number of function values, n required for A(f, ε) and cost(A, F, ε) = sup f∈F cost(A, f, ε) For our trapezoidal rule algorithm, A cost(A, F, ε) = cost(A, f, ε) = ρ 8ε = O ρ ε ∀f ∈ F = {f : f 1 ρ} This is the (up to a constant) the smallest possible cost for integrands in F = {f : f 1 ρ} , no worse than Simpson’s rule, Gauss rules, etc. 6/13
Trapezoidal Rule Is Optimal for F = {f : f 1 < ρ} Let A be any algorithm satisfying 1 0 f(x) dx − A(f, ε) ε ∀f ∈ F = {f : f 1 < ρ}, ε > 0 Observe where A(0, ε) gets function values 7/13
Trapezoidal Rule Is Optimal for F = {f : f 1 < ρ} Let A be any algorithm satisfying 1 0 f(x) dx − A(f, ε) ε ∀f ∈ F = {f : f 1 < ρ}, ε > 0 Observe where A(0, ε) gets function values Construct fbump ∈ F; note that A(0, ε) = A(fbump , ε) = A(−fbump , ε) 7/13
Trapezoidal Rule Is Optimal for F = {f : f 1 < ρ} Let A be any algorithm satisfying 1 0 f(x) dx − A(f, ε) ε ∀f ∈ F = {f : f 1 < ρ}, ε > 0 Observe where A(0, ε) gets function values Construct fbump ∈ F; note that A(0, ε) = A(fbump , ε) = A(−fbump , ε) For A to succeed for ±fbump , which have diﬀerent integrals, one must have cost(A, F, ε) = O ρ ε the same as the trapezoidal rule algorithm. 7/13
Trapezoidal Rule Is Optimal for F = {f : f 1 < ρ} Let A be any algorithm satisfying 1 0 f(x) dx − A(f, ε) ε ∀f ∈ F = {f : f 1 < ρ}, ε > 0 Observe where A(0, ε) gets function values Construct fbump ∈ F; note that A(0, ε) = A(fbump , ε) = A(−fbump , ε) For A to succeed for ±fbump , which have diﬀerent integrals, one must have cost(A, F, ε) = O ρ ε the same as the trapezoidal rule algorithm. Why do your instructors tell you that other algorithms are better than the trapezoidal rule? 7/13
Trapezoidal Rule Is Optimal for F = {f : f 1 < ρ} Let A be any algorithm satisfying 1 0 f(x) dx − A(f, ε) ε ∀f ∈ F = {f : f 1 < ρ}, ε > 0 Observe where A(0, ε) gets function values Construct fbump ∈ F; note that A(0, ε) = A(fbump , ε) = A(−fbump , ε) For A to succeed for ±fbump , which have diﬀerent integrals, one must have cost(A, F, ε) = O ρ ε the same as the trapezoidal rule algorithm. Why do your instructors tell you that other algorithms are better than the trapezoidal rule? They are only better when you assume that the integrand has more smoothness 7/13
of Our Non-Adaptive Trapezoidal Rule Algorithm Tn(f) = 1 n f(0) 2 + f 1 n + · · · + f n − 1 n + f(1) 2 A(f, ε) = Tn(f) where n = ρ 8ε Solves the integration problem for F = {f : f 1 ρ} Costs O( ρ/ε) function values Is optimal for F = {f : f 1 ρ} T50(f) = 0.3957 What more could we want? 8/13
of Our Non-Adaptive Trapezoidal Rule Algorithm Tn(f) = 1 n f(0) 2 + f 1 n + · · · + f n − 1 n + f(1) 2 A(f, ε) = Tn(f) where n = ρ 8ε Solves the integration problem for F = {f : f 1 ρ} Costs O( ρ/ε) function values Is optimal for F = {f : f 1 ρ} T50(f) = 0.3957 What more could we want? An algorithm that costs O( f 1 /ε) without needing to know f 1 ; cost depends on diﬃculty 8/13
Don’t Work for Balls, But May Work for Cones Adptive algorithms expend eﬀort comensurate with the diﬃculty of the problem If F = , then non-adaptive algorithms are optimal because is Balanced: f ∈ F =⇒ −f ∈ F, and Convex: f, g ∈ F =⇒ θF + (1 − θ)g ∈ F for 0 θ 1 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). 9/13
Don’t Work for Balls, But May Work for Cones Adptive algorithms expend eﬀort comensurate with the diﬃculty of the problem If F = , then non-adaptive algorithms are optimal because is Balanced: f ∈ F =⇒ −f ∈ F, and Convex: f, g ∈ F =⇒ θF + (1 − θ)g ∈ F for 0 θ 1 However, adaptive algorithms may be optimal for non-convex F = (f ∈ F =⇒ cf ∈ F) 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). 9/13
Don’t Work for Balls, But May Work for Cones Adptive algorithms expend eﬀort comensurate with the diﬃculty of the problem If F = , then non-adaptive algorithms are optimal because is Balanced: f ∈ F =⇒ −f ∈ F, and Convex: f, g ∈ F =⇒ θF + (1 − θ)g ∈ F for 0 θ 1 However, adaptive algorithms may be optimal for non-convex F = (f ∈ F =⇒ cf ∈ F) We developed an adaptive trapezoidal rule that optimally solves the integral for F = , a set of functions where the true curvature is not much more than its numerical approximation 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). 9/13
Don’t Work for Balls, But May Work for Cones Adptive algorithms expend eﬀort comensurate with the diﬃculty of the problem If F = , then non-adaptive algorithms are optimal because is Balanced: f ∈ F =⇒ −f ∈ F, and Convex: f, g ∈ F =⇒ θF + (1 − θ)g ∈ F for 0 θ 1 However, adaptive algorithms may be optimal for non-convex F = (f ∈ F =⇒ cf ∈ F) We developed an adaptive trapezoidal rule that optimally solves the integral for F = , a set of functions where the true curvature is not much more than its numerical approximation Common adaptive algorithms in use have no theory 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). 9/13
Don’t Work for Balls, But May Work for Cones We also developed adaptive algorithms for Univariate integration by Simpson’s rule Univariate function approximation and optimization Multivariate integration made available in the Guaranteed Automatic Integration Library (GAIL) Zhang, Y. Adaptive Quadrature. PhD thesis (Illinois Institute of Technology, 2018+). 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). H., F. 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, L. 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, Jagadeeswaran, R. & H., F. J. Automatic Bayesian Cubature. in preparation. 2018+. Choi, S.-C. T., Ding, Y., H., F. J., Jiang, L., Jiménez Rugama, L. A., Li, D., Jagadeeswaran, R., Tong, X., Zhang, K., et al. GAIL: Guaranteed Automatic Integration Library (Versions 1.0–2.2). MATLAB software. 2013–2017. http://gailgithub.github.io/GAIL_Dev/. 10/13
Be Done Find F = for which MATLAB’s integral_g or Chebfun are guaranteed to solve the integration problem Locally adaptive univariate integration Higher order function approximation Adaptive multi-level Monte Carlo methods Adaptive multivariate function approximation Other error criteria, such as relative error or a hybrid of absolute and relative error Combining our software eﬀorts with that of other research groups The MathWorks, Inc. MATLAB 9.5. (Natick, MA, 2018). Hale, N., Trefethen, L. N. & Driscoll, T. A. Chebfun Version 5.7. (2017). Giles, M. Multilevel Monte Carlo Methods. Acta Numer. 24, 259–328 (2015). Ding, Y., H., F. J. & Jiménez Rugama, L. A. An Adaptive Algorithm Employing Continuous Linear Functionals. 2018+. H., F. J., Jiménez Rugama, L. 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.) 597–619 (Springer-Verlag, 2018). doi:10.1007/978-3-319-72456-0. 11/13
for You Join the Society for Industrial and Applied Mathematics (SIAM) for free Click here Join the American Statististcal Association (ASA) for $25 per year Click here 12/13
J., Novak, E. & Rudolf, D. Solvable Integration Problems and Optimal Sample Size Selection. submitted for publication. 2018+. 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). Zhang, Y. Adaptive Quadrature. PhD thesis (Illinois Institute of Technology, 2018+). 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). H., F. 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, L. A. Reliable Adaptive Cubature Using Digital Sequences. in Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014 13/13
R. & Nuyens, D.) 163. arXiv:1410.8615 [math.NA] (Springer-Verlag, Berlin, 2016), 367–383. Jagadeeswaran, R. & H., F. J. Automatic Bayesian Cubature. in preparation. 2018+. 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/. The MathWorks, Inc. MATLAB 9.5. (Natick, MA, 2018). Hale, N., Trefethen, L. N. & Driscoll, T. A. Chebfun Version 5.7. (2017). Giles, M. Multilevel Monte Carlo Methods. Acta Numer. 24, 259–328 (2015). Ding, Y., H., F. J. & Jiménez Rugama, L. A. An Adaptive Algorithm Employing Continuous Linear Functionals. 2018+. 13/13
J., Jiménez Rugama, L. 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.) 597–619 (Springer-Verlag, 2018). doi:10.1007/978-3-319-72456-0. 13/13