Upgrade to Pro — share decks privately, control downloads, hide ads and more …

The Start of Community Supported Quasi-Monte Carlo Software

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.

Fred J. Hickernell

February 27, 2019
Tweet

More Decks by Fred J. Hickernell

Other Decks in Research

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

    View full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size 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 full-size slide