Upgrade to Pro — share decks privately, control downloads, hide ads and more …

MATRIX Week Two 2018 June

MATRIX Week Two 2018 June

Automatic cubature

Fred J. Hickernell

June 11, 2018
Tweet

More Decks by Fred J. Hickernell

Other Decks in Research

Transcript

  1. Automatic Algorithms for Multidimensional Integration
    Fred J. Hickernell
    Department of Applied Mathematics
    Center for Interdisciplinary Scientific Computation
    Illinois Institute of Technology
    [email protected] 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

    − 1

    i=1
    f(xi
    ), ^
    µ(f, ε) =
    1

    nσ+nµ
    i=nσ+1
    f(xi
    ), where

    = 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