78

# MATRIX Week Two 2018 June

Automatic cubature June 11, 2018

## Transcript

1. Automatic Algorithms for Multidimensional Integration
Fred J. Hickernell
Department of Applied Mathematics
Center for Interdisciplinary Scientiﬁc 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

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, ﬁnancial 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
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
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

− 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 Conﬁdence 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 coeﬃcients 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. Deﬁning
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 diﬀ 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 eﬀorts of several research groups
Skeleton with clearly deﬁned properties for diﬀerent 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 Conﬁdence
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