The Start of Community Supported 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 the conference organizers, the GAIL team NSF-DMS-1522687 and NSF-DMS-1638521 (SAMSI) SIAM CSE, February 27, 2019

Motivation Problem Formulation MATLAB Prototype Next References Common Questions Where can I find quality, free quasi-Monte Carlo (qMC) software? Where can I find use cases that illustrate how to employ qMC methods? How can I try that qMC method developed by X’s group for my problem? How can I try my qMC method on the example shown by Y’s group? How can my student get results without writing code from scratch? How can the code that I use benefit from recent developments? How can my work receive wider recognition? 2/14

Motivation Problem Formulation MATLAB Prototype Next References What Is Available Now Does Not Communicate with Each Other John Burkhardt Variety of qMC Software in C++, Fortran, MATLAB, and Python, people.sc.fsu.edu/~jburkardt/ Mike Giles Multi-Level Monte Carlo Software in C++, MATLAB, Python, and R, people.maths.ox.ac.uk/gilesm/mlmc/ Fred Hickernell Guaranteed Automatic Integration Library (GAIL) in MATLAB, gailgithub.github.io/GAIL_Dev/ Stephen Joe & Frances Kuo Sobol’ generators in C++, Generating vectors for lattices, web.maths.unsw.edu.au/~fkuo/ Sergei Kucherenko Sobol’ generators and related software www.broda.co.uk Pierre L’Ecuyer Random number generators, Stochastic Simulation, Lattice Builder in C/C++ and Java, simul.iro.umontreal.ca Dirk Nuyens Magic Point Shop, QMC4PDE, etc. in MATLAB, Python, and C++, people.cs.kuleuven.be/~dirk.nuyens/ Art Owen Various code, statweb.stanford.edu/~owen/code/ Christoph Schwab Partial differential equations with random coefficients MATLAB Sobol’ and Halton sequences Python Sobol’ and Halton sequences R randtoolbox Sobol’, lattice, and Halton sequences 3/14

Motivation Problem Formulation MATLAB Prototype Next References Can We Have Community qMC Software that Is ... Developed and supported by multiple research groups Used beyond the research groups that develop it Recognized as the standard Easy for the non-expert to use 4/14

Motivation Problem Formulation MATLAB Prototype Next References The Integration Problem ż T g(t) λ(t) dt looooooomooooooon original problem = ż X f(x) ρ(x) dx = ż X f(x) ν(dx) loooooooooooooooooooomoooooooooooooooooooon convenient form « ż X f(x) ^ νn(dx) = n ÿ i=1 f(xi)wi looooooooooooooooomooooooooooooooooon (quasi-)Monte Carlo approximation , n =? lo omo on stopping criterion 5/14

Motivation Problem Formulation MATLAB Prototype Next References The Integration Problem ż T g(t) λ(t) dt looooooomooooooon original problem = ż X f(x) ρ(x) dx = ż X f(x) ν(dx) loooooooooooooooooooomoooooooooooooooooooon convenient form « ż X f(x) ^ νn(dx) = n ÿ i=1 f(xi)wi looooooooooooooooomooooooooooooooooon (quasi-)Monte Carlo approximation , n =? lo omo on stopping criterion High dimensional integration problems arise in Financial risk management Statistical physics Sensitivity indices Uncertainty quantification for PDEs Bayesian inference g : T Ñ R = original integrand, e.g., option payoff λ = original weight, e.g., PDF for discrete Brownian motion 5/14

Motivation Problem Formulation MATLAB Prototype Next References The Integration Problem ż T g(t) λ(t) dt looooooomooooooon original problem = ż X f(x) ρ(x) dx = ż X f(x) ν(dx) loooooooooooooooooooomoooooooooooooooooooon convenient form « ż X f(x) ^ νn(dx) = n ÿ i=1 f(xi)wi looooooooooooooooomooooooooooooooooon (quasi-)Monte Carlo approximation , n =? lo omo on stopping criterion ν = probability measure that we can easily sample from with density ρ e.g., standard uniform or standard normal φ : X Ñ T = change of variables e.g., inverse Gaussian, Brownian bridge, PCA, importance sampling f : X Ñ R = integrand after an appropriate change of variables f(x) = g(φ(x)) det Bφ(x) Bx λ(φ(x)) ρ(x) 5/14

Motivation Problem Formulation MATLAB Prototype Next References The Integration Problem ż T g(t) λ(t) dt looooooomooooooon original problem = ż X f(x) ρ(x) dx = ż X f(x) ν(dx) loooooooooooooooooooomoooooooooooooooooooon convenient form « ż X f(x) ^ νn(dx) = n ÿ i=1 f(xi)wi looooooooooooooooomooooooooooooooooon (quasi-)Monte Carlo approximation , n =? lo omo on stopping criterion ^ νn = discrete probability distribution « ν e.g., wi = 1/n with xi being IID or low discrepancy 5/14

