MATRIX Week Two 2018 June

Automatic cubature

Transcript

1. Automatic Algorithms for Multidimensional Integration Fred J. Hickernell Department of 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
2. Introduction Algorithms So Far Example What Comes Next References When 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
3. Introduction Algorithms So Far Example What Comes Next References When 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
4. Introduction Algorithms So Far Example What Comes Next References When 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
5. Introduction Algorithms So Far Example What Comes Next References Automatic 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
6. Introduction Algorithms So Far Example What Comes Next References Automatic 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
7. Introduction Algorithms So Far Example What Comes Next References Integrands Inside and Outside the Cone 6/13
8. Introduction Algorithms So Far Example What Comes Next References Automatic 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
9. Introduction Algorithms So Far Example What Comes Next References To 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
10. Introduction Algorithms So Far Example What Comes Next References Form 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
11. Introduction Algorithms So Far Example What Comes Next References Gaussian Probability µ = [a,b] exp −1 2 tTΣ−1t (2π)d det(Σ) dt = [0,1]d−1 f(x) dx For some typical choice of a, b, Σ, d = 3; µ ≈ 0.6763 Rel. Error Median Worst 10% Worst 10% Tolerance Method Error Accuracy n Time (s) 1E−2 IID 5E−4 100% 8.1E4 0.020 1E−2 Sh. Latt. 1E−5 100% 1.0E3 0.004 1E−2 Scr. Sobol’ 4E−5 100% 1.0E3 0.005 1E−2 Bayes. Latt. 1E−5 100% 1.0E3 0.002 These algorithms are implemented in GAIL 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/. 10/13
12. Introduction Algorithms So Far Example What Comes Next References Gaussian Probability µ = [a,b] exp −1 2 tTΣ−1t (2π)d det(Σ) dt = [0,1]d−1 f(x) dx For some typical choice of a, b, Σ, d = 3; µ ≈ 0.6763 Rel. Error Median Worst 10% Worst 10% Tolerance Method Error Accuracy n Time (s) 1E−3 IID 9E−5 100% 2.0E6 0.400 1E−3 Sh. Latt. 8E−6 100% 2.0E3 0.007 1E−3 Scr. Sobol’ 2E−5 100% 2.0E3 0.006 1E−3 Bayes. Latt. 1E−5 100% 1.0E3 0.002 These algorithms are implemented in GAIL 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/. 10/13
13. Introduction Algorithms So Far Example What Comes Next References Gaussian Probability µ = [a,b] exp −1 2 tTΣ−1t (2π)d det(Σ) dt = [0,1]d−1 f(x) dx For some typical choice of a, b, Σ, d = 3; µ ≈ 0.6763 Rel. Error Median Worst 10% Worst 10% Tolerance Method Error Accuracy n Time (s) 1E−4 Sh. Latt. 4E−7 100% 8.2E3 0.014 1E−4 Scr. Sobol’ 4E−7 100% 1.6E4 0.018 1E−4 Bayes. Latt. 5E−7 100% 8.2E3 0.010 These algorithms are implemented in GAIL 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/. 10/13
14. Introduction Algorithms So Far Example What Comes Next 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 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 , T = 1/4, d = 13 here Abs. Error Median Worst 10% Worst 10% Tolerance Method 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 11/13
15. Introduction Algorithms So Far Example What Comes Next References Future 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
16. Thank you Slides available on SpeakerDeck at speakerdeck.com/fjhickernell/matrix-week-two-2018-june
17. Introduction Algorithms So Far Example What Comes Next 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. & 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
18. Introduction Algorithms So Far Example What Comes Next References 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/. (eds Cools, R. & Nuyens, D.) Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014. 163 (Springer-Verlag, Berlin, 2016). 13/13