Advances and Challenges in Quasi-Monte Carlo Software Fred J. Hickernell Department of Applied Mathematics & Center for Interdisciplinary Scientific Computation Illinois Institute of Technology [email protected] mypages.iit.edu/~hickernell Thanks to Sou-Cheng Choi, Mike McCourt, Jagadeeswaran Rathinavel, Aleksei Sorokin and the rest of the QMCPy team Thanks to all of the special session speakers Thanks to the organizers, especially during these unusual times Slides at speakerdeck.com/fjhickernell/advances-and-challenges-in-quasi-monte-carlo-software MCM 2021, August 17, 2021

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References Quasi-Monte Carlo Cubature µ := [0,1]d f(x) dx integration = E[f(X)] expectation population mean ≈ 1 n n i=1 f(Xi) sample mean =: ^ µn 2/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References Quasi-Monte Carlo Cubature [0,1]d f(x) dx integration = E[f(X)] expectation population mean ≈ 1 n n i=1 f(Xi) sample mean =: ^ µn Low Discrepancy Generators producing X1, X2, . . . LD ∼ U[0, 1]d 2/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References Quasi-Monte Carlo Cubature µ := T g(t) λ(t) dt = t=T(x) [0,1]d f(x) dx integration = E[f(X)] expectation population mean ≈ 1 n n i=1 f(Xi) sample mean =: ^ µn Low Discrepancy Generators producing X1, X2, . . . LD ∼ U[0, 1]d True Measures λ used to define the original integral, such as Lebesgue or Gaussian, along with the variable transformation (importance sampling) 2/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References Quasi-Monte Carlo Cubature µ := T g(t) λ(t) dt = t=T(x) [0,1]d f(x) dx integration = E[f(X)] expectation population mean ≈ 1 n n i=1 f(Xi) sample mean =: ^ µn Low Discrepancy Generators producing X1, X2, . . . LD ∼ U[0, 1]d True Measures λ used to define the original integral, such as Lebesgue or Gaussian, along with the variable transformation (importance sampling) Integrands g used to define original integral and the final, transformed integrand, f(x) = g(T(x))λ(T(x)) T ′(x) (preferably flat) 2/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References Quasi-Monte Carlo Cubature µ := T g(t) λ(t) dt = t=T(x) [0,1]d f(x) dx integration = E[f(X)] expectation population mean ≈ 1 n n i=1 f(Xi) sample mean =: ^ µn Low Discrepancy Generators producing X1, X2, . . . LD ∼ U[0, 1]d True Measures λ used to define the original integral, such as Lebesgue or Gaussian, along with the variable transformation (importance sampling) Integrands g used to define original integral and the final, transformed integrand, f(x) = g(T(x))λ(T(x)) T ′(x) (preferably flat) Stopping Criteria that determines how large n should be to ensure that |µ − ^ µn| ⩽ ε 2/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References Software Landscape LD Sequence Generators Multi-Level, Stopping Criteria, Applications MATLAB—Sobol’ & Halton Mike Giles—Multi-Level (Quasi-)Monte Carlo in C++, MATLAB, Python, R BRODA—Sobol’ in C, MATLAB, and Excel SciPy, PyTorch—Sobol’, Halton, more in Python Guaranteed Automatic Integration Library (GAIL)—Stopping criteria in MATLAB Tensorflow—Nets, lattices & finance in Python qrng— R package, Sobol’ and Halton UQLab—Framework for Uncertainty Quantification in MATLAB Art Owen—randomized Halton in R Pieterjan Robbe—LD sequences in Julia OpenTURNS—An Open source initiative for the Treatment of Uncertainties, Risks ’N Statistics in Python Steven Johnson—Sobol’ in Julia cuRAND—Sobol’ in CUDA Pierre L’Ecuyer—LatNet Builder and Stochastic Simulation in C/C++ and Java Dirk Nuyens—Magic Point Shop and QMC4PDE in MATLAB, Python, and C++ John Burkhardt—variety in C++, Fortran, MATLAB, & Python QMCPy—Python package incorporating and connecting the work of different groups 3/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References What We Want QMC in the hands of more users Quality QMC algorithms Platform to compare performance of competing QMC algorithms Vehicle to explore QMC for new use cases SciPy/PyTorch QMCPy cuRAND MATLAB qrng MPS/QMC4PDE LatNet/SSJ UQLab BRODA Tensorflow GAIL OpenTURNS 4/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References Low Discrepancy Generators: X1 , X2 , . . . LD ∼ U[0, 1]d Recent Advances Sobol’ and Halton into SciPy. Zeroth point included; discrepancy; transformations to other distributions Lattices, digital nets, and quantitative finance into Tensorflow Niederreiter points and other “bring your own” generators in QMCPy 5/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References Low Discrepancy Generators: X1 , X2 , . . . LD ∼ U[0, 1]d Recent Advances Sobol’ and Halton into SciPy. Zeroth point included; discrepancy; transformations to other distributions Lattices, digital nets, and quantitative finance into Tensorflow Niederreiter points and other “bring your own” generators in QMCPy Challenges Need to educate against dropping, thinning, burn-in, power of 10 sample sizes [1] Need to go beyond 232 points Nested uniform scrambling of nets [2], as an alternative to linear matrix scrambling, not yet widely available Higher order nets not widely available Which generators are significantly better For which practical problems? For higher dimensions? Is there a useful hybrid of lattices and digital sequences? 5/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References True Measures: T : T → [0, 1]d T g(t) λ(t) dt = t=T(x) [0,1]d f(x) dx Variable transformations aka importance sampling Recent Advances SciPy has several QMCPy has defaults for T = Rd or [a, b]d based on the original weight function, and you can choose your own 6/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References True Measures: T : T → [0, 1]d T g(t) λ(t) dt = t=T(x) [0,1]d f(x) dx Variable transformations aka importance sampling Recent Advances SciPy has several QMCPy has defaults for T = Rd or [a, b]d based on the original weight function, and you can choose your own Challenges Importance sampling is an art Need for adaptive importance sampling Can some importance sampling make QMC competitive with MCMC for Bayeseian inference, etc. [3]? Transformations for non-rectangular T (simplex, ball, sphere) not widely available 6/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References Integrands or Use Cases T g(t) λ(t) dt = · · · Recent Advances Tensorflow has quantitative finance UQLab and OpenTURNS include LD generators for their use cases QMCPy allows you to build your own, and has a growing number of examples New use cases such as in Robbe’s talk [4] on system reliability and Ebert’s talk combining QMC with algorithmic differentiation 7/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References Integrands or Use Cases T g(t) λ(t) dt = · · · Recent Advances Tensorflow has quantitative finance UQLab and OpenTURNS include LD generators for their use cases QMCPy allows you to build your own, and has a growing number of examples New use cases such as in Robbe’s talk [4] on system reliability and Ebert’s talk combining QMC with algorithmic differentiation Challenges Our go-to use cases need refreshing and expanding Hard to swap out different LD generators for important use cases because the use cases do not have a common interface or is not always accessible Need to demonstrate connecting a package that computes the integrand values (such as an open source PDE solver) with the package that does the QMC 7/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References Stopping Criteria, n =? to Satisfy Error Tolerance Recent Advances QMCPy includes Several kinds of stopping criteria including replications, tracking Fourier coefficients [5, 6], and Bayesian credibility intervals [7] Criteria to ensure that SOL −SOLn ⩽ max(εa , εr |SOL|), where SOL is a function of more than one integral, e.g., Bayesian inference, Sobol’ indices 8/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References Stopping Criteria, n =? to Satisfy Error Tolerance Recent Advances QMCPy includes Several kinds of stopping criteria including replications, tracking Fourier coefficients [5, 6], and Bayesian credibility intervals [7] Criteria to ensure that SOL −SOLn ⩽ max(εa , εr |SOL|), where SOL is a function of more than one integral, e.g., Bayesian inference, Sobol’ indices Challenges Need better stopping criteria for higher order nets (other than replications) Need better stopping criteria for more complex problems, such as parametric integration, density estimation, Markov Chain Monte Carlo 8/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References Final Thoughts Several multi-collaborator projects are growing in functionality and user base Can our work build a shorter path between research code and production code? Can our Python packages talk to one another? Do we ignore or encourage development in R, Julia, MATLAB, ...? Need to be careful about licensing issues as we share and adapt code SciPy/PyTorch QMCPy cuRAND MATLAB qrng MPS/QMC4PDE LatNet/SSJ UQLab BRODA Tensorflow GAIL OpenTURNS 9/10

Thank you These slides are available at speakerdeck.com/fjhickernell/ advances-and-challenges-in-quasi-monte-carlo-software

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References 1. Owen, A. B. On dropping the first Sobol’ point. 2021+. https://arxiv.org/pdf/2008.08051.pdf. 2. Owen, A. B. Randomly Permuted (t, m, s)-Nets and (t, s)-Sequences. in Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing (eds Niederreiter, H. & Shiue, P. J.-S.) 106 (Springer-Verlag, New York, 1995), 299–317. 3. Zhang, K. Adaptive (Multilevel) Quasi-Monte Carlo Methods For Bayesian Inference and Uncertainty Quantification. PhD thesis (Illinois Institute of Technology, 2021+). 4. Vromans, M. Multilevel Monte Carlo for system reliability. MA thesis (Katholieke Universiteit Leuven, 2019). 5. H., F. J. & Jiménez Rugama, L. A. 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. 10/10

Background Low Discrepancy Generators True Measures Integrands Stopping Criteria Conclusion References 6. Jiménez Rugama, L. A. & 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. 7. Jagadeeswaran, R. & H., F. J. Fast Automatic Bayesian Cubature Using Lattice Sampling. Stat. Comput. 29, 1215–1229 (2019). 10/10