Motivation Problem Formulation MATLAB Prototype Next References The Integration Problem ż T g(t) λ(t) dt looooooomooooooon original problem = ż X f(x) ρ(x) dx = ż X f(x) ν(dx) loooooooooooooooooooomoooooooooooooooooooon convenient form « ż X f(x) ^ νn(dx) = n ÿ i=1 f(xi)wi looooooooooooooooomooooooooooooooooon (quasi-)Monte Carlo approximation , n =? lo omo on stopping criterion What n guarantees ż X f(x) ν(dx) ´ n ÿ i=1 f(xi)wi ď ε with high confidence ? H., F. J. et al. Guaranteed Conservative Fixed Width Confidence 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, H., F. J. & Jiménez Rugama, Ll. 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 (Springer-Verlag, Berlin, 2016), 367–383, Jiménez Rugama, Ll. 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 (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. 5/14

Motivation Problem Formulation MATLAB Prototype Next References The Integration Problem ż T g(t) λ(t) dt looooooomooooooon original problem = ż X f(x) ρ(x) dx = ż X f(x) ν(dx) loooooooooooooooooooomoooooooooooooooooooon convenient form « ż X f(x) ^ νn(dx) = n ÿ i=1 f(xi)wi looooooooooooooooomooooooooooooooooon (quasi-)Monte Carlo approximation , n =? lo omo on stopping criterion The integral may be infinite dimensional : lim dÑ∞ ż Td gd(t) λd(t) dt = lim dÑ∞ ż Xd fd(x) νd(dx) = ż Xd1 fd1 (x) νd1 (dx) + ż Xd2 [fd2 (x) ´ fd1 (x)] νd2 (dx) + ¨ ¨ ¨ « n1 ÿ i=1 fd1 (xi,d1 )wi,d1 + n2 ÿ i=1 [fd2 (xi,d2 ) ´ fd1 (xi,d2 )]wi,d2 + ¨ ¨ ¨ + nL ÿ i=1 [fdL (xi,dL ) ´ fdL´1 (xi,dL )]wi,dL Giles, M. Multilevel Monte Carlo Methods. Acta Numer. 24, 259–328 (2015), Wasilkowski, G. W. On tractability of linear tensor product problems for ∞-variate classes of functions. J. Complexity 29, 351–369 (2013). 5/14

Motivation Problem Formulation MATLAB Prototype Next References The Integration Problem ż T g(t) λ(t) dt looooooomooooooon original problem = ż X f(x) ρ(x) dx = ż X f(x) ν(dx) loooooooooooooooooooomoooooooooooooooooooon convenient form « ż X f(x) ^ νn(dx) = n ÿ i=1 f(xi)wi looooooooooooooooomooooooooooooooooon (quasi-)Monte Carlo approximation , n =? lo omo on stopping criterion We need to be able to code perhaps only one of the following: Use cases g and λ Variable transformations ϕ : X Ñ T Discrete distributions ^ νn Stopping criteria to determine n and use what others have coded for the rest. 5/14

Motivation Problem Formulation MATLAB Prototype Next References Key Elements ż T g(t) λ(t) dt looooooomooooooon original problem = ż X f(x) ρ(x) dx = ż X f(x) ν(dx) loooooooooooooooooooomoooooooooooooooooooon convenient form « ż X f(x) ^ νn(dx) = n ÿ i=1 f(xi)wi looooooooooooooooomooooooooooooooooon (quasi-)Monte Carlo approximation , n =? lo omo on stopping criterion Use cases g and λ fun and measure classes Variable transformations ϕ : X Ñ T fun method Discrete distributions ^ νn discreteDistribution class Stopping criteria to determine n stopCriterion and accumData classes Software at github.com/QMCSoftware/QMCSoftware 6/14

Motivation Problem Formulation MATLAB Prototype Next References Keister’s Function—a Subclass of fun classdef KeisterFun < fun % Specify and generate values f(x) = πd/2 cos( x ) for x P Rd % The standard example integrates the Keister function with respect to an ... IID Gaussian distribution with variance 1/2 % B. D. Keister, Multidimensional Quadrature Algorithms, Computers in Physics, 10, pp. 119-122, 1996. methods function y = g(obj, x, coordIndex) %if the nominalValue = 0, this is efficient normx2 = sum(x.*x,2); nCoordIndex = numel(coordIndex); if (nCoordIndex ‰ obj.dimension) && (obj.nominalValue ‰ 0) normx2 = normx2 + ... (obj.nominalValue.^2) * (obj.dimension - nCoordIndex); end y = (pi.^(nCoordIndex/2)).* cos(sqrt(normx2)); end end end 7/14

Motivation Problem Formulation MATLAB Prototype Next References Keister Integration Example %% Integrating a function using our community QMC framework % An example with Keister's function integrated with respect to the uniform % distribution over the unit cube dim = 3; %dimension for the Keister example measureObj = IIDZMeanGaussian(measure, 'dimension', {dim}, 'variance', {1/2}); distribObj = IIDDistribution('trueD',stdGaussian(measure,'dimension', ... {dim})); %IID sampling stopObj = CLTStopping; %stopping criterion for IID sampling using the ... Central Limit Theorem [sol, out] = integrate(KeisterFun, measureObj, distribObj, stopObj) >> IntegrationExample sol = 2.1729 out = timeUsed: 0.0720 nSamplesUsed: 463985 confidInt: [2.1629 2.1829] 8/14

Motivation Problem Formulation MATLAB Prototype Next References Keister Integration Example stopObj.absTol = 1e-3; %decrease tolerance [sol, out] = integrate(KeisterFun, measureObj, distribObj, stopObj) sol = 2.1689 out = timeUsed: 5.2348 nSamplesUsed: 46297044 confidInt: [2.1679 2.1699] 9/14

Motivation Problem Formulation MATLAB Prototype Next References Asian Call Example, Low d %A multilevel example of Asian option pricing stopObj.absTol = 0.01; %increase tolerance stopObj.nMax = 1e8; %pushing the sample budget back up measureObj = BrownianMotion(measure,'timeVector', {1/4:1/4:1}); OptionObj = AsianCallFun(measureObj); %4 time steps distribObj = IIDDistribution('trueD',stdGaussian(measure,'dimension', ... {4})); %IID sampling [sol, out] = integrate(OptionObj, measureObj, distribObj, stopObj) sol = 6.1758 out = timeUsed: 1.1388 nSamplesUsed: 5910618 confidInt: [6.1658 6.1858] 10/14

Motivation Problem Formulation MATLAB Prototype Next References Asian Call Example, High d measureObj = BrownianMotion(measure,'timeVector', {1/64:1/64:1}); OptionObj = AsianCallFun(measureObj); %64 time steps distribObj = IIDDistribution('trueD',stdGaussian(measure,'dimension', ... {64})); %IID sampling [sol, out] = integrate(OptionObj, measureObj, distribObj, stopObj) sol = 6.2101 out = timeUsed: 31.4067 nSamplesUsed: 6073589 confidInt: [6.2001 6.2201] 11/14

Motivation Problem Formulation MATLAB Prototype Next References Asian Call Example, Multi-Level measureObj = BrownianMotion(measure,'timeVector', {1/4:1/4:1, 1/16:1/16:1, ... 1/64:1/64:1}); OptionObj = AsianCallFun(measureObj); %multi-level distribObj = IIDDistribution('trueD',stdGaussian(measure,'dimension', {4, ... 16, 64})); %IID sampling [sol, out] = integrate(OptionObj, measureObj, distribObj, stopObj) sol = 6.2064 out = timeUsed: 2.2871 nSamplesUsed: [7112556 828225 131642] confidInt: [6.1973 6.2155] 12/14

Motivation Problem Formulation MATLAB Prototype Next References Where We Are and Where We Need to Go Identified available software whose features we want to include Constructed and prototyped in MATLAB the primary software elements Undergrad Aleksei Sorokin is building a Python version 13/14

Motivation Problem Formulation MATLAB Prototype Next References Where We Are and Where We Need to Go Identified available software whose features we want to include Constructed and prototyped in MATLAB the primary software elements Undergrad Aleksei Sorokin is building a Python version Need to adapt some of the other available software to fit the prototype Need to identify which language(s) we will support Need to gain enough traction to attract more contributors, would you like to join us? 13/14

Thank you Slides available at speakerdeck.com/fjhickernell/ the-start-of-community-supported-quasi-monte-carlo-software Software available at github.com/QMCSoftware/QMCSoftware

Motivation Problem Formulation MATLAB Prototype Next References H., F. J., Jiang, L., Liu, Y. & Owen, A. B. Guaranteed Conservative Fixed Width Confidence 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. & Jiménez Rugama, Ll. 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 (Springer-Verlag, Berlin, 2016), 367–383. Jiménez Rugama, Ll. 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 (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 14/14

Motivation Problem Formulation MATLAB Prototype Next References birthday of Ian Sloan (eds Dick, J., Kuo, F. Y. & Woźniakowski, H.) (Springer-Verlag, 2018), 597–619. Giles, M. Multilevel Monte Carlo Methods. Acta Numer. 24, 259–328 (2015). Wasilkowski, G. W. On tractability of linear tensor product problems for ∞-variate classes of functions. J. Complexity 29, 351–369 (2013). 14/14