The Start of Community Supported Quasi-Monte Carlo Software

Several groups of researchers in quasi-Monte Carlo software have developed code that implements their ideas, but these different code libraries do not work with each other well. This talk describes a framework and prototype of a quasi-Monte Carlo software to be supported by the research community, where different components can be replaced with alternatives seamlessly.

Transcript

1. 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
   hickernell@iit.edu 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
   

   View Slide

2. 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
   

   View Slide

3. 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
   

   View Slide

4. 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
   

   View Slide

5. 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
   

   View Slide

6. 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
   

   View Slide

7. 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
   

   View Slide

8. 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
   

   View Slide

9. 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
   

   View Slide

10. 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
    

    View Slide

11. 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
    

    View Slide

12. 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
    

    View Slide

13. 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
    

    View Slide

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
    

    View Slide

15. 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
    

    View Slide

16. 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
    

    View Slide

17. 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
    

    View Slide

18. 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
    

    View Slide

19. 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
    

    View Slide

20. 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
    

    View Slide

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

    View Slide

22. 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
    

    View Slide

23. 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
    

    View Slide