Homework 2

Due: end of week 7

Monte Carlo, Multi-level Monte Carlo, and Multi-order Monte Carlo Sampling

Note: The codes for this homework are available in the public Bitbucket repository: https://bitbucket.org/motamed/uq2017 within the subdirectory hw2. You can clone the whole repository to your local repository, for instance by typing in a shell terminal:

$ git clone https://motamed@bitbucket.org/motamed/uq2017.git

Description: In this assignment you will consider a stochastic ODE problem and compare the performance of the following three Monte Carlo sampling techniques:

  • Classical Monte Carlo (MC)
  • Multi-level Monte Carlo (MOMC)
  • Multi-order Monte Carlo (MOMC)

For a vector-valued stochastic process \({\bf u}(t,y) = [u_1(t,y),u_2(t,y)]^{\top}\), consider the stochastic ODE problem:

\begin{eqnarray*} \left\{ \begin{array}{l} {\bf u}_t(t,y) = A(y) \, {\bf u}(t,y), \qquad t \in [0,T]\\ \\ {\bf u}(0,y) = [1, 0]^{\top}, \hskip 2.2cm t=0 \end{array} \right. \end{eqnarray*}

where \(A(y)\) is a stochastic matrix given in terms of one uniform random variable \(y\):

$$ A(y) = \left[ \begin{array}{c c} 0 & \, \, 1+y\\ -1-y & \, \, 0 \end{array} \right], \qquad y \sim U[0,1]. $$

Consider the following stochastic quantity of interest (QoI):

$$Q(y) = u_1(T, y).$$

We want to approximate the expected value of \(Q\) (i.e. \({\mathbb E}[Q]\)) by MC, MLMC, and MOMC.

The stochastic problem can indeed be solved analytically. The analytical solution is given by \(u_1 = \sin (t (1+y))\) and \(u_2 = \cos (t (1+y))\). Hence, the exact expected value of the QoI is

$${\mathbb E}[Q] = (\cos(T)-\cos(2T))/T.$$

This value will be used to measure the total error. Throughout this homework, we take \(T=1\).

  1. In the case of MLMC and MOMC, consider a failure probability for the statistical error

    $$P \Bigl( \varepsilon_{II} \ > \ C_{\alpha} \, \bigl(\sum_{\ell =0}^L \frac{V_{\ell}}{M_{\ell}} \bigr)^{1/2} \ > \ \theta \,\varepsilon_{TOL} \Bigr) \ \le \ \alpha.$$

    Here, \(\varepsilon_{II}\) is the statistical error in either MLMC or MOMC, \(C_{\alpha}\) is the confidence parameter, \(V_{\ell}\) is the variance of the difference \(Q_{\ell} - Q_{\ell-1}\) at level \(\ell\) (setting \(Q_{-1} = 0\)), and \(M_{\ell}\) is the number of samples at level \(\ell\). Find an appropriate confidence parameter \(C_{\alpha}\) for two cases: \(\alpha = 5 \%\) and \(\alpha =0.1 \%\). We will use these two values of \(C_{\alpha}\) and observe their effect on the error later.

  2. Download the following MATLAB files from the subdirectory hw2 in the public Bitbucket repository https://bitbucket.org/motamed/uq2017:

    • MC codes: mc_conv.m and mc.m
    • MLMC codes: mlmc_conv.m and mlmc.m and mlmc_ode.m
    • MOMC codes: momc_conv.m and momc.m and momc_ode.m

    All three methods will call and use the MATLAB file ode_taylor.m, where the deterministic solver (the \(q\) -th order accurate Taylor’s method) is implemented. In total you will need to download nine MATLAB files. For the following taks, you will need to modify only the three main codes: mc_conv.m, mlmc_conv.m, and momc_conv.m. Do not modify the other six files. Each of these three files consists of 3 parts: 1) Input parameters; 2) Computations; 3) Plots.

    Consider a set of tolerances \(\varepsilon_{TOL} = 20 \times 10^{-4}, 19 \times 10^{-4}, 18 \times 10^{-4}, \dotsc, 1 \times 10^{-4}\) for all three main codes. In the codes this is written as Eps=(20:-1:1)*(1e-4)}. Use the following parameter values:

    • mc_conv.m:

      \(q=2\) (a fixed order of accuracy of deterministic solver)

    • mlmc_conv.m:

      \(q=2\) (a fixed order of accuracy of deterministic solver)

      \(h_0 = 0.5\) (time step for deterministic solver at level \(\ell = 0\))

      \(\beta = 2\) (mesh refinement ratio, with the mesh hierarchy \(h_{\ell} = h_0 \, \beta^{-\ell}\))

      \(\theta = 0.9\) (splitting parameter)

      \(C_{\alpha} =\) the value corresponding to \(\alpha = 5 \%\) obtained in task 1

    • momc_conv.m

      \(h=0.25\) (a fixed time step for deterministic solver)

      \(q_0 = 2\) (order of accuracy of deterministic solver at level \(\ell = 0\))

      \(\beta = 2\) (order increment parameter, with the order hierarchy \(q_{\ell} = q_0 + \beta \ \ell\))

      \(\theta = 0.9\) (splitting parameter)

      \(C_{\alpha} =\) the value corresponding to \(\alpha = 5 \%\) obtained in task 1.

    Generate two figures:

    • Figure 1: plot the total errors in MC, MLMC, and MOMC (y-axis) versus tolerances \(\varepsilon_{TOL}\) (x-asis) in the same figure. Use loglog to get a logarithmic scale on both axes. Use different markers/colors for the errors generated by different methods. The plot should have labels on each axis and a legend to describe each error. Use large fonts. Also plot the line \(\varepsilon_{TOL}\) versus \(\varepsilon_{TOL}\). From your figure you should be able to see if the total error generated by the three methods remain below \(\varepsilon_{TOL}\) or not. If not, why? Can you explain?
    • Figure 2: plot the total computational cost of MC, MLMC, and MOMC (y-axis) versus tolerances \(\varepsilon_{TOL}\) (x-asis) in the same figure. Use loglog to get a logarithmic scale on both axes. Use different markers/colors for the errors generated by different methods. The plot should have labels on each axis and a legend. Use large fonts. From your figure you should be able to see which method is faster. Comment and discuss on the convergence rate of the computational costs with respect to tolerance \(\varepsilon_{TOL}\).

  1. Repeat task 2 with \(C_{\alpha} =\) the value corresponding to \(\alpha = 0.1 \%\) obtained in task 1. Describe and explain the change that you may observe compared to task 2.

  1. Next, consider only mlmc_conv.m and momc_conv.m with a tolerance \(\varepsilon_{TOL} = 10^{-3}\) and \(C_{\alpha} =\) the value corresponding to \(\alpha = 0.1 \%\).
    1. Numerically play with parameters \(q, h_0, \beta, \theta\) in mlmc_conv.m and empirically find a combination of parameters that produces the smallest possible computational cost of MLMC. Discuss!
    2. Numerically play with parameters \(h, q_0, \beta, \theta\) in mlmc_conv.m and empirically find a combination of parameters that produces the smallest possible computational cost of MOMC. Discuss!
    3. Compare the minimum possible costs of MLMC and MOMC.