SIAM UQ 2020-Bayesian Cubature

Transcript

1. The Successes and Challenges of Automatic Bayesian Cubature
   Fred J. Hickernell
   Department of Applied Mathematics
   Center for Interdisciplinary Scientific Computation
   Illinois Institute of Technology
   hickernell@iit.edu mypages.iit.edu/~hickernell
   with Jagadeeswaran R. and the GAIL team
   partially supported by NSF-DMS-1522687 and NSF-DMS-1638521 (SAMSI)
   Thanks to the special session organizers
   SIAM-UQ 2020, March 2020
   

   View Slide

2. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
   The Guaranteed Automatic Integration Library (GAIL) and QMCPy Teams
   Sou-Cheng Choi (Chief Data Scientist, Kamakura)
   Yuhan Ding (IIT PhD ’15, Lecturer, IIT)
   Lan Jiang (IIT PhD ’16, Compass)
   Lluís Antoni Jiménez Rugama (IIT PhD ’17, UBS)
   Jagadeeswaran Rathinavel (IIT PhD ’19, Wi-Tronix)
   Aleksei Sorokin (IIT BS + MAS ’21 exp.)
   Tong Xin (IIT MS, UIC PhD ’20 exp.)
   Kan Zhang (IIT PhD ’20 exp.)
   Yizhi Zhang (IIT PhD ’18, Jamran Int’l)
   Xuan Zhou (IIT PhD ’15, JP Morgan)
   and others
   Adaptive software libraries GAIL and QMCPy
   2/15
   

   View Slide

3. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
   Problem f
   Want fixed tolerance cubature ALG : L1[0, 1]d × (0, ∞) → R such that
   [0,1]d
   f(x) dx − ALG(f, ε) ε ∀ε > 0, for reasonable f
   3/15
   

   View Slide

4. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
   Problem f
   Want fixed tolerance cubature ALG : L1[0, 1]d × (0, ∞) → R such that
   Pf
   [0,1]d
   f(x) dx − ALG(f, ε) ε 99% ∀ε > 0, f ∼ GP(m, s2Cθ)
   3/15
   

   View Slide

5. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
   Problem f
   Want fixed tolerance cubature ALG : L1[0, 1]d × (0, ∞) → R such that
   Pf
   [0,1]d
   f(x) dx − ALG(f, ε) ε 99% ∀ε > 0, f ∼ GP(m, s2Cθ)
   design or node array X = (x1
   , . . . , xn
   )T ∈ [0, 1]n×d, function data f = f(X) ∈ Rn
   c0θ =
   [0,1]d×[0,1]d
   Cθ(t, x) dtdx > 0, cθ =
   [0,1]d
   Cθ(X, t) dt ∈ [0, 1]n, Cθ = Cθ(X, X) ∈ [0, 1]n×n
   [0,1]d
   f(x) dx (f = y) ∼ N m[1 − cT
   θ
   C−1
   θ
   1] + cT
   θ
   C−1y
   ALG(f,ε)
   , s2(c0θ − cT
   θ
   C−1
   θ
   cθ)
   Choosing n large enough to make 2.58s c0θ − cT
   θ
   C−1
   θ
   cθ ε would seem to achieve our goal
   3/15
   

   View Slide

6. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
   Problem f
   Want fixed tolerance cubature ALG : L1[0, 1]d × (0, ∞) → R such that
   Pf
   [0,1]d
   f(x) dx − ALG(f, ε) ε 99% ∀ε > 0, f ∼ GP(m, s2Cθ)
   design or node array X = (x1
   , . . . , xn
   )T ∈ [0, 1]n×d, function data f = f(X) ∈ Rn
   c0θ =
   [0,1]d×[0,1]d
   Cθ(t, x) dtdx > 0, cθ =
   [0,1]d
   Cθ(X, t) dt ∈ [0, 1]n, Cθ = Cθ(X, X) ∈ [0, 1]n×n
   [0,1]d
   f(x) dx (f = y) ∼ N m[1 − cT
   θ
   C−1
   θ
   1] + cT
   θ
   C−1y
   ALG(f,ε)
   , s2(c0θ − cT
   θ
   C−1
   θ
   cθ)
   Choosing n large enough to make 2.58s c0θ − cT
   θ
   C−1
   θ
   cθ ε would seem to achieve our goal
   Issues requiring attention
   Inference of m, s, θ, and Cθ
   , which urn f comes from
   Ill-conditioning and numerical cost of vector-matrix calculations
   Whether Gaussian process is a reasonable assumption
   3/15
   

   View Slide

7. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
   Problem f
   Want fixed tolerance cubature ALG : L1[0, 1]d × (0, ∞) → R such that
   Pf
   [0,1]d
   f(x) dx − ALG(f, ε) ε 99% ∀ε > 0, f ∼ GP(m, s2Cθ)
   design or node array X = (x1
   , . . . , xn
   )T ∈ [0, 1]n×d, function data f = f(X) ∈ Rn
   c0θ =
   [0,1]d×[0,1]d
   Cθ(t, x) dtdx > 0, cθ =
   [0,1]d
   Cθ(X, t) dt ∈ [0, 1]n, Cθ = Cθ(X, X) ∈ [0, 1]n×n
   [0,1]d
   f(x) dx (f = y) ∼ N m[1 − cT
   θ
   C−1
   θ
   1] + cT
   θ
   C−1y
   ALG(f,ε)
   , s2(c0θ − cT
   θ
   C−1
   θ
   cθ)
   Choosing n large enough to make 2.58s c0θ − cT
   θ
   C−1
   θ
   cθ ε would seem to achieve our goal
   Issues requiring attention
   Inference of m, s, θ, and Cθ
   , which urn f comes from
   Ill-conditioning and numerical cost of vector-matrix calculations
   Whether Gaussian process is a reasonable assumption
   3/15
   

   View Slide

8. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
   Inferring Gaussian Process Parameters for GP(m, s2Cθ
   ) f
   Using empirical Bayes
   mEB
   =
   1TC−1
   θ
   y
   1TC−1
   θ
   1
   , s2
   EB
   =
   1
   n
   yT C−1
   θ
   −
   C−1
   θ
   11TC−1
   θ
   1TC−1
   θ
   1
   y,
   θEB
   = argmin
   θ
   log yT C−1
   θ
   −
   C−1
   θ
   11TC−1
   θ
   1TC−1
   θ
   1
   y +
   1
   n
   log(det(Cθ)) ,
   ALG(f, ε) =
   (1 − 1TC−1
   θ
   cθ)1
   1TC−1
   θ
   1
   + cθ
   T
   C−1
   θ
   y when 2.58sEB
   c0θ − cT
   θ
   C−1
   θ
   cθ ε
   Jagadeeswaran, R. & H., F. J. Fast Automatic Bayesian Cubature Using Lattice Sampling. Stat. Comput. 29, 1215–1229
   (2019). 4/15
   

   View Slide

9. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
   Inferring Gaussian Process Parameters for GP(m, s2Cθ
   ) f
   Using empirical Bayes
   mEB
   =
   1TC−1
   θ
   y
   1TC−1
   θ
   1
   , s2
   EB
   =
   1
   n
   yT C−1
   θ
   −
   C−1
   θ
   11TC−1
   θ
   1TC−1
   θ
   1
   y,
   θEB
   = argmin
   θ
   log yT C−1
   θ
   −
   C−1
   θ
   11TC−1
   θ
   1TC−1
   θ
   1
   y +
   1
   n
   log(det(Cθ)) ,
   ALG(f, ε) =
   (1 − 1TC−1
   θ
   cθ)1
   1TC−1
   θ
   1
   + cθ
   T
   C−1
   θ
   y when 2.58sEB
   c0θ − cT
   θ
   C−1
   θ
   cθ ε
   Ill-conditioning and numerical cost of vector-matrix calculations
   Whether Gaussian process is a reasonable assumption
   Jagadeeswaran, R. & H., F. J. Fast Automatic Bayesian Cubature Using Lattice Sampling. Stat. Comput. 29, 1215–1229
   (2019). 4/15
   

   View Slide

10. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Fast Bayesian Transforms in General
    Find a kernel Cθ
    to match the design X so that
    Cθ =
    1
    n
    VΛθ
    VH, VH = nV−1, V = (v1
    , . . . , vn
    )T = (V1
    , . . . , Vn
    ) known analytically
    v1
    = V1
    = 1, cθ = 1, b := VHb requires only O(n log(n)) operations ∀b.
    Cθ
    is a fast Bayesian transform kernel and b → VHb a fast Bayesian transform (FBT)
    Then by empirical Bayes
    y = FBT of function data y, λθ = diag(Λθ) = (λθ,1
    , . . . , λθ,n
    )T = Cθ,1
    = FBT of first column of Cθ
    θEB
    = argmin
    θ
    log
    n
    i=2
    |yi
    |2
    λθ,i
    +
    1
    n
    n
    i=1
    log(λθ,i
    )
    ALG(f, ε) =
    y1
    n
    =
    1
    n
    n
    i=1
    yi
    = sample mean when
    2.58
    n
    n
    i=2
    |yi
    |2
    λθ,i
    1 −
    n
    λθ,1
    ε
    Cost is O(n log(n)) times the number of iterations for optimizing θ
    5/15
    

    View Slide

11. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Types of FBT Kernels
    Cθ(t, x) = Kθ(x t), {xi
    }2m
    i=1
    = affine shift of a group under ⊕ for m = 0, 1, . . .
    6/15
    

    View Slide

12. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Types of FBT Kernels
    Cθ(t, x) = Kθ(x t), {xi
    }2m
    i=1
    = affine shift of a group under ⊕ for m = 0, 1, . . .
    Shifted Lattice Nodes, ⊕ = addition mod1 Scrambled Sobol’ Nodes, ⊕ = bitwise addition
    FBT = Fast Fourier Transform FBT = Fast Walsh Transform
    6/15
    

    View Slide

13. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Flexible FBT Kernel for Lattice Node Designs
    Cθ(t, x) = Kθ(x t) must be positive definite, where ⊕ = addition mod1. Common example is
    Kθ(x) =
    d
    j=1
    [1 + a B2r
    (xj
    )] =
    d
    j=1
    [1 + aκr
    (xj
    )], θ = (a, r) ∈ (0, ∞) × (1, ∞), r = 2r
    κr
    (x) :=
    |k| 1
    exp(2π
    √
    −1kx)
    |k|r
    B2r
    closed form, but r ∈ N; κr
    defined for r > 1, but infinite sum
    7/15
    

    View Slide

14. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Flexible FBT Kernel for Lattice Node Designs
    Cθ(t, x) = Kθ(x t) must be positive definite, where ⊕ = addition mod1. Common example is
    Kθ(x) =
    d
    j=1
    [1 + a B2r
    (xj
    )] =
    d
    j=1
    [1 + aκr
    (xj
    )], θ = (a, r) ∈ (0, ∞) × (1, ∞), r = 2r
    κr
    (x) :=
    |k| 1
    exp(2π
    √
    −1kx)
    |k|r
    B2r
    closed form, but r ∈ N; κr
    defined for r > 1, but infinite sum
    But all we need to compute the integral and credible interval is y and
    λθ = Cθ,1
    , where Cθ,1
    = Kθ(xi
    x1
    )
    n
    i=1
    , which only depends on κr
    ( /n) for = 0, . . . , n − 1
    7/15
    

    View Slide

15. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Flexible FBT Kernel for Lattice Node Designs
    Cθ(t, x) = Kθ(x t) must be positive definite, where ⊕ = addition mod1. Common example is
    Kθ(x) =
    d
    j=1
    [1 + a B2r
    (xj
    )] =
    d
    j=1
    [1 + aκr
    (xj
    )], θ = (a, r) ∈ (0, ∞) × (1, ∞), r = 2r
    κr
    (x) :=
    |k| 1
    exp(2π
    √
    −1kx)
    |k|r
    B2r
    closed form, but r ∈ N; κr
    defined for r > 1, but infinite sum
    But all we need to compute the integral and credible interval is y and
    λθ = Cθ,1
    , where Cθ,1
    = Kθ(xi
    x1
    )
    n
    i=1
    , which only depends on κr
    ( /n) for = 0, . . . , n − 1
    Moreover, κr
    ( /n) n−1
    =0
    can be computed using one FFT details
    7/15
    

    View Slide

16. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Keister’s Example
    Rd
    cos( t ) exp(− t 2) dt =
    [0,1]d
    f(x) dx
    ε = 5 × 10−3, d = 3
    Method MC Lattice Sobol BayesLat BayesSobol
    Absolute Error 1.40 × 10−3 5.20 × 10−4 5.40 × 10−4 4.10 × 10−4 6.80 × 10−4
    Tolerance Met 100% 100% 100% 100% 100%
    n 2 600 000 4100 3900 510 1900
    Time (seconds) 0.1400 0.0110 0.0091 0.0034 0.0410
    Algorithms are implemented in GAIL and soon QMCPy
    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/.
    Choi, S.-C. T., H., F. J., McCourt, M. & Sorokin, A. QMCPy: A quasi-Monte Carlo Python Library. 2020+.
    https://github.com/QMCSoftware/QMCSoftware.
    Keister, B. D. Multidimensional Quadrature Algorithms. Computers in Physics 10, 119–122 (1996). 8/15
    

    View Slide

17. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Gaussian Probability
    (a,b)
    exp −1
    2
    tTΣ−1t
    (2π)d det(Σ)
    dt =
    [0,1]d−1
    f(x) dx by Genz’s transformation
    ε = 1 × 10−4, d = 5, Σ = 0.4I + 0.611T, a = (−∞, . . . , −∞), b ∼
    √
    dU[0, 1]d
    Method MC Lattice Sobol’ BayesLat BayesSobol
    Absolute Error 2.00 × 10−5 5.00 × 10−6 4.10 × 10−6 9.20 × 10−6 3.40 × 10−6
    Tolerance Met 100% 100% 100% 100% 100%
    n 62 000 000 4100 4100 2000 4100
    Time (seconds) 17.0000 0.0110 0.0097 0.0880 0.0950
    Algorithms are implemented in GAIL and soon QMCPy
    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/.
    Choi, S.-C. T., H., F. J., McCourt, M. & Sorokin, A. QMCPy: A quasi-Monte Carlo Python Library. 2020+.
    https://github.com/QMCSoftware/QMCSoftware.
    Genz, A. Comparison of Methods for the Computation of Multivariate Normal Probabilities. Computing Science and Statistics
    25, 400–405 (1993). 9/15
    

    View Slide

18. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Is f a Typical Instance of a Gaussian Process?
    If f ∈ GP(m, s2Cθ), then Z := s−1 n−1VΛ−1/2
    θ
    VH
    C−1/2
    θ
    (f − m1) ∼ N(0, I). Generate the data
    z =
    1
    nsEB
    VΛ−1/2
    θEB
    VH(y − mEB
    1) =
    n
    i=2
    |yi
    |2
    λθEB
    ,i
    −1
    VΛ−1/2
    θEB
    (y − y1
    e1
    )
    Q-Q plots of z vs. standard Gaussian
    10/15
    

    View Slide

19. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Is f a Typical Instance of a Gaussian Process?
    If f ∈ GP(m, s2Cθ), then Z := s−1 n−1VΛ−1/2
    θ
    VH
    C−1/2
    θ
    (f − m1) ∼ N(0, I). Generate the data
    z =
    1
    nsEB
    VΛ−1/2
    θEB
    VH(y − mEB
    1) =
    n
    i=2
    |yi
    |2
    λθEB
    ,i
    −1
    VΛ−1/2
    θEB
    (y − y1
    e1
    )
    Q-Q plots of z vs. standard Gaussian
    Random function fitting kernel with given r, θ
    f(x) = f0
    +
    k∈{1,...,N}d
    fc
    (k) cos(2πkTx) + fs
    (k) sin(2πkTx)
    f0
    , fc
    (k), fs
    (k) IID
    ∼ N 0, a k 0 bd− k 0
    kj=0
    k−r
    j
    θ = a/b, N = 256
    10/15
    

    View Slide

20. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Is f a Typical Instance of a Gaussian Process?
    If f ∈ GP(m, s2Cθ), then Z := s−1 n−1VΛ−1/2
    θ
    VH
    C−1/2
    θ
    (f − m1) ∼ N(0, I). Generate the data
    z =
    1
    nsEB
    VΛ−1/2
    θEB
    VH(y − mEB
    1) =
    n
    i=2
    |yi
    |2
    λθEB
    ,i
    −1
    VΛ−1/2
    θEB
    (y − y1
    e1
    )
    Q-Q plots of z vs. standard Gaussian
    Random function fitting kernel with given r, θ
    f(x) = f0
    +
    k∈{1,...,N}d
    fc
    (k) cos(2πkTx) + fs
    (k) sin(2πkTx)
    f0
    , fc
    (k), fs
    (k) IID
    ∼ N 0, a k 0 bd− k 0
    kj=0
    k−r
    j
    θ = a/b, N = 256
    10/15
    

    View Slide

21. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Is f a Typical Instance of a Gaussian Process?
    If f ∈ GP(m, s2Cθ), then Z := s−1 n−1VΛ−1/2
    θ
    VH
    C−1/2
    θ
    (f − m1) ∼ N(0, I). Generate the data
    z =
    1
    nsEB
    VΛ−1/2
    θEB
    VH(y − mEB
    1) =
    n
    i=2
    |yi
    |2
    λθEB
    ,i
    −1
    VΛ−1/2
    θEB
    (y − y1
    e1
    )
    Q-Q plots of z vs. standard Gaussian
    Random function fitting kernel with given r, θ
    f(x) = f0
    +
    k∈{1,...,N}d
    fc
    (k) cos(2πkTx) + fs
    (k) sin(2πkTx)
    f0
    , fc
    (k), fs
    (k) IID
    ∼ N 0, a k 0 bd− k 0
    kj=0
    k−r
    j
    θ = a/b, N = 256
    10/15
    

    View Slide

22. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Is f a Typical Instance of a Gaussian Process?
    If f ∈ GP(m, s2Cθ), then Z := s−1 n−1VΛ−1/2
    θ
    VH
    C−1/2
    θ
    (f − m1) ∼ N(0, I). Generate the data
    z =
    1
    nsEB
    VΛ−1/2
    θEB
    VH(y − mEB
    1) =
    n
    i=2
    |yi
    |2
    λθEB
    ,i
    −1
    VΛ−1/2
    θEB
    (y − y1
    e1
    )
    Q-Q plots of z vs. standard Gaussian
    Random function fitting kernel with given r, θ
    f(x) = f0
    +
    k∈{1,...,N}d
    fc
    (k) cos(2πkTx) + fs
    (k) sin(2πkTx)
    f0
    , fc
    (k), fs
    (k) IID
    ∼ N 0, a k 0 bd− k 0
    kj=0
    k−r
    j
    θ = a/b, N = 256
    10/15
    

    View Slide

23. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Is f a Typical Instance of a Gaussian Process?
    If f ∈ GP(m, s2Cθ), then Z := s−1 n−1VΛ−1/2
    θ
    VH
    C−1/2
    θ
    (f − m1) ∼ N(0, I). Generate the data
    z =
    1
    nsEB
    VΛ−1/2
    θEB
    VH(y − mEB
    1) =
    n
    i=2
    |yi
    |2
    λθEB
    ,i
    −1
    VΛ−1/2
    θEB
    (y − y1
    e1
    )
    Q-Q plots of z vs. standard Gaussian
    Random function fitting kernel with given r, θ
    f(x) = f0
    +
    k∈{1,...,N}d
    fc
    (k) cos(2πkTx) + fs
    (k) sin(2πkTx)
    f0
    , fc
    (k), fs
    (k) IID
    ∼ N 0, a k 0 bd− k 0
    kj=0
    k−r
    j
    θ = a/b, N = 256
    10/15
    

    View Slide

24. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Is f a Typical Instance of a Gaussian Process?
    If f ∈ GP(m, s2Cθ), then Z := s−1 n−1VΛ−1/2
    θ
    VH
    C−1/2
    θ
    (f − m1) ∼ N(0, I). Generate the data
    z =
    1
    nsEB
    VΛ−1/2
    θEB
    VH(y − mEB
    1) =
    n
    i=2
    |yi
    |2
    λθEB
    ,i
    −1
    VΛ−1/2
    θEB
    (y − y1
    e1
    )
    Q-Q plots of z vs. standard Gaussian
    Random function fitting kernel with given r, θ
    f(x) = f0
    +
    k∈{1,...,N}d
    fc
    (k) cos(2πkTx) + fs
    (k) sin(2πkTx)
    f0
    , fc
    (k), fs
    (k) IID
    ∼ N 0, a k 0 bd− k 0
    kj=0
    k−r
    j
    θ = a/b, N = 256
    10/15
    

    View Slide

25. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Is f a Typical Instance of a Gaussian Process?
    If f ∈ GP(m, s2Cθ), then Z := s−1 n−1VΛ−1/2
    θ
    VH
    C−1/2
    θ
    (f − m1) ∼ N(0, I). Generate the data
    z =
    1
    nsEB
    VΛ−1/2
    θEB
    VH(y − mEB
    1) =
    n
    i=2
    |yi
    |2
    λθEB
    ,i
    −1
    VΛ−1/2
    θEB
    (y − y1
    e1
    )
    Q-Q plots of z vs. standard Gaussian
    Random function fitting kernel with given r, θ
    f(x) = f0
    +
    k∈{1,...,N}d
    fc
    (k) cos(2πkTx) + fs
    (k) sin(2πkTx)
    f0
    , fc
    (k), fs
    (k) IID
    ∼ N 0, a k 0 bd− k 0
    kj=0
    k−r
    j
    θ = a/b, N = 256
    10/15
    

    View Slide

26. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Is f a Typical Instance of a Gaussian Process?
    If f ∈ GP(m, s2Cθ), then Z := s−1 n−1VΛ−1/2
    θ
    VH
    C−1/2
    θ
    (f − m1) ∼ N(0, I). Generate the data
    z =
    1
    nsEB
    VΛ−1/2
    θEB
    VH(y − mEB
    1) =
    n
    i=2
    |yi
    |2
    λθEB
    ,i
    −1
    VΛ−1/2
    θEB
    (y − y1
    e1
    )
    Q-Q plots of z vs. standard Gaussian
    f from Keister’s example
    True r and θ unknown
    10/15
    

    View Slide

27. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Is f a Typical Instance of a Gaussian Process?
    If f ∈ GP(m, s2Cθ), then Z := s−1 n−1VΛ−1/2
    θ
    VH
    C−1/2
    θ
    (f − m1) ∼ N(0, I). Generate the data
    z =
    1
    nsEB
    VΛ−1/2
    θEB
    VH(y − mEB
    1) =
    n
    i=2
    |yi
    |2
    λθEB
    ,i
    −1
    VΛ−1/2
    θEB
    (y − y1
    e1
    )
    Q-Q plots of z vs. standard Gaussian
    f from Keister’s example
    True r and θ unknown
    10/15
    

    View Slide

28. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Is f a Typical Instance of a Gaussian Process?
    If f ∈ GP(m, s2Cθ), then Z := s−1 n−1VΛ−1/2
    θ
    VH
    C−1/2
    θ
    (f − m1) ∼ N(0, I). Generate the data
    z =
    1
    nsEB
    VΛ−1/2
    θEB
    VH(y − mEB
    1) =
    n
    i=2
    |yi
    |2
    λθEB
    ,i
    −1
    VΛ−1/2
    θEB
    (y − y1
    e1
    )
    Q-Q plots of z vs. standard Gaussian
    f from Keister’s example
    True r and θ unknown
    10/15
    

    View Slide

29. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Is f a Typical Instance of a Gaussian Process?
    If f ∈ GP(m, s2Cθ), then Z := s−1 n−1VΛ−1/2
    θ
    VH
    C−1/2
    θ
    (f − m1) ∼ N(0, I). Generate the data
    z =
    1
    nsEB
    VΛ−1/2
    θEB
    VH(y − mEB
    1) =
    n
    i=2
    |yi
    |2
    λθEB
    ,i
    −1
    VΛ−1/2
    θEB
    (y − y1
    e1
    )
    Q-Q plots of z vs. standard Gaussian
    f from Keister’s example
    True r and θ unknown
    10/15
    

    View Slide

30. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Is f a Typical Instance of a Gaussian Process?
    If f ∈ GP(m, s2Cθ), then Z := s−1 n−1VΛ−1/2
    θ
    VH
    C−1/2
    θ
    (f − m1) ∼ N(0, I). Generate the data
    z =
    1
    nsEB
    VΛ−1/2
    θEB
    VH(y − mEB
    1) =
    n
    i=2
    |yi
    |2
    λθEB
    ,i
    −1
    VΛ−1/2
    θEB
    (y − y1
    e1
    )
    Q-Q plots of z vs. standard Gaussian
    f from Keister’s example
    True r and θ unknown
    10/15
    

    View Slide

31. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Is f a Typical Instance of a Gaussian Process?
    If f ∈ GP(m, s2Cθ), then Z := s−1 n−1VΛ−1/2
    θ
    VH
    C−1/2
    θ
    (f − m1) ∼ N(0, I). Generate the data
    z =
    1
    nsEB
    VΛ−1/2
    θEB
    VH(y − mEB
    1) =
    n
    i=2
    |yi
    |2
    λθEB
    ,i
    −1
    VΛ−1/2
    θEB
    (y − y1
    e1
    )
    Q-Q plots of z vs. standard Gaussian
    f from Keister’s example
    True r and θ unknown
    10/15
    

    View Slide

32. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    The Successes of Bayesian Cubature f
    Credible intervals lead to adaptive cubature algorithms that provide answers to the
    desired tolerances
    Low discrepancy sampling provides cubatures requiring fewer samples that IID sampling
    Empirical Bayes and other methods can infer reasonable Gaussian process parameters
    for a family of covairance kernels—choosing a reasonable urn for f
    Covariance kernels that match the low discrepancy sampling facilitate fast Bayeian
    transforms for fitting the Gaussian process parameters, computing the cubature, and
    constructing the credible intervals
    Extra effort provides a richer family of covariance kernels
    11/15
    

    View Slide

33. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    The Challenges of Bayesian Cubature f
    Inferring the Gaussian process parameters and allowing a larger family of kernels takes
    extra computational time
    Overfitting as well as underfitting the kernel are possible
    The Gaussian assumption is not always justified by the data
    Which periodizing transformations of f are appropriate for lattice sampling and its
    matching kernels?
    12/15
    

    View Slide

34. Thank you
    These slides are available at
    speakerdeck.com/fjhickernell/siam-uq-2020-bayesian-cubature
    

    View Slide

35. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    References
    Jagadeeswaran, R. & H., F. J. Fast Automatic Bayesian Cubature Using Lattice Sampling. Stat.
    Comput. 29, 1215–1229 (2019).
    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/.
    Choi, S.-C. T., H., F. J., McCourt, M. & Sorokin, A. QMCPy: A quasi-Monte Carlo Python Library.
    2020+. https://github.com/QMCSoftware/QMCSoftware.
    Keister, B. D. Multidimensional Quadrature Algorithms. Computers in Physics 10, 119–122 (1996).
    Genz, A. Comparison of Methods for the Computation of Multivariate Normal Probabilities.
    Computing Science and Statistics 25, 400–405 (1993).
    14/15
    

    View Slide

36. Introduction Fast Bayesian Transforms Numerical Examples Gaussian Process Diagnostics Summary References
    Computing κ(x) for x = 0, 1/n . . . , 1 − 1/n
    ζ(s, a) :=
    ∞
    m=0
    (a + m)−s Hurwitz zeta function
    κr
    ( /n) =
    |k| 1
    exp(2π
    √
    −1k /n)
    |k|r
    =
    2
    nr
    ∞
    m=1
    1
    |m|r
    +
    n−1
    k=1
    −∞
    m=−∞
    exp(2π
    √
    −1k /n)
    |k + mn|r
    =
    1
    nr
    2ζ(r) +
    n−1
    k=1
    exp(2π
    √
    −1k /n)[ζ(r, k/n) + ζ(r, 1 − k/n)]
    κr
    = (κ( /n))n−1
    =0
    = WHκr
    κr
    := n−r(2ζ(r), ζ(r, 1/n) + ζ(r, (n − 1)/n), . . . , ζ(r, (n − 1)/n) + ζ(r, 1/n))T
    W = (exp(2π
    √
    −1k /n))n−1
    k, =0
    return
    15/15
    

    View Slide