Understand what quasi-Monte Carlo (qMC) methods are • Try qMC in place of simple (or IID) Monte Carlo • Use qMC properly • Feel free to interrupt and ask questions
Understand what quasi-Monte Carlo (qMC) methods are • Try qMC in place of simple (or IID) Monte Carlo • Use qMC properly • Feel free to interrupt and ask questions • Join our friendly qMC research community
Understand what quasi-Monte Carlo (qMC) methods are • Try qMC in place of simple (or IID) Monte Carlo • Use qMC properly • Feel free to interrupt and ask questions • Join our friendly qMC research community Try the computations at https://tinyurl.com/QMCTutorialNotebook
x0 , x1 , … • Discrepancy (quality) measures for x0 , x1 , … • Choosing so that n |μ − ̂ μn | ≤ ε • Making our original problem look like the above • Ongoing research μ := expectation 𝔼 [f( X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn
[f( X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn f(X) = option payo ff underground water pressure with random rock porosity pixel intensity from random ray option price average water pressure average pixel intensity = μ
[f( X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn f(X) = option payo ff underground water pressure with random rock porosity pixel intensity from random ray option price average water pressure average pixel intensity = μ may be dozens or hundreds d
• Discrepancy (quality) measures for • Choosing so that • Making our original problem look like the above • Ongoing research x0 , x1 , … x0 , x1 , … n |μ − ̂ μn | ≤ ε μ := expectation 𝔼 [f( X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn
:= expectation 𝔼 [f( X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn Grids do not fi ll space well Hard to extend
i0 z0 ⊕ i1 z1 ⊕ i2 z2 ⊕ ⋯ ⊕ Δ, i = 0,1,…,2m − 1, m ∈ ℕ0 ⊕ = bitwise addition, zm ∈ [0,1)d, Δ ∼ 𝒰 [0,1]d, zm random Every tile has the same # of points How small a tile can be depends on d 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1
i0 z0 ⊕ i1 z1 ⊕ i2 z2 ⊕ ⋯ ⊕ Δ, i = 0,1,…,2m − 1, m ∈ ℕ0 ⊕ = bitwise addition, zm ∈ [0,1)d, Δ ∼ 𝒰 [0,1]d, zm random Every tile has the same # of points How small a tile can be depends on d 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 Random scrambles & shifts make unbiased move the fi rst point away from the origin ̂ μn
drop or skip points • Use randomized LD sequences to avoid points on the boundaries and to make your answers unbiased; you may also gain in convergence rate n = 2m
• Discrepancy (quality) measures for • Choosing so that • Making our original problem look like the above • Ongoing research x0 , x1 , … x0 , x1 , … n |μ − ̂ μn | ≤ ε μ := expectation 𝔼 [f( X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn
decay practically for small • decay practically for large • Discrepancy increases with 𝒪 (n−1+δ) ∀d 𝒪 (n−1+δ) d 𝒪 (n−1/2) d d |μ − ̂ μn | ≤ discrepancy({xi }n−1 i=0 ) variation(f )
∫ [0,1]d f(x) dx = ∫ [0,1]2d K(t, x) dt dx = ( 13 12) d/2 • Dividing discrepancy by helps a bit • LD still (usually) beats IID, although not in convergence order discrepancy(∅)
13 12 ) d must be > ( 13 12 ) d − 2α d ∏ ℓ=1 [1 + 1 2 ( x0ℓ − 1/2 − x0ℓ − 1/2 2 )] +α2 d ∏ ℓ=1 [1 + x0ℓ − 1/2 ] = sup ∥f∥ℋ ≤1 |μ−α ̂ μ1 |2 for optimal α Is this the right and corresponding discrepancy? K
tractable if the work required to solve them to error grows slower than exponentially in • Tractability requires shrinking • Hopefully, the encountered in practice have moderate even with decaying coordinate weights ≤ ε ε−1 {f : variation(f) ≤ 1} f variation(f) |μ − ̂ μn | ≤ discrepancy({xi }n−1 i=0 ) variation(f )
but … discrepancy({xn−1 i=0 }) 𝒪 (dn2) • Computing for shifted lattices and nets requires only operations to compute, because … 𝔼 [(discrepancy2({xshift i }n−1 i=0 ))] 𝒪 (dn) 𝔼 [(discrepancy2({xshift i }n−1 i=0 ; K))] = discrepancy2({xi }n−1 i=0 ; ˜ K ) = 1 n n−1 ∑ i=0 ˜ K (xi , x0 ) − ∫ [0,1]d ˜ K (x, x0 ) dx where is a (digital) shift invariant version of , e.g., for lattices ˜ K K ˜ K (t, x) = ∫ [0,1]d K(t + Δ mod 1, x + Δ mod 1) dΔ
may allow faster decay of the error—see lattices for periodic with higher order smoothness and higher order nets for with higher order smoothness f K f f Is worst case error analysis too pessimistic?
may allow faster decay of the error—see lattices for periodic with higher order smoothness and higher order nets for with higher order smoothness f K f f • For randomized , {xn−1 i=0 } Swapping the order may lead to an extra 𝒪 (n−1/2) Is worst case error analysis too pessimistic? 𝔼 {xi }n−1 i=0 sup variation(f)≤1 |μ − ̂ μn |2 = 𝔼 {xi }n−1 i=0 discrepancy2({xi }n−1 i=0 ) ≥ sup variation(f)≤1 𝔼 {xi }n−1 i=0 |μ − ̂ μn |2
may allow faster decay of the error—see lattices for periodic with higher order smoothness and higher order nets for with higher order smoothness f K f f • For randomized , {xn−1 i=0 } Swapping the order may lead to an extra 𝒪 (n−1/2) • [H18] surveys worst case, randomized, and Bayesian error analyses Is worst case error analysis too pessimistic? 𝔼 {xi }n−1 i=0 sup variation(f)≤1 |μ − ̂ μn |2 = 𝔼 {xi }n−1 i=0 discrepancy2({xi }n−1 i=0 ) ≥ sup variation(f)≤1 𝔼 {xi }n−1 i=0 |μ − ̂ μn |2
• Discrepancy (quality) measures for • Choosing so that • Making our original problem look like the above • Ongoing research x0 , x1 , … x0 , x1 , … n |μ − ̂ μn | ≤ ε μ := expectation 𝔼 [f( X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn
:= expectation 𝔼 [f( X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn • is for low discrepancy points versus for IID discrepancy({xi }n−1 i=0 ) 𝒪 (n−1) 𝒪 (n−1/2) • There is an explicit formula for discrepancy, but the variation is hard to compute in practice |μ − ̂ μn | ≤ discrepancy({xi }n−1 i=0 ) variation(f ) want |μ − ̂ μn | ≤ ε
:= expectation 𝔼 [f( X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn • is for low discrepancy points versus for IID discrepancy({xi }n−1 i=0 ) 𝒪 (n−1) 𝒪 (n−1/2) • There is an explicit formula for discrepancy, but the variation is hard to compute in practice • One method is random replications plus the Student’s con fi dence interval [LENOT24] t |μ − ̂ μn | ≤ discrepancy({xi }n−1 i=0 ) variation(f ) want |μ − ̂ μn | ≤ ε
X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn , want |μ − ̂ μn | ≤ ε • For lattices, , , ̂ f(k) := ∫ [0,1]d f(x) 𝖾 −2π −1k′  x dx f(x) = ∑ k∈ℤd ̂ f(k) 𝖾 2π −1k′  x μ = ̂ f(0) • due to aliasing ( constant for all ) μ − ̂ μn 𝖾 2π −1k′  xi i • Approximate by a one-dimensional FFT of ̂ f(k) {f(xi )}n−1 i=0 • If the decay in a reasonable manner, their FFT approximations can be used to provide a rigorous data-driven bound on ̂ f(k) μ − ̂ μn
X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn , want |μ − ̂ μn | ≤ ε • For lattices, , , ̂ f(k) := ∫ [0,1]d f(x) 𝖾 −2π −1k′  x dx f(x) = ∑ k∈ℤd ̂ f(k) 𝖾 2π −1k′  x μ = ̂ f(0) • due to aliasing ( constant for all ) μ − ̂ μn 𝖾 2π −1k′  xi i • Approximate by a one-dimensional FFT of ̂ f(k) {f(xi )}n−1 i=0 • If the decay in a reasonable manner, their FFT approximations can be used to provide a rigorous data-driven bound on ̂ f(k) μ − ̂ μn • Similarly for digital nets
X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn , want |μ − ̂ μn | ≤ ε • If is assumed to be a Gaussian stochastic process with covariance kernel , where the hyper-parameters are properly tuned, then one may construct a Bayesian credible interval for the error • If is chosen to be shift invariant and lattice/digital sequences are used, then the computation normally required is reduced to f K K 𝒪 (n3) 𝒪 (n log n)
• Discrepancy (quality) measures for • Choosing so that • Making our original problem look like the above • Ongoing research x0 , x1 , … x0 , x1 , … n |μ − ̂ μn | ≤ ε μ := expectation 𝔼 [f( X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn
something wonderful ⏞ ⋯ = expectation 𝔼 [f( X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn If for , then • is not unique • Want to be as small as possible (an art) • Often choose , but not necessary Z = Ψ(X) X ∼ 𝒰 [0,1]d f(x) = g(Ψ(x)) λ(Ψ(x)) ∂Ψ ∂x Ψ variation(f) λ(Ψ(x)) ∂Ψ ∂x = 1
wonderful ⏞ ⋯ = expectation 𝔼 [f( X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn • Variable transforms are akin to importance sampling • May be di ff i cult for Bayesian posterior mean problems • Control variates may also be used • Acceptance-rejection sampling does not work well
[f( X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn • The cost to evaluate is typically , so the cost to obtain is typically f(xi ) 𝒪 (d) |μ − ̂ μn | ≤ ε 𝒪 (dε−1−δ)
[f( X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn • The cost to evaluate is typically , so the cost to obtain is typically f(xi ) 𝒪 (d) |μ − ̂ μn | ≤ ε 𝒪 (dε−1−δ) • If one can approximate by lower dimensional approximations, , then f fs : [0,1]s → ℝ μ = 𝔼 [fs1 (X1:s1 )] μ(1) + 𝔼 [fs2 (X1:s2 ) − fs1 (X1:s1 )] μ(2) + ⋯ + 𝔼 [f(X1:d ) − fsL−1 (X1:sL−1 )] μ(L) • Balance the cost to approximate well and the total cost to obtain may be as small as as 𝒪 (nl sl ) μ(s) |μ − ̂ μn | ≤ ε 𝒪 (ε−1−δ) d, ε−1 → ∞
probability density, , for , where ϱ Y = f(X) X ∼ 𝒰 [0,1]d y = f(x) ⟺ x1 = g(y; x2:d ) Then ϱ(y) = ∫ [0,1]d−1 g′  (y; x) dx If one can identify , such that g
• Discrepancy (quality) measures for • Choosing so that • Making our original problem look like the above • Ongoing research x0 , x1 , … x0 , x1 , … n |μ − ̂ μn | ≤ ε μ := expectation 𝔼 [f( X ⏟ ∼ 𝒰 ([0,1]d) )] = integral ∫ [0,1]d f(x) dx ≈ sample mean 1 n n−1 ∑ i=0 f(xi ) =: ̂ μn
• Error estimation — Art Owen — and stopping rules, especially for multilevel methods • Multi- fi delity models and Bayesian approaches — Chris Oates • Connections with machine learning — Frances Kuo • Applications beyond integration • Good software — QMCPy and others
Understand what quasi-Monte Carlo (qMC) methods are • Try qMC in place of simple (or IID) Monte Carlo • Use qMC properly • Feel free to interrupt and ask questions • Join our friendly qMC research community
ENSAE Peter Glynn, Stanford U Roshan Joseph, Georgia Tech Christiane Lemieux, U Waterloo Matt Pharr, NVIDIA Veronika Rockova, U Chicago Uros Seljak, U California, Berkeley Michaela Szölgyenyi, U Klagenfurt
M. McCourt, and A. Sorokin, QMCPy: A quasi- Monte Carlo Python library (versions 1–1.5), 2024. [DKP22] J. Dick, P. Kritzer, and F. Pillichshammer, Lattice rules: Numerical integration, approximation, and discrepancy, Springer Series in Computational Mathematics, Springer Cham, 2022. [DP10] J. Dick and F. Pillichshammer, Digital nets and sequences: Discrepancy theory and quasi-Monte Carlo integration, Cambridge University Press, Cambridge, 2010. GKS23] A. D. Gilbert, F. Y. Kuo, and I. H. Sloan, Analysis of preintegration followed by quasi-Monte Carlo integration for distribution functions and densities, SIAM J. Numer. Anal. 61 (2023), 135–166. [G13] M. Giles, Multilevel Monte Carlo methods, Monte Carlo and Quasi-Monte Carlo Methods 2012 (J. Dick, F. Y. Kuo, G. W. Peters, and I. H. Sloan, eds.), Springer Proceedings in Mathematics and Statistics, vol. 65, Springer-Verlag, Berlin, 2013. [G04] P. Glasserman, Monte Carlo methods in fi nancial engineering, Applications of Mathematics, vol. 53, Springer-Verlag, New York, 2004.
accuracy of quasi-Monte Carlo quadrature?, Monte Carlo and Quasi-Monte Carlo Methods 1998 (H. Niederreiter and J. Spanier, eds.), Springer-Verlag, Berlin, 2000, pp. 16–55. [H18] __________, The trio identity for quasi-Monte Carlo error analysis, Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Stanford, USA, August 2016 (P. Glynn and A. Owen, eds.), Springer Proceedings in Mathematics and Statistics, Springer-Verlag, Berlin, 2018, pp. 3–27. [LEPBA22] P. L'Ecuyer, F. Puchhammer, and A. Ben Abdellah, Monte Carlo and quasi-Monte Carlo density estimation via conditioning, INFORMS J. Comput. 34 (2022), no. 3, 1729–1748. [NW10] E. Novak and H. Woźniakowski, Tractability of multivariate problems Volume II: Standard information for functionals, EMS Tracts in Mathematics, no. 12, European Mathematical Society, Zürich, 2010. [LENOT24] P. L'Ecuyer, M. K. Nakayama, A. B. Owen, and B. Tu ffi n, Con fi dence intervals for randomized quasi-Monte Carlo estimators, WSC ’23: Proceedings of the Winter Simulation Conference, 2024, pp. 445–456.
M. Hyman, F. J. Hickernell, and N.~W. Hengartner, Computationally e ff i cient and error aware surrogate construction for numerical solutions of subsurface fl ow through porous media, 2024+, arXiv:2310.13765.
fjhickernell
kuehara
dmattsun
sankaku0724
hpprc
langstat
naok615
fujiihideo
ayaco0
joisino
o_ob
abap34
ikuyamada
bluesmoon
holman
chriscoyier
62gerente
addyosmani
quramy
panda_program
matthewcrist
denniskardys
rmw
sachag
eileencodes
![Integration lattices [DKP 2022] 0.00 0.25 0.50 0.75 1.00 0.0](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_26.jpg){kind=link}
![Integration lattices [DKP 2022] xi = ih/n (mod 1), i](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_27.jpg){kind=link}
![Integration lattices [DKP 2022] xi = ih/n (mod 1), i](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_28.jpg){kind=link}
![Digital nets and sequences [DP 2010] xi = i0 z0](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_34.jpg){kind=link}
![Digital nets and sequences [DP 2010] xi = i0 z0](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_35.jpg){kind=link}
![Discrepancy measures the quality of ? [H00] x0 , x1](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_41.jpg){kind=link}
![Discrepancy measures the quality of ? [H00] x0 , x1](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_42.jpg){kind=link}
![Discrepancy measures the quality of ? [H00] x0 , x1](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_43.jpg){kind=link}
![Decaying coordinate weights make integration tractable [NW10] |μ − ̂](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_56.jpg){kind=link}
![Decaying coordinate weights make integration tractable [NW10] • Problems are](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_57.jpg){kind=link}
![Ex., option pricing [G04] μ := 𝔼 [payoff(Brownian motion(t1 ,](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_82.jpg){kind=link}
![Multilevel methods reduce computational cost [G12] μ := expectation 𝔼](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_84.jpg){kind=link}
![Multilevel methods reduce computational cost [G12] μ := expectation 𝔼](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_85.jpg){kind=link}
![Multilevel methods reduce computational cost [G12] μ := expectation 𝔼](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_86.jpg){kind=link}
![Multilevel methods reduce computational cost [G12] μ := expectation 𝔼](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_87.jpg){kind=link}
![Gaussian Process Regression for PDEs with random fields [SPOMHHH24]](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_88.jpg){kind=link}
![Conditional Monte Carlo for Density Estimation [LEPBA22, GKS24] Estimate the](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_89.jpg){kind=link}
![Conditional Monte Carlo for Density Estimation [LEPBA22, GKS24] Estimate the](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_90.jpg){kind=link}
![Conditional Monte Carlo for Density Estimation [LEPBA22, GKS24] Estimate the](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_91.jpg){kind=link}
![References [QMCPy] S.-C. T. Choi, F. J. Hickernell, R. Jagadeeswaran,](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_96.jpg){kind=link}
![References [H00] F. J. Hickernell, What a ff ects the](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_97.jpg){kind=link}
![References [SPDOHHH24] A. G. Sorokin, A. Pachalieva, D. O'Malley, J.](https://files.speakerdeck.com/presentations/279b6e13d6464995bc305e7ad43bac3a/slide_98.jpg){kind=link}