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
   

   View Slide

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
   

   View Slide

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
   

   View Slide

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
   

   View Slide

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
   

   View Slide

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
   

   View Slide

7. Introduction Algorithms So Far Example What Comes Next References
   Integrands Inside and Outside the Cone
   6/13
   

   View Slide

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
   

   View Slide

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
   

   View Slide

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
    

    View Slide

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
    

    View Slide

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
    

    View Slide

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
    

    View Slide

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
    

    View Slide

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
    

    View Slide

16. Thank you
    Slides available on SpeakerDeck at
    speakerdeck.com/fjhickernell/matrix-week-two-2018-june
    

    View Slide

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
    

    View Slide

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
    

    View Slide