Applied Mathematics Center for Interdisciplinary Scientific Computation Illinois Institute of Technology hickernell@iit.edu mypages.iit.edu/~hickernell Thanks to the GAIL team, esp. Lan Jiang, Lluís Antoni Jiménez Rugama, and Jagadees Rathinavel NSF-DMS-1522687 and NSF-DMS-1638521 (SAMSI) Matrix Workshop on the Frontiers of High Dimensional Computation, June 11, 2018
Do We Stop? Compute an integral µ(f) = Rd f(x) (x) dx Bayesian inference, financial risk, statistical physics, ... Desired solution: An adaptive algorithm, ^ µ(·, ·) of the form ^ µ(f, ε) = w0,n + n i=1 wi,n f(xi ), e.g., w0,n = 0, w1,n = · · · = wn,n = 1 n where n, {xi }∞ i=1 , w0,n , and w = (wi,n )n i=1 are chosen to guarantee µ(f) − ^ µ(f, ε) ε with high probability ∀ε > 0, reasonable f 2/13
Do We Stop? µ(f) = Rd f(x) (x) dx, ^ µ(f, ε) = w0,n + n i=1 wi,n f(xi ), Want µ(f) − ^ µ(f, ε) ε with high probability ∀ε > 0, reasonable f Possible solutions Drawbacks IID sampling: n = 2.58^ σ ε 2 Does CLT hold yet? How to determine ^ σ? Low discrepancy sampling: choose n to make disc({xi }n i=1 ) Var(f) ε Don’t know Var(f); computing disc({xi }n i=1 ) may be expensive Randomized low discrepancy sampling: look at spread of sample means from random replications of low discrepancy points How many replications? What measure of spread? 3/13
Do We Stop? µ(f) = Rd f(x) (x) dx, ^ µ(f, ε) = w0,n + n i=1 wi,n f(xi ), Want µ(f) − ^ µ(f, ε) ε with high probability ∀ε > 0, reasonable f Possible solutions Drawbacks IID sampling: n = 2.58^ σ ε 2 Does CLT hold yet? How to determine ^ σ? Low discrepancy sampling: choose n to make disc({xi }n i=1 ) Var(f) ε Don’t know Var(f); computing disc({xi }n i=1 ) may be expensive Randomized low discrepancy sampling: look at spread of sample means from random replications of low discrepancy points How many replications? What measure of spread? Identify a data-driven measure of error Identify a cone of nice integrands for which that data-driven error can be proven to hold Is integrand really in the cone? Well, sometimes there are necessary conditions that can be checked 3/13
IID Sampling Replace the Central Limit Theorem by Berry-Esseen inequalities: µ := X f(x) (x)dx =: E[f(X)], C = {f ∈ L4(X, ) : kurt(f(X)) κup } For nσ determined by κup , and x1 , x2 , . . . IID let σ2 up = (1.2)2 nσ − 1 nσ i=1 f(xi ), ^ µ(f, ε) = 1 nµ nσ+nµ i=nσ+1 f(xi ), where nµ = min n ∈ N : Φ − √ nε/σup + 1 √ n min 0.3328(κ3/4 up + 0.429), 18.1139κ3/4 up 1 + √ nδ/σup 3 0.25% Then P{|µ(f) − ^ µ(f, ε)| ε} 99%. 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. 4/13
Quasi-Monte Carlo Sampling µ = [0,1]d f(x) dx = f(0), f(x) = k f(k)φk(x), f(k) = [0,1]d f(x)φk (x) dx φk are complex exponentials (for lattices) or Walsh functions (for nets) C = integrands whose true Fourier coefficients do not decay erratically The error bound is derived in terms of discrete transform, fn (k) = 1 n n i=1 f(xi )φk (xi ) ∀k, can be computed in O(n log(n)) operations C(n) moderatek fn (k) ε =⇒ |µ(f) − ^ µ(f, ε)| ε H., F. J. & Ll. A. Jiménez Rugama. 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, Ll. A. Jiménez Rugama & 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, H., F. J. et al. Adaptive Quasi-Monte Carlo Methods for Cubature. in Contemporary Computational Mathematics — a celebration of the 80th birthday of Ian Sloan (eds Dick, J. et al.) (Springer-Verlag, 2018), 597–619. doi:10.1007/978-3-319-72456-0. 5/13
QMC Sampling Assuming a Random f Assume f ∼ GP(m, s2Cθ), a sample from a Gaussian process. Defining c = [0,1]d [0,1]d Cθ(t, x) dtdx, c = [0,1]d Cθ(t, xi ) dt n i=1 , C = Cθ(xi , xj ) n i,j=1 and choosing the weights as w0 = m[1 − cTC−11], w = C−1c, ^ µ(f, ε) = w0 + wTf, f = f(xi ) n i=1 . yields an unbiased approximation: µ(f) − ^ µ(f, ε) f = y ∼ N 0, s2(c − cTC−1c) Choosing n large enough to make 2.58s c − cTC−1c ε, assures us that Pf [|µ(f) − ^ µ(f, ε)| ε] 99%. Issues requiring attention: choosing parameters of the covariance kernel, expensive matrix operations 7/13
Make This Approach Practical Use MLE to estimate parameters Choose covariance kernels that match the low discrepancy design c = [0,1]d [0,1]d Cθ(t, x) dtdx = 1, c = [0,1]d Cθ(t, xi ) dt n i=1 = 1 C = Cθ(xi , xj ) n i,j=1 = 1 n VΛVH, Λ = diag(λ), V1 = v1 = 1 VTz is O(n log n), λ = VTC1 , C−11 = 1 λ1 Dealing with the fast transformed data, ^ y = VTy, where y = f(xi ) n i=1 , it follows that ^ µ(f, ε) = ^ y1 n = 1 n n i=1 f(xi ), Pf [|µ(f) − ^ µ(f, ε)| ε] 99%., provided 2.58 1 − n λ1 1 n2 n i=2 |yi |2 λi ε Jagadeeswaran, R. & H., F. J. Automatic Bayesian Cubature. in preparation. 2018+. 8/13
of Matching Covariance Kernels Typically the domain of f is [0, 1)d, and C(t, x) = C(x − t mod 1) integration lattices C(x ⊕ t) Sobol’ sequences, ⊕ means digitwise addition modulo 2 E.g., for integration lattices C(x, t) = d k=1 [1 − θ1 (−1)θ2 B2θ2 (xk − tk mod 1)], θ1 > 0, θ2 ∈ N 9/13
Work Bayesian cubature with higher order nets and smoother kernels, control variates Stopping criteria for multilevel and multivariate decomposition methods Community QMC software that combines the efforts of several research groups Skeleton with clearly defined properties for different kinds of objects Plug-and-play functionality %% Example stopObj = CLTStopping %stopping criterion distribObj = IIDDistribution; %IID sampling with uniform distribution [sol, out] = integrate(KeisterFun, distribObj, stopObj) stopObj.absTol = 1e−3; %decrease tolerance [sol, out] = integrate(KeisterFun, distribObj, stopObj) >> IntegrationExample sol = 0.428567222603452 sol = 0.425203913946775 12/13
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. & Ll. A. Jiménez Rugama. 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. Ll. A. Jiménez Rugama & 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. H., F. J., Jiménez Rugama, Ll. A. & Li, D. Adaptive Quasi-Monte Carlo Methods for Cubature. in Contemporary Computational Mathematics — a celebration of the 80th birthday of Ian Sloan (eds Dick, J., Kuo, F. Y. & Woźniakowski, H.) (Springer-Verlag, 2018), 597–619. doi:10.1007/978-3-319-72456-0. Jagadeeswaran, R. & H., F. J. Automatic Bayesian Cubature. in preparation. 2018+. 13/13
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/. (eds Cools, R. & Nuyens, D.) Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014. 163 (Springer-Verlag, Berlin, 2016). 13/13