# R/mlmc.R In louisaslett/mlmc: Multi-Level Monte Carlo

#### Documented in mlmc

# This code is derived and adapted from the original GPL-2 Matlab version by
# Mike Giles.  See http://people.maths.ox.ac.uk/~gilesm/mlmc/

#' Multi-level Monte Carlo estimation
#'
#' This function is the Multi-level Monte Carlo driver which will sample from the levels of user specified function.
#'
#' Multilevel Monte Carlo Method method originated in works Giles (2008) and Heinrich (1998).
#'
#' Consider  a sequence \eqn{P_0, P_1, \ldots,} which approximates \eqn{P_L} with
#' increasing accuracy, but also increasing cost, we have the simple identity
#' \deqn{E[P_L] = E[P_0] + \sum_{l=1}^L E[P_l-P_{l-1}],}
#' and therefore we can use the following unbiased estimator for \eqn{E[P_L]},
#' \deqn{N_0^{-1} \sum_{n=1}^{N_0} P_0^{(0,n)} +
#' \sum_{l=1}^L \{
#' N_l^{-1} \sum_{n=1}^{N_l} (P_l^{(l,n)} - P_{l-1}^{(l,n)})
#' \}}
#' with the inclusion of the level \eqn{l} in the superscript \eqn{(l,n)} indicating
#' that the samples used at each level of correction are independent.
#'
#' Set \eqn{C_0}, and \eqn{V_0} to be the cost and variance of one sample of \eqn{P_0},
#' and \eqn{C_l, V_l} to be the cost and variance of one sample of \eqn{P_l - P_{l-1}},
#' then the overall cost and variance of the multilevel estimator is
#' \eqn{\sum_{l=0}^L N_l C_l\}}
#' and
#' \eqn{\sum_{l=0}^L N_l^{-1} V_l},
#' respectively.
#'
#' The idea begind the method, is that provided that the product \eqn{V_l C_l}
#' decreases with \eqn{l}, i.e. the cost increases with level slower than the
#' variance decreases, then one can achieve significant computational savings,
#' which can be formalised as in Theorem 1 of Giles (2015).
#'
#' For further information on multilevel Monte Carlo methods, see the webpage
#' \url{http://people.maths.ox.ac.uk/gilesm/mlmc_community.html}
#' which lists the research groups working in the area, and their
#' main publications.
#'
#' This function is based on GPL-2 'Matlab' code by Mike Giles.
#'
#' @author Louis Aslett <[email protected]>
#' @author Mike Giles <[email protected]>
#' @author Tigran Nagapetyan <[email protected]>
#'
#' @references
#' M.B. Giles. Multilevel Monte Carlo path simulation. \emph{Operations Research}, 56(3):607-617, 2008.
#'
#' M.B. Giles. Multilevel Monte Carlo methods. \emph{Acta Numerica}, 24:259-328, 2015.
#'
#' S. Heinrich. Monte Carlo complexity of global solution of integral equations. \emph{Journal of Complexity}, 14(2):151-175, 1998.
#'
#' @param Lmin the minimum level of refinement.  Must be \eqn{\ge 2}.
#' @param Lmax the maximum level of refinement.  Must be \eqn{\ge} Lmin.
#' @param N0 initial number of samples which are used for the first 3 levels and
#'   for any subsequent levels which are automatically added.  Must be
#'   \eqn{> 0}.
#' @param eps the target accuracy of the estimate.  Must be \eqn{> 0}.
#' @param mlmc_l a user supplied function which provides the estimate for level
#'   l
#' @param alpha the weak error, \eqn{O(2^{-alpha*l})}.  If \code{NA} then
#'   \code{alpha} will be estimated.
#' @param beta the variance, \eqn{O(2^{-beta*l})}.  If \code{NA} then
#'   \code{beta} will be estimated.
#' @param gamma the sample cost, \eqn{O(2^{gamma*l})}.  Must be \eqn{> 0}.
#' @param parallel if an integer is supplied, R will fork \code{parallel} parallel
#'   processes an compute each level estimate in parallel.
#' @param ... additional arguments which are passed on when the user supplied
#'   \code{mlmc_l} function is called
#'
#' @return A list containing: \describe{
#'   \item{\code{P}}{The MLMC estimate;}
#'   \item{\code{Nl}}{A vector of the number of samples performed on each level.}
#' }
#'
#' @examples
#' mlmc(2, 6, 1000, 0.01, opre_l, gamma=1, option=1)
#'
#' mlmc(2, 10, 1000, 0.01, mcqmc06_l, gamma=1, option=1)
#'
#' @importFrom parallel mcmapply
#' @importFrom stats lm
#' @export
mlmc <- function(Lmin, Lmax, N0, eps, mlmc_l, alpha=NA, beta=NA, gamma, parallel=NA, ...) {
# check parameters are acceptable
if(Lmin<2) {
stop("Lmin must be >= 2.")
}
if(Lmax<Lmin) {
stop("must have Lmax >= Lmin.")
}
if(N0<=0 || eps<=0 || gamma <= 0){
stop("N0, eps and gamma must all be greater than zero.")
}
if(!is.na(alpha) && alpha<=0) {
stop("if specified, alpha must be greater than zero.  Set alpha to NA to automatically estimate.")
}
if(!is.na(beta) && beta<=0) {
stop("if specified, beta must be greater than zero.  Set beta to NA to automatically estimate.")
}

# initialise the MLMC run
est.alpha <- is.na(alpha)
alpha <- ifelse(is.na(alpha), 0, alpha)
est.beta  <- is.na(beta)
beta <- ifelse(is.na(beta), 0, beta)

theta <- 0.25

L <- Lmin

Nl <- rep(0, L+1)
suml <- matrix(0, nrow=2, ncol=L+1)
dNl <- rep(N0, L+1)

while(sum(dNl) > 0) {
# update sample sums from each level
if(is.na(parallel)) {
for(l in 0:L) {
if(dNl[l+1] > 0) {
sums        <- mlmc_l(l, dNl[l+1], ...)
Nl[l+1]     <- Nl[l+1] + dNl[l+1]
suml[1,l+1] <- suml[1,l+1] + sums[1]
suml[2,l+1] <- suml[2,l+1] + sums[2]
}
}
} else if(is.numeric(parallel)) {
par.vars <- data.frame(l=0:L, dNl=dNl)[!(dNl==0),]
sums <- mcmapply(function(l, dNl, ...) {
mlmc_l(l, dNl, ...)
}, l = par.vars$l, dNl = par.vars$dNl, ..., mc.preschedule = FALSE, mc.cores = parallel)
Nl <- Nl+dNl
suml[,!(dNl==0)] <- suml[,!(dNl==0)] + sums[1:2,]
}

# compute absolute average and variance
ml <- abs(    suml[1,]/Nl)
Vl <- pmax(0, suml[2,]/Nl - ml^2)

# fix to cope with possible zero values for ml and Vl
# (can happen in some applications when there are few samples)
for(l in 3:(L+1)) {
ml[l] <- max(ml[l], 0.5*ml[l-1]/2^alpha)
Vl[l] <- max(Vl[l], 0.5*Vl[l-1]/2^beta)
}

# estimate alpha and beta by regression if specified
if(est.alpha) {
alpha <- max(0.5, -lm(y~x, data.frame(y=log2(ml[-1]), x=1:L))$coefficients["x"]) } if(est.beta) { beta <- max(0.5, -lm(y~x, data.frame(y=log2(Vl[-1]), x=1:L))$coefficients["x"])
}
# choose optimal number of additional samples
Cl  <- 2^(gamma*(0:L))
Ns  <- ceiling(sqrt(Vl/Cl) * sum(sqrt(Vl*Cl)) / ((1-theta)*eps^2))
dNl <- pmax(0, Ns-Nl)
# if (almost) converged, estimate remaining error and decide
# whether a new level is required
if( !any(dNl > 0.01*Nl) ) {
rem <- ml[L+1] / (2^alpha - 1)

if(rem > sqrt(theta)*eps) {
if(L==Lmax) {
warning("failed to achieve weak convergence.")
} else {
L       <- L+1
Vl[L+1] <- Vl[L] / 2^beta
Nl[L+1] <- 0
suml <- cbind(suml, 0)

Cl <- 2^(gamma*(0:L))
Ns  <- ceiling(sqrt(Vl/Cl) * sum(sqrt(Vl*Cl)) / ((1-theta)*eps^2))
dNl <- pmax(0, Ns-Nl)
}
}
}
}

# finally, evaluate multilevel estimator
P <- sum(suml[1,]/Nl)
list(P=P, Nl=Nl)
}
louisaslett/mlmc documentation built on May 21, 2017, 2:39 p.m.