Elliptic PDE¶
import copy
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.special import gamma, kv
from qmcpy.integrand import Integrand
from qmcpy.util.data import Data
import qmcpy as qp
# matplotlib options
rc_fonts = {
"text.usetex": True,
"font.size": 14,
"mathtext.default": "regular",
"axes.titlesize": 14,
"axes.labelsize": 14,
"legend.fontsize": 14,
"xtick.labelsize": 12,
"ytick.labelsize": 12,
"figure.titlesize": 16,
"font.family": "serif",
"font.serif": "computer modern roman",
}
mpl.rcParams.update(rc_fonts)
# set random seed for reproducibility
np.random.seed(9999)
We will apply various multilevel Monte Carlo and multilevel quasi-Monte Carlo methods to approximate the expected value of a quantity of interest derived from the solution of a one-dimensional partial differential equation (PDE), where the diffusion coefficient of the PDE is a lognormal Gaussian random field. This example problem serves as an important benchmark problem for various methods in the uncertainty quantification and quasi-Monte Carlo literature. It is often referred to as the fruitfly problem of uncertainty quantification.
1. Problem definition¶
Let $Q$ be a quantity of interest derived from the solution $u(x, \omega)$ of the one-dimensional partial differential equation (PDE)
$$-\frac{d}{dx}\bigg(a(x, \omega) \frac{d}{dx} u(x, \omega) \bigg) = f(x), \quad 0 \leq x \leq 1,$$ $$u(0, \cdot) = u_0,$$ $$u(1, \cdot) = u_1.$$
The notation $u(x, \omega)$ is used to indicate that the solution depends on both the spatial variable $x$ and the uncertain parameter $\omega$. This uncertainty is present because the diffusion coefficient, $a(x, \omega)$, is given by a lognormal Gaussian random field with given covariance function. A common choice for the covariance function is the so-called Matérn covariance function
$$c(x, y) = \hat{c}\bigg(\frac{\|x - y\|}{\lambda}\bigg)\quad \quad \hat{c}(r) = \frac{2^{1-\nu}}{\Gamma(\nu)} r^\nu K_\nu(r)$$
with $\Gamma$ the gamma function and $K_\nu$ the Bessel function of the second kind. This covariance function has two parameters: $\lambda$, the length scale, and $\nu$, the smoothness parameter.
We begin by defining the Matérn covariance function Matern(x, y):
def Matern(x, y, smoothness=1, lengthscale=1):
distance = abs(x - y)
r = distance/lengthscale
prefactor = 2**(1-smoothness)/gamma(smoothness)
term1 = r**smoothness
term2 = kv(smoothness, r)
np.fill_diagonal(term2, 1)
cov = prefactor * term1 * term2
np.fill_diagonal(cov, 1)
return cov
Let's take a look at the covariance matrix obtained by evaluating the covariance function in n=25 equidistant points in [0, 1].
def get_covariance_matrix(pts, smoothness=1, lengthscale=1):
X, Y = np.meshgrid(pts, pts)
return Matern(X, Y, smoothness, lengthscale)
n = 25
pts = np.linspace(0, 1, num=n)
fig, ax = plt.subplots(figsize=(6, 5))
ax.pcolor(get_covariance_matrix(pts).T)
ax.invert_yaxis()
ax.set_ylabel(r"$x$")
ax.set_xlabel(r"$y$")
ax.set_title(f"Matern kernel")
plt.show()
A lognormal Gaussian random field $a(x, \omega)$ can be expressed as $a(x, \omega) = \exp(b(x, \omega))$, where $b(x, \omega)$ is a Gaussian random field. Samples of the Gaussian random field $b(x, \omega)$ can be computed from a factorization of the covariance matrix. Specifically, suppose we have a spectral (eigenvalue) expansion of the covariance matrix $C$ as
$$C = V W V^T,$$
then samples of the Gaussian random field can be computed as
$$\boldsymbol{b} = S \boldsymbol{x},$$
where $S = V W^{1/2}$ and $\boldsymbol{x}$ is a vector of standard normal independent and identically distributed random variables. This is easy to see, since
$$\mathbb{E}[\boldsymbol{b}] = \mathbb{E}[S \boldsymbol{x}] = S\mathbb{E}[\boldsymbol{x}] = \boldsymbol{0}$$ $$\mathbb{E}[\boldsymbol{b} \boldsymbol{b}^T] = \mathbb{E}[S \boldsymbol{x} \boldsymbol{x}^T S^T] = S \mathbb{E}[\boldsymbol{x} \boldsymbol{x}^T] S^T = SS^T = VWV^T = C.$$
First, let's compute an eigenvalue decomposition of the covariance matrix.
def get_eigenpairs(n, smoothness=1, lengthscale=1):
h = 1/(n-1)
pts = np.linspace(h/2, 1 - h/2, num=n - 1)
cov = get_covariance_matrix(pts, smoothness, lengthscale)
w, v = np.linalg.eig(cov)
# ensure all eigenvectors are correctly oriented
for col in range(v.shape[1]):
if v[0, col] < 0:
v[:, col] *= -1
return pts, w, v
Next, we plot the eigenfunctions for different values of $n$, the number of grid points.
n = [8, 16, 32] # list of number of gridpoints to plot
m = 5 # number of eigenfunctions to plot
fig, axes = plt.subplots(m, len(n), figsize=(8, 6))
for j, k in enumerate(n):
x, w, v = get_eigenpairs(k)
for i in range(m):
axes[i, j].plot(x, v[:, i])
axes[i, j].set_xlim(0, 1)
axes[i, j].get_yaxis().set_ticks([])
if i < m - 1:
axes[i, j].get_xaxis().set_ticks([])
if i == 0:
axes[i, j].set_title(r"$n = " + repr(k) + r"$")
plt.tight_layout()
With this eigenvalue decomposition, we can compute samples of the Gaussian random field $b(x, \omega)$, and hence, also of the lognormal Gaussian random field $a(x, \omega) = \exp(b(x, \omega))$, since $\boldsymbol{b} = V W^{1/2} \boldsymbol{x}$.
def evaluate(w, v, y=None):
if y is None:
y = np.random.randn(len(w) - 1)
m = len(y)
return v[:, :m] @ np.diag(np.sqrt(w[:m])) @ y
Let's plot a couple of realizations of the Gaussian random field $b(x, \omega)$.
n = 64
x, w, v = get_eigenpairs(n)
fig, ax = plt.subplots(figsize=(6, 4))
for _ in range(10):
ax.plot(x, evaluate(w, v))
ax.set_xlim(0, 1)
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$\log(a(x, \cdot))$")
plt.show()
Now that we are able to compute realizations of the Gaussian random field, a next step is to compute a numerical solution of the PDE
$$-\frac{d}{dx}\bigg(a(x, \omega) \frac{d}{dx} u(x, \omega) \bigg) = f(x), \quad 0 \leq x \leq 1.$$
Using a straightforward finite-difference approximation, it is easy to show that the numerical solution $\boldsymbol{u}$ is the solution of a tridiagonal system. The solutions of such a tridiagonal system can be easily obtained in $O(n)$ (linear) time using the tridiagonal matrix algorithm (also known as the Thomas algorithm). More details can be found here.
def thomas(a, b, c, d):
n = len(b)
x = np.zeros(n)
for i in range(1, n):
w = a[i-1]/b[i-1]
b[i] -= w*c[i-1]
d[i] -= w*d[i-1]
x[n-1] = d[n-1]/b[n-1]
for i in reversed(range(n-1)):
x[i] = (d[i] - c[i]*x[i+1])/b[i]
return x
For the remainder of this notebook, we will assume that the source term $f(x)=1$ and Dirichlet boundary conditions $u(0) = u(1) = 0$.
def pde_solve(a):
n = len(a)
b = np.full((n-1, 1), 1/n**2)
x = thomas(-a[1:n-1], a[:n-1] + a[1:], -a[1:n-1], b)
return np.insert(x, [0, n-1], [0, 0])
Let's compute and plot a couple of solutions $u(x, \omega)$.
n = 64
_, w, v = get_eigenpairs(n)
x = np.linspace(0, 1, num=n)
fig, ax = plt.subplots(figsize=(6, 4))
for _ in range(10):
a = np.exp(evaluate(w, v))
u = pde_solve(a)
ax.plot(x, u)
ax.set_xlim(0, 1)
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$u(x, \cdot)$")
plt.show()
In the multilevel Monte Carlo method, we will rely on the ability to generate "correlated" solutions of the PDE with varying mesh sizes. Such a correlated solutions can be used as efficient control variates to reduce the variance (or statistical error) in the approximation of the expected value $\mathbb{E}[Q]$. Since we are using a factorization of the covariance matrix to generate realizations of the Gaussian random field, it is quite easy to obtain correlated samples: when sampling from the "coarse" solution level, use the same set of random numbers used to sample from the "fine" solution level, but truncated to the appropriate size. Since the eigenvalue decomposition will reveal the most important modes in the covariance matrix, that same eigenvalue decomposition on a "coarse" approximation level will contain the same eigenfunctions, represented on the coarse grid. Let's illustrate this property on an example using n = 16 grid points for the fine solution level and n = 8 grid points for the coarse solution level.
nf = 16
nc = nf//2
_, wf, vf = get_eigenpairs(nf)
_, wc, vc = get_eigenpairs(nc)
xf = np.linspace(0, 1, num=nf)
xc = np.linspace(0, 1, num=nc)
fig, ax = plt.subplots(figsize=(6, 4))
for _ in range(10):
yf = np.random.randn(nf - 1)
af = np.exp(evaluate(wf, vf, y=yf))
uf = pde_solve(af)
ax.plot(xf, uf)
yc = yf[:nc - 1]
ac = np.exp(evaluate(wc, vc, y=yc))
uc = pde_solve(ac)
ax.plot(xc, uc, color=ax.lines[-1].get_color(), linestyle="dashed", dash_capstyle="round")
ax.set_xlim(0, 1)
ax.set_ylim(bottom=0)
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$u(x, \cdot)$")
plt.show()
The better the coarse solution matches the fine grid solution, the more efficient the multilevel methods in Section 3 will perform.
2. Single-level methods¶
Let's begin by using the single-level Monte Carlo and quasi-Monte Carlo methods to compute the expected value $\mathbb{E}[Q]$. As quantity of interest $Q$ we take the solution of the PDE at $x=1/2$, i.e., $Q = u(1/2, \cdot)$.
To integrate the elliptic PDE problem into QMCPy, we construct a simple class as follows:
class EllipticPDE(Integrand):
def __init__(self, sampler, smoothness=1, lengthscale=1):
self.parameters = ["smoothness", "lengthscale", "n"]
self.smoothness = smoothness
self.lengthscale = lengthscale
self.n = sampler.gen_samples(n=1).shape[-1] + 1
self.compute_eigenpairs()
self.sampler = sampler
self.true_measure = qp.Gaussian(self.sampler)
super(EllipticPDE, self).__init__(dimension_indv=(),dimension_comb=(),parallel=False)
def compute_eigenpairs(self):
_, w, v = get_eigenpairs(self.n)
self.eigenpairs = w, v
def g(self, x):
y = np.zeros(x.shape[:-1])
for i in np.ndindex(y.shape):
y[i] = self.__g(x[i])
return y
def __g(self, x):
w, v = self.eigenpairs
a = np.exp(evaluate(w, v, y=x))
u = pde_solve(a)
return u[len(u)//2]
def _spawn(self, level, sampler):
return EllipticPDE(sampler, smoothness=self.smoothness, lengthscale=self.lengthscale)
# Custom print function
def print_data(data):
for key, val in vars(data).items():
kv = getattr(data, key)
if hasattr(kv, "parameters"):
print(f"{key}: {type(val).__name__}")
for param in kv.parameters:
print(f"\t{param}: {getattr(kv, param)}")
for param in data.parameters:
print(f"{param}: {getattr(data, param)}")
# Main function to test different methods
def test(problem, sampler, stopping_criterion, abs_tol=5e-3, verbose=True, **kwargs):
integrand = problem(sampler)
solution, data = stopping_criterion(integrand, abs_tol=abs_tol, **kwargs).integrate()
if verbose:
print(data)
print("\nComputed solution %.3f in %.2f s"%(solution, data.time_integrate))
Next, let's apply simple Monte Carlo to approximate the expected value $\mathbb{E}[Q]$. The Monte Carlo estimator for $\mathbb{E}[Q]$ is simply the sample average over a finite set of samples, i.e.,
$$\mathcal{Q}_N^\text{MC} := \frac{1}{N} \sum_{n=0}^{N-1} Q^{(n)},$$
where $Q^{(n)} := u(1/2, \boldsymbol{x}^{(n)})$ and we explicitly denote the dependency of $Q$ on the standard normal random numbers $\boldsymbol{x}$ used to sample from the Gaussian random field. We will continue to increase the number of samples $N$ until a certain error criterion is satisfied.
# MC
test(EllipticPDE, qp.IIDStdUniform(32), qp.CubMCCLT)
Data (Data)
solution 0.189
bound_low 0.185
bound_high 0.194
bound_diff 0.010
n_total 22341
time_integrate 6.148
CubMCCLT (AbstractStoppingCriterion)
abs_tol 0.005
rel_tol 0
n_init 2^(10)
n_limit 2^(30)
inflate 1.200
alpha 0.010
EllipticPDE (AbstractIntegrand)
smoothness 1
lengthscale 1
n 33
Gaussian (AbstractTrueMeasure)
mean 0
covariance 1
decomp_type PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
d 2^(5)
replications 1
entropy 183041573331055986930938415720323172091
Computed solution 0.189 in 6.15 s
The solution should be $\approx 0.189$.
Similarly, the quasi-Monte Carlo estimator for $\mathbb{E}[Q]$ is defined as
$$\mathcal{Q}_N^\text{QMC} := \frac{1}{N} \sum_{n=0}^{N-1} Q^{(n)},$$
where $Q^{(n)} := u(1/2, \boldsymbol{t}^{(n)})$ with $\boldsymbol{t}^{(n)}$ the $n$th low-discrepancy point transformed to the distribution of interest. For our elliptic PDE, this means that the quasi-Monte Carlo points, generated inside the unit cube $[0, 1)^d$, are mapped to $\mathbb{R}^d$.
Because the quasi-Monte Carlo estimator doesn't come with a reliable error estimator, we run $K$ different quasi-Monte Carlo estimators in parallel. The sample variance over these $K$ different estimators can then be used as an error estimator.
# QMC
test(EllipticPDE, qp.Lattice(32, replications=8), qp.CubQMCCLT, n_init=32)
Data (Data)
solution 0.185
comb_bound_low 0.182
comb_bound_high 0.188
comb_bound_diff 0.006
comb_flags 1
n_total 1024
n 1024
n_rep 2^(7)
time_integrate 0.276
CubQMCRepStudentT (AbstractStoppingCriterion)
inflate 1
alpha 0.010
abs_tol 0.005
rel_tol 0
n_init 2^(5)
n_limit 2^(30)
EllipticPDE (AbstractIntegrand)
smoothness 1
lengthscale 1
n 33
Gaussian (AbstractTrueMeasure)
mean 0
covariance 1
decomp_type PCA
Lattice (AbstractLDDiscreteDistribution)
d 2^(5)
replications 2^(3)
randomize SHIFT
gen_vec_source kuo.lattice-33002-1024-1048576.9125.txt
order NATURAL
n_limit 2^(20)
entropy 195885308766153261444310305163423607705
Computed solution 0.185 in 0.28 s
3. Multilevel methods¶
Implicit to the Monte Carlo and quasi-Monte Carlo methods above is a discretization parameter used in the numerical solution of the PDE. Let's denote this parameter by $\ell$, $0 \leq \ell \leq L$. Multilevel methods are based on a telescopic sum expansion for the expected value $\mathbb{E}[Q_L]$, as follows:
$$\mathbb{E}[Q_L] = \mathbb{E}[Q_0] + \mathbb{E}[Q_1 - Q_0] + ... + \mathbb{E}[Q_L - Q_{L-1}].$$
Using a Monte Carlo method for each of the terms on the right hand side yields a multilevel Monte Carlo method. Similarly, using a quasi-Monte Carlo method for each term on the right hand side yields a multilevel quasi-Monte Carlo method.
3.1 Multilevel (quasi-)Monte Carlo¶
Our class EllipticPDE needs some changes to be integrated with the multilevel methods in QMCPy.
class MLEllipticPDE(Integrand):
def __init__(self, sampler, smoothness=1, lengthscale=1, _level=None):
self.l = _level
self.parameters = ["smoothness", "lengthscale", "n", "nb_of_levels"]
self.smoothness = smoothness
self.lengthscale = lengthscale
dim = sampler.d + 1
self.nb_of_levels = int(np.log2(dim + 1))
self.n = [2**(l+1) + 1 for l in range(self.nb_of_levels)]
self.compute_eigenpairs()
self.sampler = sampler
self.true_measure = qp.Gaussian(self.sampler)
self.cost = nf
d_out = () if _level is None else (2,)
super(MLEllipticPDE, self).__init__(dimension_indv=d_out,dimension_comb=d_out,parallel=False)
def _spawn(self, level, sampler):
return MLEllipticPDE(sampler, smoothness=self.smoothness, lengthscale=self.lengthscale, _level=level)
def compute_eigenpairs(self):
self.eigenpairs = {}
for l in range(self.nb_of_levels):
_, w, v = get_eigenpairs(self.n[l])
self.eigenpairs[l] = w, v
def g(self, x): # This function is called by keyword reference for the level parameter "l"!
Qf = np.zeros(x.shape[:-1])
for i in np.ndindex(Qf.shape):
Qf[i] = self.__g(x[i], self.l)
if self.l==0:
Qc = np.zeros_like(Qf)
else:
Qc = np.zeros(x.shape[:-1])
for i in np.ndindex(Qf.shape):
Qc[i] = self.__g(x[i], self.l-1)
return np.stack([Qc,Qf],axis=0)
def __g(self, x, l):
w, v = self.eigenpairs[l]
n = self.n[l]
a = np.exp(evaluate(w, v, y=x[:n-1]))
u = pde_solve(a)
return u[len(u)//2]
def _dimension_at_level(self, l):
return self.n[l]
Let's apply multilevel Monte Carlo to the elliptic PDE problem.
test(MLEllipticPDE, qp.IIDStdUniform(32), qp.CubMLMC)
Data (Data)
solution 0.191
n_total 86220
levels 3
n_level [37594 3456 2444]
mean_level [1.906e-01 1.620e-04 4.757e-05]
var_level [0.058 0. 0. ]
cost_per_sample [16. 16. 16.]
alpha 1.768
beta 2.000
gamma 2^(-1)
time_integrate 2.754
CubMLMC (AbstractStoppingCriterion)
rmse_tol 0.002
n_init 2^(8)
levels_min 2^(1)
levels_max 10
theta 2^(-1)
MLEllipticPDE (AbstractIntegrand)
smoothness 1
lengthscale 1
n [ 3 5 9 17 33]
nb_of_levels 5
Gaussian (AbstractTrueMeasure)
mean 0
covariance 1
decomp_type PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
d 2^(5)
replications 1
entropy 307775516609041022608648907501235120384
Computed solution 0.191 in 2.75 s
Now it's easy to switch to multilevel quasi-Monte Carlo. Just change the discrete distribution from IIDStdUniform to Lattice.
test(MLEllipticPDE, qp.Lattice(32,replications=8), qp.CubMLQMCCont, n_init=32)
Data (Data)
solution 0.189
n_total 8704
levels 3
n_level [1024 32 32]
mean_level [ 0.19 -0. -0. ]
var_level [1.545e-06 1.602e-06 1.390e-07]
bias_estimate 5.95e-04
time_integrate 0.526
CubMLQMCCont (AbstractStoppingCriterion)
rmse_tol 0.002
n_init 2^(5)
n_limit 10000000000
replications 2^(3)
levels_min 2^(1)
levels_max 10
n_tols 10
inflate 1.668
theta_init 2^(-1)
theta 0.094
MLEllipticPDE (AbstractIntegrand)
smoothness 1
lengthscale 1
n [ 3 5 9 17 33]
nb_of_levels 5
Gaussian (AbstractTrueMeasure)
mean 0
covariance 1
decomp_type PCA
Lattice (AbstractLDDiscreteDistribution)
d 2^(5)
replications 2^(3)
randomize SHIFT
gen_vec_source kuo.lattice-33002-1024-1048576.9125.txt
order NATURAL
n_limit 2^(20)
entropy 138093628998779790374931977379982164469
Computed solution 0.189 in 0.53 s
3.2 Continuation multilevel (quasi-)Monte Carlo¶
In the continuation multilevel (quasi-)Monte Carlo method, we run the standard multilevel (quasi-)Monte Carlo method for a sequence of larger tolerances to obtain better estimates of the algorithmic parameters. The continuation multilevel heuristic will generally compute the same solution just a bit faster.
test(MLEllipticPDE, qp.IIDStdUniform(32), qp.CubMCMLCont)
Data (Data)
solution 0.190
n_total 29976
levels 3
n_level [16665 1220 256 256 256]
mean_level [1.894e-01 3.968e-04 1.736e-04]
var_level [5.481e-02 3.258e-04 8.380e-06]
cost_per_sample [16. 16. 16.]
alpha 1.192
beta 5.281
gamma 2^(-1)
time_integrate 1.193
CubMLMCCont (AbstractStoppingCriterion)
rmse_tol 0.002
n_init 2^(8)
levels_min 2^(1)
levels_max 10
n_tols 10
inflate 1.668
theta_init 2^(-1)
theta 0.010
MLEllipticPDE (AbstractIntegrand)
smoothness 1
lengthscale 1
n [ 3 5 9 17 33]
nb_of_levels 5
Gaussian (AbstractTrueMeasure)
mean 0
covariance 1
decomp_type PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
d 2^(5)
replications 1
entropy 220267298921554485839372423247101207645
Computed solution 0.190 in 1.19 s
test(MLEllipticPDE, qp.Lattice(32,replications=8), qp.CubMLQMCCont, n_init=32)
Data (Data)
solution 0.190
n_total 8960
levels 3
n_level [1024 64 32]
mean_level [ 0.189 0.001 -0. ]
var_level [1.488e-06 5.978e-07 2.139e-07]
bias_estimate 2.44e-05
time_integrate 0.558
CubMLQMCCont (AbstractStoppingCriterion)
rmse_tol 0.002
n_init 2^(5)
n_limit 10000000000
replications 2^(3)
levels_min 2^(1)
levels_max 10
n_tols 10
inflate 1.668
theta_init 2^(-1)
theta 0.010
MLEllipticPDE (AbstractIntegrand)
smoothness 1
lengthscale 1
n [ 3 5 9 17 33]
nb_of_levels 5
Gaussian (AbstractTrueMeasure)
mean 0
covariance 1
decomp_type PCA
Lattice (AbstractLDDiscreteDistribution)
d 2^(5)
replications 2^(3)
randomize SHIFT
gen_vec_source kuo.lattice-33002-1024-1048576.9125.txt
order NATURAL
n_limit 2^(20)
entropy 308444004426483247734546558914189653748
Computed solution 0.190 in 0.56 s
4. Convergence tests¶
Finally, we will run some convergence tests to see how these methods behave as a function of the error tolerance.
# Main function to test convergence for given problem
def test_convergence(problem, sampler, stopping_criterion, abs_tol=1e-3, verbose=True, smoothness=1, lengthscale=1, **kwargs):
integrand = problem(sampler, smoothness=smoothness, lengthscale=lengthscale)
stopping_crit = stopping_criterion(integrand, abs_tol=abs_tol, **kwargs)
stopping_crit.data = stopping_crit._construct_data()
# manually call "integrate()"
tol = []
n_samp = []
for t in range(stopping_crit.n_tols):
stopping_crit.rmse_tol = stopping_crit.inflate**(stopping_crit.n_tols-t-1)*stopping_crit.target_tol # update tol
stopping_crit._integrate(stopping_crit.data) # call _integrate()
tol.append(copy.copy(stopping_crit.rmse_tol))
n_samp.append(copy.copy(stopping_crit.data.n_level))
if verbose:
print("tol = {:5.3e}, number of samples = {}".format(tol[-1], n_samp[-1]))
return tol, n_samp
# Execute the convergence test
def execute_convergence_test(smoothness=1, lengthscale=1):
# Convergence test for MLMC
tol_mlmc, n_samp_mlmc = test_convergence(MLEllipticPDE, qp.IIDStdUniform(32), qp.CubMLMCCont, verbose=False)
# Convergence test for MLQMC
tol_mlqmc, n_samp_mlqmc = test_convergence(MLEllipticPDE, qp.Lattice(32,replications=8), qp.CubMLQMCCont, verbose=False, n_init=32)
# Compute cost per level
max_levels = max(max([len(n_samp) for n_samp in n_samp_mlmc]), max([len(n_samp) for n_samp in n_samp_mlqmc]))
cost_per_level = np.array([2**level + int(2**(level-1)) for level in range(max_levels)])
cost_per_level = cost_per_level/cost_per_level[-1]
# Compute total cost for each tolerance and store the result
cost = {}
cost["mc"] = (tol_mlmc, [n_samp_mlmc[tol][0] for tol in range(len(tol_mlmc))]) # where we assume V[Q_0] = V[Q_L]
cost["qmc"] = (tol_mlqmc, [n_samp_mlqmc[tol][0] for tol in range(len(tol_mlqmc))]) # where we assume V[Q_0] = V[Q_L]
cost["mlmc"] = (tol_mlmc, [sum([n_samp*cost_per_level[j] for j, n_samp in enumerate(n_samp_mlmc[tol])]) for tol in range(len(tol_mlmc))])
cost["mlqmc"] = (tol_mlqmc, [sum([n_samp*cost_per_level[j] for j, n_samp in enumerate(n_samp_mlqmc[tol])]) for tol in range(len(tol_mlqmc))])
return cost
# Plot the result
def plot_convergence(cost):
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(cost["mc"][0], cost["mc"][1], marker="o", label="MC")
ax.plot(cost["qmc"][0], cost["qmc"][1], marker="o", label="QMC")
ax.plot(cost["mlmc"][0], cost["mlmc"][1], marker="o", label="MLMC")
ax.plot(cost["mlqmc"][0], cost["mlqmc"][1], marker="o", label="MLQMC")
ax.legend(frameon=False)
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel(r"error tolerance $\varepsilon$")
ax.set_ylabel(r"equivalent \# model evaluations at finest level")
plt.show()
This command takes a while to execute (about 1 minute on my laptop):
plot_convergence(execute_convergence_test())
The benefit of the low-discrepancy point set depends on the smoothness of the random field: the smoother the random field, the better. Here's an example for a Gaussian random field with a smaller smoothness $\nu=1/2$ and smaller length scale $\lambda=1/3$.
smoothness = 1/2
lengthscale = 1/3
n = 256
x, w, v = get_eigenpairs(n, smoothness=smoothness, lengthscale=lengthscale)
fig, ax = plt.subplots(figsize=(6, 4))
for _ in range(10):
ax.plot(x, evaluate(w, v))
ax.set_xlim(0, 1)
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$\log(a(x, \cdot))$")
plt.show()
plot_convergence(execute_convergence_test(lengthscale=lengthscale, smoothness=smoothness))
While the multilevel quasi-Monte Carlo method is still the fastest method, the asymptotic cost complexity of the QMC-based methods reduces to approximately the same rate as the MC-based methods.
The benefits of the multilevel methods over single-level methods will be even larger for two- or three-dimensional PDE problems, since it will be even more computationally efficient to take samples on a coarse grid.