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

SIAM UQ 2020-Bayesian Cubature

SIAM UQ 2020-Bayesian Cubature

Fred J. Hickernell

April 19, 2020
Tweet

More Decks by Fred J. Hickernell

Other Decks in Research

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
    [email protected] 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.

    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