#  A replica exchange preconditioned Crank-Nicolson Langevin dynamic MCMC
method for Bayesian inverse problems

Na Ou School of Mathematics and Statistics, Changsha University of Science and
Technology, Changsha, 410114, China<EMAIL_ADDRESS>Zecheng Zhang
Department of Mathematics, Carnegie Mellon University, Pittsburgh, PA 15217,
USA. Guang Lin Department of Mathematics $\&$ School of Mechanical
Engineering, Purdue University, West Lafayete, IN 47907-2067,
<EMAIL_ADDRESS>
ABSTRACT

The preconditioned Crank-Nicolson Markov Chain Monte Carlo (pCN-MCMC) method
is one of the efficient approaches for large-scale Bayesian inference
problems. It is constructed by discretizing the stochastic partial
differential equation (SPDE); when the gradient information is coupled in the
SPDE, the corresponding SPDE is also referred to as the Langevin diffusion. In
this work, we present the replica exchange preconditioned Crank-Nicolson
Langevin dynamic (repCNLD) Monte Carlo method to speed up the convergence of
the pCN-MCMC and tackle the challenge of multimodal distribution simulation.
The replica exchange method comprises two preconditioned Langevin diffusion
chains with low and high temperatures. More specifically, the diffusion chain
with high temperature aims at global exploration and the particles from the
low-temperature diffusion are for the sake of local exploitation. The
particles from the two chains are allowed to exchange with a specific swapping
rate. We provide brief proof that the replica exchange preconditioned Langevin
diffusion (repLD) has the same invariant distribution as the ones without
swapping; the acceleration effect is also confirmed. We use the Crank-Nicolson
scheme to discretize the repLD, and set the preconditioned matrix as the
covariance matrix of the Gaussian prior. The discretization error is analyzed.
Moreover, we use the multi-variance strategy to save the computational cost
for the proposed method. That is, we adopt a high-fidelity method to evaluate
the forward model in the low-temperature chain while employing a low-fidelity
method to compute the forward model in the high-temperature chain. Therefore,
perturbations are induced into the energy function, leading to the biased
swapping rate. To correct the bias caused by the errors in the energy
function, we provide an unbiased estimator derived from some assumptions on
the forward solvers. Furthermore, we provide the adjoint-based repCNLD
algorithm under the discrete frame to speed up the simulations, where the
gradient of the log-likelihood is calculated efficiently by the adjoint
method. The proposed method is tested by sampling Gaussian mixture
distribution and solving PDE-based nonlinear inverse problems numerically.

Keywords: Replica exchange Monte Carlo, Langevin diffusion, preconditioned
Crank-Nicolson scheme, Bayesian inverse problem

##  1 Introduction

Markov chain Monte Carlo (MCMC) approaches are widely used in many fields to
study the statistical properties of complex systems and sample the posterior
distribution in the Bayesian framework [30, 25]. The posterior usually
involves high dimensionality, nonlinearity, and multimodality. The standard
MCMC approaches, such as the random walk Metropolis-Hastings (RWMH), can mix
arbitrarily slowly as the dimension of the unknowns increase [2, 26, 8], which
means the number of forwarding model simulations may tend to infinity.
Moreover, the nonlinearity of the map between the unknowns and output makes it
computationally prohibitive for the Markov chain construction.

A family of dimension-independent MCMC algorithms was presented in [3, 8], of
which the proposal distribution was constructed by the preconditioned Crank-
Nicolson (pCN) discretization of a stochastic partial differential equation
(SPDE) that preserves the reference or target measure. When the proposal is
generalized to the Langevin setting, the Metropolis-adapted Langevin algorithm
(MALA) [31, 29] is resulted, which belongs to Langevin Monte Carlo (LMC)
algorithms. Another branch of the LMC algorithms is the unadjusted Langevin
algorithm (ULA)[31, 12, 9], where the Markov chain is constructed by naturally
discretizing the Langevin diffusion [21, 4]. The ULA is computationally
efficient for sampling high-dimensional probabilities due to the absence of
the Metropolis-Hastings acceptance step [17]. In this paper, we will use the
ULA to build the Markov chains efficiently.

Note that the stochastic noise in the Langevin diffusion involved with the
pCN-MCMC approach is scaled by a temperature parameter and a positive-definite
preconditioned matrix; they greatly affect the efficiency of the ULA method;
we focus on the effect of the temperature parameter in this paper. It plays a
crucial role in global exploration and local exploitation trade. When the
temperature is low, the particle can exploit the local structure of the
likelihood landscape; however, it may trap in the local mode and hence fail to
explore the global domain. On the other hand, when the temperature is high,
the particle can explore the global domain. However, it may not converge to
the local optima.

To exploit the temperature effects and better use the advantages of the high
and low temperatures, the authors in [5] propose the replica exchange Langevin
diffusion (reLD) method. It comprises two Langevin diffusion chains with two
temperatures. By allowing swapping the particle between two chains, the
authors show that the low-temperature chain will converge to the same
stationary distribution with fewer iteration steps.

Authors in [11] later extend the reLD as an optimization algorithm to find the
optima of the neural network. In practice, similar to the gradient descent
neural network optimization algorithms, it is impossible to send all samples
to the GPU; the authors hence also assume a mini-batch setting. Specifically,
one can only obtain an estimation of the gradient with a stochastic batch of
samples. As a result, they call their algorithms replica exchange stochastic
gradient LD, or reSGLD.

It is known that the reLD can speed up the convergence, i.e., one may need
fewer iterations to reach the stationary distribution. However, since the
algorithm employs two chains, the per-iteration cost doubles. Consequently,
the total computation cost may be even higher than the single-chain
algorithms. The multi-variance replica exchange method is developed in[24] to
save computational cost. The method is then applied to handle noise data in
the operator learning optimization [23]. The idea is to use two different
estimations in two chains. Specifically, one can use a low-accuracy estimation
in the high-temperature chain, which explores the region while adopting the
normal estimator in the low-temperature chain. One hence can save costs in the
high-temperature chain without sacrificing the convergence too much.

We propose a replica exchange pCN Langevin dynamic (repCNLD) MCMC method. The
replica exchange method comprises two preconditioned Langevin diffusion chains
(repLD); these Langevin diffusion algorithms are assigned different
temperature parameters. The particles from the two diffusion chains are
allowed to swap with some specific probability. We prove that the repLD with
swapping has the same invariant distribution with the ones without swapping,
the acceleration effect of the repLD method is also confirmed in theory, the
convergence of the diffusion chains is affected by the temperature parameter,
preconditioned matrix, and the swapping intensity.

To discretize the Langevin diffusion, we use the Crank-Nicolson algorithm [8].
It is a linear-implicit scheme [19, 3], the implicitness parameter is
incorporated only in the linear part of the drift to allow for analytical
tractability. The Crank-Nicolson algorithm is the special case for the linear-
implicit scheme introduced in [3]. Since we set the preconditioned matrix to
be the covariance matrix of the Gaussian prior, we refer to it as the
preconditioned Crank-Nicolson scheme following [8]. We analyze the
discretization error; the mean squared error of discretization grows linearly
against the step size.

As the temperature parameter can be treated as an adjustment factor to the
variance of the observed data given the predicted output, we adopt the multi-
variance strategy [24] to save the doubled cost in two chains. The forward
model in the chain with low temperature is solved with the high-fidelity
method, while the chain with high temperature is assigned with a low-fidelity
forward solver. This leads to the swapping rate bias, we correct the bias from
the perspective of forwarding model approximation error. The modified factor
in the swapping rate provides a guide for the choice of the low-fidelity model
used in the second chain.

To further speed up the evaluation of the forward models, we provide the
framework of the discrete adjoint method [1] to calculate the Jacobian of the
output for nonlinear inverse problems. The derivative of the log-likelihood
function is calculated through the forward model and the corresponding costate
equation, which requires much less computational cost. The replica exchange
methods are successfully used to escape from the local optima and explore
distributions with multiple modes. The efficiency improvement and convergence
acceleration effects are confirmed by some numerical examples. We summarize
our contributions as follows.

  1. 1.

We propose a replica exchange preconditioned Langevin diffusion (repLD) method
to handle the multimodal distribution simulation and speed up the sampling of
the posterior density. The preconditioned matrix plays an important part in
acceleration.

  2. 2.

We use the Crank-Nicolson algorithm to discrete the repLD, and set the
preconditioned matrix as the covariance of the Gaussian prior. We also present
the discretization error analysis.

  3. 3.

We propose an unbiased estimator of the swapping rate for the multi-variance
repCNLD method, which provides a guide for the choice of the low-fidelity
model used in the second chain.

  4. 4.

The framework of the discrete adjoint method is provided to efficiently
calculate the gradient for PDE nonlinear inverse problems.

The paper is organized as follows. In Section 2, we describe the setup of the
Bayesian inference and the Crank-Nicolson scheme and present the repLDs and
repCNLD in detail, some results are demonstrated. Section 3 considers the
methods for simulation speeding up: the multi-variance strategy and the
discrete adjoint method. Some numerical experiments are presented to
demonstrate the efficiency of the proposed method in Section 4. Section 5
concludes with a summary discussion and tissues for future research.

##  2 Replica exchange pCN Langevin dynamics MCMC

###  2.1 Bayesian inference

Denote the unknown parameter as $\boldsymbol{\xi}\in{\mathbb{R}}^{n}$, let
${\bf G}:{\boldsymbol{\xi}}\rightarrow{\mathbb{R}}^{n_{d}}$ be the forward
map, where $n_{d}$ is the dimension of the observed data. We assume that the
unknown parameter and the measurement data ${\bf d}$ are related by

${\bf d}={\bf G}({\boldsymbol{\xi}})+\varepsilon,$ (2.1)

where $\varepsilon\in{\mathbb{R}}^{n_{d}}$ is a vector of random noise and
$\varepsilon\sim{\cal N}(0,\sigma_{o}^{2}I)$. Both ${\bf d}$ and
$\boldsymbol{\xi}$ are random variables under the framework of Bayesian
inference, as the Bayes’ rule states, the posterior distribution of
$\boldsymbol{\xi}$ is

$\pi({\boldsymbol{\xi}}|{\bf d})\propto\pi({\bf
d}|\boldsymbol{\xi})\pi_{pri}(\boldsymbol{\xi}).$ (2.2)

The data enters the posterior through the likelihood function $\pi({\bf
d}|\boldsymbol{\xi})$, it is given by

$\pi({\bf
d}|\boldsymbol{\xi})=(2\pi\sigma_{o}^{2})^{-\frac{n_{d}}{2}}\exp\bigg{(}-\frac{\|{\bf
d}-{\bf G}(\boldsymbol{\xi})\|_{2}^{2}}{2\sigma_{o}^{2}}\bigg{)},$

where $\|\cdot\|_{2}$ refers to the Euclidean norm. $\pi({\bf
d}|\boldsymbol{\xi})$ is the likelihood function and
$\pi_{pri}(\boldsymbol{\xi})$ is prior information of the unknown parameter
$\boldsymbol{\xi}$. We focus on the Gaussian prior in this paper, which is
widely used in statistical applications [32, 16]. Without loss of generality,
we assume the prior is ${\cal N}(m,B)$, i.e.,

$\pi_{pri}(\boldsymbol{\xi})\propto\exp\bigg{(}-\frac{(\boldsymbol{\xi}-m)^{T}B^{-1}(\boldsymbol{\xi}-m)}{2}\bigg{)},$

where $m\in{\mathbb{R}}^{n}$ is the mean, and the positive and symmetric
matrix $B\in\mathbb{R}^{n\times n}$ is the covariance matrix of the Gaussian
prior.

Let

$\psi(\boldsymbol{\xi})=-\log\pi({\bf d}|\boldsymbol{\xi})$

be the potential and

$U(\boldsymbol{\xi})=\frac{1}{2}(\boldsymbol{\xi}-m)^{T}B^{-1}(\boldsymbol{\xi}-m)+\psi(\boldsymbol{\xi})$
(2.3)

be the energy function, the posterior in the equation $(\ref{posterior})$ can
then be expressed in the form

$\pi({\boldsymbol{\xi}}|{\bf d})\propto\exp(-U(\boldsymbol{\xi})).$

From here on, we omit the conditioning on data ${\bf d}$ and write the
posterior density as $\pi({\boldsymbol{\xi}})$.

For some cases, when the forward map ${\bf G}$ is linear, the adoption of
Gaussian priors lead to Gaussian posteriors [15], the statistical estimators,
e.g., the mean, the maximum a posteriori (MAP), covariance, credible intervals
can be obtained easily. However, the forward map ${\bf G}$ is often nonlinear,
resulting in the lack of an analytical expression of the posterior
distribution, and poses a significant challenge to estimating the quantity of
interests.

###  2.2 The preconditioned Crank-Nicolson algorithms

MCMC methods are often used to sample from the posterior distribution that
arises from Bayesian inference. The preconditioned Crank-Nicolson MCMC (pCN-
MCMC) method is one of the efficient approaches to dealing with large-scale
inference problems. It is constructed by a discretization of stochastic
partial differential equations (SPDEs)

${\rm
d}\boldsymbol{\xi}_{t}=-L\left(B^{-1}(\boldsymbol{\xi}_{t}-m)+\alpha{\nabla\psi(\boldsymbol{\xi}_{t})}\right){\rm
d}t+\sqrt{2L\tau}{\rm d}\boldsymbol{W}_{t},$ (2.4)

where $\nabla\psi$ is the negative derivative of the log-likelihood function,
$L$ is an arbitrary symmetric, positive-definite matrix,
${\boldsymbol{W}}_{t}\in\mathbb{R}^{n}$ is a standard $n$-dimensional Brown
motion and $\tau>0$ is the temperature parameter. The SPDE preserves the
reference measure ${\cal N}(m,\tau B)$ for $\alpha=0$ or the target measure

$\pi({\boldsymbol{\xi}})\propto\exp\left(-\frac{U(\boldsymbol{\xi})}{\tau}\right)$

for $\alpha=1$ [33, 3, 20].

Our goal is to obtain the target measure through the natural discretization of
the Langevin diffusion or SPDE. Hence we set $\alpha=1$ in the SPDE
$(\ref{cn1})$.  We use the liner-implicit scheme [3] to obtain the discrete-
time version of the continuous process, yielding,

${\boldsymbol{\xi}}_{k+1}-{\boldsymbol{\xi}}_{k}=-\left((1-\theta)LB^{-1}{\boldsymbol{\xi}}_{k}+\theta
LB^{-1}{\boldsymbol{\xi}}_{k+1}-LB^{-1}m+L\nabla\psi({\boldsymbol{\xi}}_{k})\right)\delta+\sqrt{2L\delta\tau}{\boldsymbol{w}}_{k},$

for some $\theta\in(0,1)$, ${\boldsymbol{w}}_{k}\sim{\cal N}(0,I)$,
$I\in\mathbb{R}^{n\times n}$ is the identity matrix and $\delta$ is the time
step. The formula can be rearranged as

$\left(I+\theta\delta
LB^{-1}\right){\boldsymbol{\xi}}_{k+1}=\left(I-(1-\theta)\delta
LB^{-1}\right){\boldsymbol{\xi}}_{k}+LB^{-1}\delta\left(m-B{\nabla\psi({\boldsymbol{\xi}}_{k})}\right)+\sqrt{2L\delta\tau}{\boldsymbol{w}}_{k},$

since the matrix $(I+\theta\delta LB^{-1})$ is symmetric positive-definite,
the uniqueness of the location ${\boldsymbol{\xi}}_{k+1}$ is guaranteed given
the realization of ${\boldsymbol{\xi}}_{k}$. Convergence of such implicit
methods for SPDEs is investigated in [28].

When the implicit parameter is set as $\theta=1/2$, we have the Crank-Nicolson
scheme,

$\left(I+\frac{1}{2}\delta
LB^{-1}\right){\boldsymbol{\xi}}_{k+1}=\left(I-\frac{1}{2}\delta
LB^{-1}\right){\boldsymbol{\xi}}_{k}+LB^{-1}\delta\left(m-B{\nabla\psi({\boldsymbol{\xi}}_{k})}\right)+\sqrt{2L\delta\tau}{\boldsymbol{w}}_{k}.$
(2.5)

Moreover, we set the preconditioned matrix $L$ to be the covariance matrix $B$
of the Gaussian prior, the equation $(\ref{discrete})$ becomes

$\left(1+\frac{1}{2}\delta\right){\boldsymbol{\xi}}_{k+1}=\left(1-\frac{1}{2}\delta\right){\boldsymbol{\xi}}_{k}+\delta\left(m-B\nabla\psi({\boldsymbol{\xi}}_{k})\right)+\sqrt{2B\delta\tau}{\boldsymbol{w}}_{k}.$
(2.6)

Let

$\beta=\frac{2\sqrt{2\delta}}{2+\delta},$

then for $\delta\in(0,2)$, we have $\beta\in[0,1]$. The formula of the
preconditioned Crank-Nicolson Langevin dynamics (pCNLD) is resulted,

${\boldsymbol{\xi}}_{k+1}=\sqrt{1-\beta^{2}}{\boldsymbol{\xi}}_{k}+(1-\sqrt{1-\beta^{2}})(m-B\nabla\psi({\boldsymbol{\xi}}_{k}))+\beta\sqrt{B\tau}{\boldsymbol{w}}_{k}.$
(2.7)

Or equivalently, by replacing the term
$B^{-1}({\boldsymbol{\xi}}_{k}-m)+\nabla\psi({\boldsymbol{\xi}}_{k})$ with
$\nabla U({\boldsymbol{\xi}}_{k})$ in the equation $(\ref{pCN-delta})$, it can
be expressed as

${\boldsymbol{\xi}}_{k+1}={\boldsymbol{\xi}}_{k}-(1-\sqrt{1-\beta^{2}})B\nabla
U\left({\boldsymbol{\xi}}_{k}\right)+\beta\sqrt{B\tau}\boldsymbol{w}_{k}.$
(2.8)

The pCNLD method belongs to the ULA, the MCMC samples generated by this method
produce a biased approximation of the target distribution [31, 13], where the
bias is stemmed from the discretization error. The asymptotic bias can be
reduced when the step size is sufficiently small.

######  Remark 2.1.

For the Langevin dynamics discussed above, the Gaussian prior is incorporated
into the posterior density; when the hybrid prior information is imposed,
e.g., the Laplace and the Gaussian prior are jointly suggested, it is
equivalent to set

$\psi(\boldsymbol{\xi})=-\log\pi(\bf
d|\boldsymbol{\xi})+\gamma|\boldsymbol{\xi}|$

where $\gamma$ is the hyperparameter in the central Laplace distribution and
$|\cdot|$ is the $l_{1}$ norm of the vector.

###  2.3 Replica exchange preconditioned Langevin diffusion

We note that the temperature parameter $\tau$ and preconditioned matrix $L$
play crucial roles in the Langevin diffusion, we focus on the effect of the
temperature parameter in this paper. When the temperature is low, the particle
$\boldsymbol{\xi}_{t}$ can exploit the local structure and converge locally
fast, however, it may get trapped in the local optima and does not travel
globally. As a result, when the solution landscape is multi-mode, the particle
may not converge to the global optima in practice. On the other hand, when the
temperature parameter $\tau$ is large, the convergence of the diffusion can be
accelerated since the particles are able to explore the entire state space.

Inspired by the success of replica exchange Langevin diffusion presented in
paper [5], we use two chains and allow swapping to tackle the multimodal
problems and improve the pCNLD convergence efficiency. Consider a pair of
particles governed by two preconditioned Langevin diffusion chains as defined
in $(\ref{cn1})$ with $\alpha=1$ and temperatures $\tau_{1}<\tau_{2}$,

$\begin{split}\mathrm{d}\boldsymbol{\xi}^{(1)}_{t}=-L\left(B^{-1}(\boldsymbol{\xi}^{(1)}_{t}-m)+\nabla\psi(\boldsymbol{\xi}^{(1)}_{t})\right){\rm
d}t+\sqrt{2L\tau_{1}}\rm{d}\boldsymbol{W}^{(1)}_{t},\\\
\mathrm{d}\boldsymbol{\xi}^{(2)}_{t}=-L\left(B^{-1}(\boldsymbol{\xi}^{(2)}_{t}-m)+\nabla\psi(\boldsymbol{\xi}^{(2)}_{t})\right){\rm
d}t+\sqrt{2L\tau_{2}}\rm{d}\boldsymbol{W}^{(2)}_{t},\end{split}$

they can also be expressed equivalently as,

$\begin{split}\mathrm{d}\boldsymbol{\xi}^{(1)}_{t}=-L\nabla
U(\boldsymbol{\xi}^{(1)}_{t}){\rm d}t+\sqrt{2L\tau_{1}}{\rm
d}\boldsymbol{W}^{(1)}_{t},\\\ \mathrm{d}\boldsymbol{\xi}^{(2)}_{t}=-L\nabla
U(\boldsymbol{\xi}^{(2)}_{t}){\rm d}t+\sqrt{2L\tau_{2}}{\rm
d}\boldsymbol{W}^{(2)}_{t},\end{split}$ (2.9)

where $\boldsymbol{W}^{(1)}$ is independent with $\boldsymbol{W}^{(2)}$. The
equation $(\ref{cns-U})$ converges to the invariant distribution

$\mu(\boldsymbol{\xi}^{(1)},\boldsymbol{\xi}^{(2)})=\pi(\boldsymbol{\xi}^{(1)},\boldsymbol{\xi}^{(2)})\rm{d}\boldsymbol{\xi}^{(1)}\rm{d}\boldsymbol{\xi}^{(2)},$

with the density

$\pi(\boldsymbol{\xi}^{(1)},\boldsymbol{\xi}^{(2)})\propto
e^{-\frac{U(\boldsymbol{\xi}^{(1)})}{\tau_{1}}-\frac{U(\boldsymbol{\xi}^{(2)})}{\tau_{2}}}.$
(2.10)

When the two particles are allowed to swap from position
$(\boldsymbol{\xi}_{t}^{(1)},\boldsymbol{\xi}_{t}^{(2)})$ at time $t$ to
position $(\boldsymbol{\xi}_{t+dt}^{(2)},\boldsymbol{\xi}_{t+dt}^{(1)})$ with
a swapping rate $cs(\boldsymbol{\xi}_{t}^{(1)},\boldsymbol{\xi}_{t}^{(2)})dt$,
where $c\geq 0$ is the swapping intensity, and

$s(\boldsymbol{\xi}_{t}^{(1)},\boldsymbol{\xi}_{t}^{(2)})=\min\bigg{\\{}1,\exp\left(\tau_{\delta}\left(U(\boldsymbol{\xi}_{t}^{(1)})-U(\boldsymbol{\xi}_{t}^{(2)})\right)\right)\bigg{\\}},$
(2.11)

where $\tau_{\delta}=\frac{1}{\tau_{1}}-\frac{1}{\tau_{2}}$. We have the
following result:

######  Theorem 2.1 (Reversibility and invariance).

For any $c$ with $cdt<1$ satisfied, the infinitesimal generator of the
Langevin jump process $\\{\boldsymbol{\xi}_{t}\\}_{t\geq 0}$ defined by
equations $(\ref{cns-U})$ and $(\ref{swap})$ is

$\begin{array}[]{l}\mathcal{L}^{c}\left(f(x,y)\right)=\underbrace{-[L\nabla_{x}U(x)]\cdot\nabla_{x}f(x,y)+\tau_{1}L:\nabla_{x}\nabla_{x}f(x,y)}_{\mathcal{L}^{c}_{1}\left(f(x,y)\right)}\\\
\underbrace{-[L\nabla_{y}U(y)]\cdot\nabla_{y}f(x,y)+\tau_{2}L:\nabla_{y}\nabla_{y}f(x,y)}_{\mathcal{L}^{c}_{2}\left(f(x,y)\right)}+\underbrace{cs(x,y)\left(f(y,x)-f(x,y)\right)}_{\mathcal{L}^{c}_{s}\left(f(x,y)\right)},\end{array}$

the domain of the generator $\mathcal{D}(\mathcal{L}^{c})$ is the space of all
twice-differentiable functions with compact support. The symbol $\cdot$ and
$:$ represent the Frobenius inner products, i.e., for any
$x,y\in\mathbb{R}^{n}$,

$\displaystyle[L\nabla_{x}U(x)]\cdot\nabla_{x}f(x,y)$ $\displaystyle=$
$\displaystyle\sum_{i}\left[L\nabla_{x}U(x)\right]_{i}\left(\nabla_{x}f(x,y)\right)_{i}$
$\displaystyle L:\nabla_{x}\nabla_{x}f(x,y)$ $\displaystyle=$
$\displaystyle\sum_{i,j}\left[L\right]_{i,j}\frac{\partial^{2}f(x,y)}{\partial
x_{i}\partial x_{j}}.$

Furthermore, $\\{\boldsymbol{\xi}_{t}\\}_{t\geq 0}$ is reversible and its
invariant distribution has density $(\ref{invariant})$.

The proof of the Theorem 2.1 is analog to the proof of the Lemma 3.2 in paper
[5], the term $\mathcal{L}_{s}^{c}$ arises from the swap has the same form as
the one in [5], the symmetry of the preconditioned matrix $L$ makes the
equality formula involved with $\mathcal{L}_{1}^{c}$ and $\mathcal{L}_{2}^{c}$
established, a detailed proof is provided in the Appendix.

In the following,  we will prove that in addition to the temperature
parameters and the swapping, the preconditioned matrix also contributes to the
acceleration of the replica exchange preconditioned Langevin diffusion chains
(repLD). Let $\nu_{t}$ be the distribution of the repLD
$\\{\boldsymbol{\xi}_{t}\\}_{t\geq 0}$ at time $t$, $\mu$ be the invariant
distribution. According to the analysis in [5], the derivative of
$\chi^{2}(\nu_{t}\|\mu)$ with time $t$ can be characterized by the Dirichlet
form, i.e.,

$\frac{\mathrm{d}}{\mathrm{d}t}\chi^{2}\left(\nu_{t}\|\mu\right)=-2\mathcal{E}^{c}\left(\frac{\mathrm{d}\nu_{t}}{\mathrm{d}\mu}\right),$
(2.12)

where the $\chi^{2}$-divergence is used to quantify the discrepancy between
$\nu_{t}$ and $\mu$, it is defined as

$\chi^{2}\left(\nu_{t}\|\mu\right)=\int\left(\frac{\mathrm{d}\nu_{t}}{\mathrm{d}\mu}-1\right)^{2}\mathrm{d}\mu.$

The Dirichlet form is defined as

$\mathcal{E}^{c}(f)=\int\Gamma^{c}(f)\mathrm{d}\mu,$

where

$\Gamma^{c}(f)=\frac{1}{2}\left(\mathcal{L}^{c}\left(f^{2}\right)-2f\mathcal{L}^{c}(f)\right)$

is the Carré du Champ operator, $\mu$ and $\mathcal{L}^{c}$ is the invariant
distribution and infinitesimal generator of $\\{\boldsymbol{\xi}_{t}\\}_{t\geq
0}$, respectively.

######  Theorem 2.2 (Acceleration).

For any fixed function $f$, $\mathcal{E}^{c}(f)$ has the form

$\displaystyle{\mathcal{E}}^{c}(f)$
$\displaystyle=\underbrace{\int\left(\tau_{1}\left\|\sqrt{L}\nabla_{x}f(x,y)\right\|^{2}+\tau_{2}\left\|\sqrt{L}\nabla_{y}f(x,y)\right\|^{2}\right)\mathrm{d}\mu\left(x,y\right)}_{{\mathcal{E}}^{0}\geq
0}$ (2.13)
$\displaystyle+\underbrace{\int\frac{c}{2}s\left(x,y\right)\left(f(y,x)-f(x,y)\right)^{2}\mathrm{d}\mu\left(x,y\right)}_{\geq
0},$

The derivation of $\mathcal{E}^{c}(f)$ is provided in the Appendix. The first
term $\mathcal{E}^{0}$ in $(\ref{Dirichlet})$ corresponds to the repLD without
swapping, larger temperatures $\tau_{1}$, $\tau_{2}$ and appropriate matrix
$L$ may result in a larger $\mathcal{E}^{0}$, which means a faster rate of
convergence. The nonnegativity of the second term in $(\ref{Dirichlet})$ can
further promotes the Dirichlet form as analysed in [5] and therefore leads to
faster evolution of $\chi^{2}\left(\nu_{t}\|\mu\right)$.

###  2.4 Replica exchange pCN Langevin dynamics

We discretize the repLD with the pCN scheme introduced in the equation
$(\ref{pCN-langevin})$ or $(\ref{pCN-U})$, the pair of particles are then
updated following

$\begin{split}{{\boldsymbol{\xi}}}^{(1)}_{k+1}=\sqrt{1-\beta^{2}}{{\boldsymbol{\xi}}}^{(1)}_{k}+\left(1-\sqrt{1-\beta^{2}}\right)\left(m-B\nabla\psi({{\boldsymbol{\xi}}}^{(1)}_{k})\right)+\beta\sqrt{B\tau_{1}}\boldsymbol{w}^{(1)}_{k},\\\
{{\boldsymbol{\xi}}}^{(2)}_{k+1}=\sqrt{1-\beta^{2}}{{\boldsymbol{\xi}}}^{(2)}_{k}+\left(1-\sqrt{1-\beta^{2}}\right)\left(m-B\nabla\psi({{\boldsymbol{\xi}}}^{(2)}_{k})\right)+\beta\sqrt{B\tau_{2}}\boldsymbol{w}^{(2)}_{k},\end{split}$
(2.14)

or

$\begin{split}{\boldsymbol{\xi}}^{(1)}_{k+1}={\boldsymbol{\xi}}^{(1)}_{k}-\left(1-\sqrt{1-\beta^{2}}\right)B\nabla
U\left({\boldsymbol{\xi}}^{(1)}_{k}\right)+\beta\sqrt{B\tau_{1}}\boldsymbol{w}^{(1)}_{k},\\\
{{\boldsymbol{\xi}}}^{(2)}_{k+1}={{\boldsymbol{\xi}}}^{(2)}_{k}-\left(1-\sqrt{1-\beta^{2}}\right)B\nabla
U\left({\boldsymbol{\xi}}^{(1)}_{k}\right)+\beta\sqrt{B\tau_{2}}\boldsymbol{w}^{(2)}_{k},\end{split}$
(2.15)

where the tunable parameter $\beta$ depends on the time step $\delta$,
${\boldsymbol{w}}^{(i)}_{k}\sim{\cal N}(0,I)$, for $i=1,2$. We call the
schemes above the replica exchange preconditioned Crank-Nicolson Langevin
dynamics (repCNLD) method. In order to analyse the discretization error of the
repCNLD method, we relax the terms $\delta$ and $\sqrt{\delta}$ from the
$\beta$ involved parameter in the equation $(\ref{repCNLD-U})$, and obtain the
equivalent form

$\begin{split}{\boldsymbol{\xi}}^{(1)}_{k+1}={\boldsymbol{\xi}}^{(1)}_{k}-\eta\delta
B\nabla
U\left({\boldsymbol{\xi}}^{(1)}_{k}\right)+\eta\sqrt{2B\tau_{1}\delta}\boldsymbol{w}^{(1)}_{k},\\\
{{\boldsymbol{\xi}}}^{(2)}_{k+1}={{\boldsymbol{\xi}}}^{(2)}_{k}-\eta\delta
B\nabla
U\left({\boldsymbol{\xi}}^{(1)}_{k}\right)+\eta\sqrt{2B\tau_{2}\delta}\boldsymbol{w}^{(2)}_{k},\end{split}$
(2.16)

where $\eta=\frac{2}{2+\delta}$.

We have some assumptions for the energy function. Firstly, we assume that the
energy function is $K$-smooth, that is, for all $x,y\in{\mathbb{R}}^{d}$,
there exists a constant $K<+\infty$ such that

$\|\nabla U(x)-\nabla U(y)\|\leq K\|x-y\|,$ (2.17)

we also assume that there exist positive constants $\alpha_{1}$ and
$\alpha_{2}$, such that for all $x\in{\mathbb{R}}^{d}$,

$[B\nabla U(x)]\cdot x\geq\alpha_{1}\|x\|^{2}-\alpha_{2}.$ (2.18)

The repLD can be characterized as

$\mathrm{d}\boldsymbol{\xi}_{t}=-H\nabla
D(\boldsymbol{\xi}_{t})\mathrm{d}t+\sqrt{H}\Sigma_{t}\mathrm{d}\boldsymbol{W}_{t},$
(2.19)

where
$\boldsymbol{\xi}_{t}=\left[\boldsymbol{\xi}_{t}^{(1)};\boldsymbol{\xi}_{t}^{(2)}\right]\in{\mathbb{R}}^{2n}$,
$D(\boldsymbol{\xi}_{t})=\left[U(\boldsymbol{\xi}_{t}^{(1)});U(\boldsymbol{\xi}_{t}^{(2)})\right]$,
${\boldsymbol{W}}_{t}\in{\mathbb{R}}^{2n}$ is the Brown motion, the matrix
$H=\left(\begin{array}[]{cc}B&0\\\ 0&B\end{array}\right)$,
$\sqrt{H}=\left(\begin{array}[]{cc}\sqrt{B}&0\\\
0&\sqrt{B}\end{array}\right)$, here $B$ is the covariance matrix of the
Gaussian prior, since the preconditioned matrix is set as $B$. $\Sigma_{t}$ is
a random matrix which switches between
$\left(\begin{array}[]{cc}\sqrt{2\tau_{1}I}&0\\\
0&\sqrt{2\tau_{2}I}\end{array}\right)$ and
$\left(\begin{array}[]{cc}\sqrt{2\tau_{2}I}&0\\\
0&\sqrt{2\tau_{1}I}\end{array}\right)$ with the probability
$cs\left(\boldsymbol{\xi}_{t}^{(1)},\boldsymbol{\xi}_{t}^{(2)}\right)\mathrm{d}t$.

Denote
$\boldsymbol{\xi}^{\delta}(k)=\left[\boldsymbol{\xi}^{\delta,1}(k);\boldsymbol{\xi}^{\delta,2}(k)\right]$
as the solution of the discretization scheme $(\ref{repCNLDs})$ or
$(\ref{repCNLD-U})$ at the $k$-th iteration, i.e.,
$\boldsymbol{\xi}^{\delta}(k)$ satisfies

$\boldsymbol{\xi}^{\delta}(k+1)={\boldsymbol{\xi}}^{\delta}(k)-\eta\delta
H\nabla
D\left({\boldsymbol{\xi}}^{\delta}(k)\right)+\eta\sqrt{H\delta}\Sigma^{\delta}(k)\boldsymbol{w}(k),$

where $\boldsymbol{w}(k)\in{\mathbb{R}}^{2d}$ is a random vector generated
from the standard Gaussian distribution. $\Sigma^{\delta}(k)$ is defined the
same as $\Sigma_{t}$ but switches with the probability
$cs\left(\boldsymbol{\xi}^{\delta,1}(k),\boldsymbol{\xi}^{\delta,2}(k)\right)\delta$.
Let $\\{\boldsymbol{\xi}^{\delta}_{t}\\}_{t\geq 0}$ be the continuous-time
interpolation of $\\{\boldsymbol{\xi}^{\delta}(k)\\}_{k\geq 0}$ which is
defined as

$\boldsymbol{\xi}_{t}^{\delta}=\boldsymbol{\xi}_{0}-\eta\int_{0}^{t}H\nabla
D\left({\boldsymbol{\xi}}^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}\right){\rm
d}s+\eta\int_{0}^{t}\sqrt{H}\Sigma^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}{\rm d}\boldsymbol{W}_{s},$ (2.20)

here $\lfloor\cdot\rfloor$ denotes the floor function. It can be verified that
given $\boldsymbol{\xi}_{0}={\boldsymbol{\xi}}^{\delta}(0)$, we then have

$\boldsymbol{\xi}_{t}^{\delta}=\boldsymbol{\xi}^{\delta}(k),$

for all integers $k>0$ and $t=k\delta$.

######  Theorem 2.3 (Discretization Error).

When the energy function satisfies assumptions $(\ref{smooth})$ and
$(\ref{B-dissipative})$, there exists a constant
$\epsilon(n,\tau_{1},\tau_{2},c,K,\alpha_{1},\alpha_{2},T,B)$ that depends on
the dimension $n$, temperature parameters $\tau_{1}$ and $\tau_{2}$, swapping
intensity $c$, smoothness constant $K$ and constants $(\alpha_{1},\alpha_{2})$
in the inequality $(\ref{B-dissipative})$, length of the time interval $T$,
and the covariance matrix of the Gaussian prior $B$, such that for all
$t\in[0,T]$,

${\mathbb{E}}\left[\|\boldsymbol{\xi}_{t}-\boldsymbol{\xi}_{t}^{\delta}\|^{2}\right]\leq\epsilon(n,\tau_{1},\tau_{2},c,K,\alpha_{1},\alpha_{2},T,B)\delta.$

provided the step size $\delta\in(0,2)$ satisfies
$0<\frac{2\delta}{2+\delta}<\frac{\alpha_{1}}{2\|B\|_{F}^{2}K^{2}}$.

The proof of Theorem 2.3 adapts from the proof framework of Theorem 3.6 in
paper [5], more details about the proof is demonstrated in the Appendix.

##  3 Simulation speeding up

Compared with the single-chain method, the simulation of the forwarding model
for the replica exchange method doubles. Meanwhile, the gradient of the log-
likelihood requires much computational resources. We provide two schemes to
speed up the simulation of the proposed method: the multi-variance strategy
and the discrete adjoint method.

###  3.1 Multi-variance repCNLD

In practice, one only has an estimation of the energy function $U(\cdot)$ or
its gradient. In the optimization literature, this is due to the mini-batch
setting. Specifically, the limited computing memory only allows one to compute
the gradient by a random (stochastic) batch of samples. In the inverse PDE
literature, errors caused by the forward solver and errors in the observation
all contribute to the inaccurate estimation of the energy function [14, 7].
Consequently, investigators are interested in developing Langevin dynamic
algorithms which assume a perturbation in the gradient estimation. Following
the convention in the optimization society, researchers usually call such
kinds of algorithms stochastic gradient Langevin dynamic (SGLD). The authors
of paper [24] proposed the multi-variance replica exchange stochastic gradient
Langevin dynamics (m-reSGLD) method to accelerate the sampling process.

Similarly, for the proposed repCNLD method, we use the high-fidelity method to
solve the forward model for the chain with low temperature, while in the high-
temperature chain, the forward model is solved with the low-fidelity method.
We refer to it as m-repCNLD method for the convenience of description. As
discussed in paper [11, 24], the biased swapping rate could be modified when
perturbations are induced into the energy function. Denote the energy function
involved term in calculating the swapping rate as

$S(\boldsymbol{\xi}^{(1)},\boldsymbol{\xi}^{(2)})=\exp\left(\tau_{\delta}\left(U(\boldsymbol{\xi}^{(1)})-U(\boldsymbol{\xi}^{(2)})\right)\right),$

it is then approximated by

$\displaystyle\widetilde{S}({{\boldsymbol{\xi}}}^{(1)},{{\boldsymbol{\xi}}}^{(2)})$
$\displaystyle=e^{\tau_{\delta}\left(U({{\boldsymbol{\xi}}}^{(1)})-\widetilde{U}({{\boldsymbol{\xi}}}^{(2)})\right)}$
$\displaystyle=e^{\tau_{\delta}\left(U({{\boldsymbol{\xi}}}^{(1)})-U({{\boldsymbol{\xi}}}^{(2)})\right)}e^{\tau_{\delta}\left(U({{\boldsymbol{\xi}}}^{(2)})-\widetilde{U}({{\boldsymbol{\xi}}}^{(2)})\right)},$
$\displaystyle=S({{\boldsymbol{\xi}}}^{(1)},{{\boldsymbol{\xi}}}^{(2)})e^{\tau_{\delta}\left(U({{\boldsymbol{\xi}}}^{(2)})-\widetilde{U}({{\boldsymbol{\xi}}}^{(2)})\right)},$

here $U(\cdot)$ and $\widetilde{U}(\cdot)$ are calculated by the accurate and
approximated forward model, respectively. The expectation of
$\widetilde{S}({{\boldsymbol{\xi}}}^{(1)},{{\boldsymbol{\xi}}}^{(2)})$ is

$\mathbb{E}\left[\widetilde{S}({{\boldsymbol{\xi}}}^{(1)},{{\boldsymbol{\xi}}}^{(2)})\right]=S({{\boldsymbol{\xi}}}^{(1)},{{\boldsymbol{\xi}}}^{(2)})\mathbb{E}\left[e^{\tau_{\delta}\left(U({{\boldsymbol{\xi}}}^{(2)})-\widetilde{U}({{\boldsymbol{\xi}}}^{(2)})\right)}\right].$
(3.21)

Let ${\bf\widetilde{G}}:{\boldsymbol{\xi}}\rightarrow{\mathbb{R}}^{n_{d}}$ be
the approximated forward output for the second chain, ${\bf G}$ be the
accurate forward output for the first chain, we assume the normality holds for
the accurate and approximated model, i.e.,

${\bf\widetilde{G}}({\boldsymbol{\xi}})\sim{\cal N}({\bf
G}({\boldsymbol{\xi}}),\widetilde{\sigma}^{2}I),\quad\text{for}\
\forall{\boldsymbol{\xi}}.$

The relationship between the observed data and the forward model output can be
induced by equation $(\ref{likelihood})$,

${\bf d}\sim{\cal N}({\bf
G}({\boldsymbol{\xi}}),\sigma_{o}^{2}I),\quad\text{for}\
\forall{\boldsymbol{\xi}}.$

We then have the difference between the accurate and approximated energy.

$\displaystyle U({{\boldsymbol{\xi}}})-{\widetilde{U}}({{\boldsymbol{\xi}}})$
$\displaystyle=\psi({{\boldsymbol{\xi}}})-{\widetilde{\psi}}({{\boldsymbol{\xi}}})$
$\displaystyle=\frac{1}{2\sigma_{o}^{2}}\left(\|{\bf d}-{\bf
G}({\boldsymbol{\xi}})\|^{2}-\|{\bf
d}-{\bf\widetilde{G}}({\boldsymbol{\xi}})\|^{2}\right)$
$\displaystyle=\frac{1}{2\sigma_{o}^{2}}\left(\|{\bf
G}({\boldsymbol{\xi}})-{\bf\widetilde{G}}({\boldsymbol{\xi}})\|^{2}-2\left({\bf
d}-{\bf G}({\boldsymbol{\xi}})\right)^{T}\left({\bf
G}({\boldsymbol{\xi}})-\bf{\widetilde{G}}({\boldsymbol{\xi}})\right)\right).$

The expectation in the right hand side of the equation $(\ref{bias})$ can be
calculated as

$\mathbb{E}\left[e^{\tau_{\delta}\left(U({{\boldsymbol{\xi}}}^{(2)})-\widetilde{U}({{\boldsymbol{\xi}}}^{(2)})\right)}\right]=\mathbb{E}\left[e^{\frac{\tau_{\delta}}{2\sigma_{o}^{2}}\left(Z^{T}Z-2X^{T}Z\right)}\right],$

where

$Z={\bf
G}({{\boldsymbol{\xi}}^{(2)}})-{\bf\widetilde{G}}({{\boldsymbol{\xi}}^{(2)}}),\quad
X={\bf d}-{\bf G}({{\boldsymbol{\xi}}^{(2)}}),$

they are random variables independent with each other for any fixed
${\boldsymbol{\xi}}^{(2)}$, and $Z\sim{\cal N}(0,\widetilde{\sigma}^{2}I)$,
$X\sim{\cal N}(0,\sigma_{o}^{2}I)$.

######  Theorem 3.1.

Let $r_{\sigma}=\frac{\widetilde{\sigma}^{2}}{\sigma_{o}^{2}}$, then when
$r_{\sigma}<\frac{1}{\tau_{\delta}^{2}+\tau_{\delta}}$ satisfied,

$\widetilde{S}_{m}({{\boldsymbol{\xi}}}^{(1)},{{\boldsymbol{\xi}}}^{(2)})=\left[1-(\tau_{\delta}+\tau_{\delta}^{2})r_{\sigma}\right]^{\frac{n_{d}}{2}}e^{\tau_{\delta}(U({{\boldsymbol{\xi}}}^{(1)})-\widetilde{U}({{\boldsymbol{\xi}}}^{(2)}))},$
(3.22)

is an unbiased estimator of
$S(\boldsymbol{\xi}^{(1)},\boldsymbol{\xi}^{(2)})$.

The Theorem 3.1 can be easily proved, according to the definition of the
expectation, we have

$\displaystyle\mathbb{E}\left[e^{\frac{\tau_{\delta}}{2\sigma_{o}^{2}}\left(Z^{T}Z-2X^{T}Z\right)}\right]$
$\displaystyle=\int\int(2\pi\widetilde{\sigma}^{2})^{-\frac{n_{d}}{2}}e^{-\frac{\sigma_{o}^{2}-\widetilde{\sigma}^{2}(\tau_{\delta}^{2}+\tau_{\delta})}{2\sigma_{o}^{2}\widetilde{\sigma}^{2}}z^{T}z}(2\pi\sigma_{o}^{2})^{-\frac{n_{d}}{2}}e^{-\frac{(x+\tau_{\delta}z)^{T}(x+\tau_{\delta}z)}{2\sigma_{o}^{2}}}{\rm
d}x{\rm d}z$
$\displaystyle=\int(2\pi\widetilde{\sigma}^{2})^{-\frac{n_{d}}{2}}e^{-\frac{\sigma_{o}^{2}-\widetilde{\sigma}^{2}(\tau_{\delta}^{2}+\tau_{\delta})}{2\sigma_{o}^{2}\widetilde{\sigma}^{2}}z^{T}z}{\rm
d}z$
$\displaystyle=\left[1-(\tau_{\delta}+\tau_{\delta}^{2})r_{\sigma}\right]^{-\frac{n_{d}}{2}},$

which indicates that the term
$\widetilde{S}({{\boldsymbol{\xi}}}^{(1)},{{\boldsymbol{\xi}}}^{(2)})$ is a
biased estimator of
$S({{\boldsymbol{\xi}}}^{(1)},{{\boldsymbol{\xi}}}^{(2)})$. It can be verified
that

$\displaystyle\mathbb{E}\left[\widetilde{S}_{m}({{\boldsymbol{\xi}}}^{(1)},{{\boldsymbol{\xi}}}^{(2)})\right]$
$\displaystyle=\left[1-(\tau_{\delta}+\tau_{\delta}^{2})r_{\sigma}\right]^{\frac{n_{d}}{2}}\mathbb{E}\left[e^{\tau_{\delta}(U({{\boldsymbol{\xi}}}^{(1)})-\widetilde{U}({{\boldsymbol{\xi}}}^{(2)}))}\right]$
$\displaystyle=\left[1-(\tau_{\delta}+\tau_{\delta}^{2})r_{\sigma}\right]^{\frac{n_{d}}{2}}\mathbb{E}\left[\widetilde{S}({{\boldsymbol{\xi}}}^{(1)},{{\boldsymbol{\xi}}}^{(2)})\right]$
$\displaystyle=S({{\boldsymbol{\xi}}}^{(1)},{{\boldsymbol{\xi}}}^{(2)}).$

When applying the m-repCNLD method, the term
$\widetilde{S}_{m}({{\boldsymbol{\xi}}}^{(1)},{{\boldsymbol{\xi}}}^{(2)})$
would be used to calculate the swapping rate.

###  3.2 Discrete adjoint method for gradient calculation

We note that the map between the unknowns ${\boldsymbol{\xi}}$ and output of
the PDEs is nonlinear, the gradient of the log-likelihood function with
respect to the unknown parameters is required, we use the adjoint method [1]
to calculate the gradient efficiently. Let the output of the forward model be
expressed in the following form.

${\bf G}({\boldsymbol{\xi}})={\bf CA}({\boldsymbol{\xi}})^{-1}{\bf
F}({\boldsymbol{\xi}}),$

where ${\bf C}\in{\mathbb{R}}^{n_{d}\times N}$ defines the linear state-to-
observation map and ${\bf A}({\boldsymbol{\xi}})^{-1}{\bf
F}({\boldsymbol{\xi}})$ is the solution of the discretized PDE, where ${\bf
A}({\boldsymbol{\xi}})\in{\mathbb{R}}^{N\times N}$ is assumed to be invertible
and Fréchet differentiable [1], ${\bf F}({\boldsymbol{\xi}})$ is also assumed
to be Fréchet differentiable here.

Denote ${\bf r}(\boldsymbol{\xi})={\bf G}({\boldsymbol{\xi}})-{\bf d}$ and
${\bf J}(\boldsymbol{\xi})$ the Jacobian of ${\bf G}({\boldsymbol{\xi}})$,
recall that

$\psi(\boldsymbol{\xi})=-\log\pi({\bf d}|\boldsymbol{\xi})=\frac{n_{d}}{2}\log
2\pi\sigma_{o}^{2}+\frac{\|{\bf d}-{\bf
G}(\boldsymbol{\xi})\|_{2}^{2}}{2\sigma_{o}^{2}},$

we then have

$\nabla\psi(\boldsymbol{\xi})=\frac{1}{\sigma_{o}^{2}}{\bf
J}(\boldsymbol{\xi})^{T}{\bf r}(\boldsymbol{\xi}),$ (3.23)

notice that the $i$th row of ${\bf J}(\boldsymbol{\xi})^{T}{\bf
r}(\boldsymbol{\xi})$ can be expressed as

$\displaystyle{[{\bf J}(\boldsymbol{\xi})^{T}{\bf r}(\boldsymbol{\xi})]_{i}}$
$\displaystyle=\frac{\partial{\bf
G}({\boldsymbol{\xi}})}{\partial\xi_{i}}\cdot{\bf r}(\boldsymbol{\xi})$ (3.24)
$\displaystyle=-\left[{\bf CA}({\boldsymbol{\xi}})^{-1}\left({\bf
F}_{i}^{\prime}({\boldsymbol{\xi}})-{\bf
A}_{i}^{\prime}({\boldsymbol{\xi}}){\bf A}({\boldsymbol{\xi}})^{-1}{\bf
F}({\boldsymbol{\xi}})\right)\right]\cdot{\bf r}(\boldsymbol{\xi}),$

where the symbol $\cdot$ represents the Frobenius inner products, ${\bf
A}_{i}^{\prime}({\boldsymbol{\xi}}):=\frac{d}{d\tau}{\bf
A}({\boldsymbol{\xi}}+\tau{\bf e_{i}})|_{\tau=0}$, ${\bf
F}_{i}^{\prime}({\boldsymbol{\xi}}):=\frac{d}{d\tau}{\bf
F}({\boldsymbol{\xi}}+\tau{\bf e_{i}})|_{\tau=0}$. The second equality in the
equation $(\ref{irow})$ holds because

$\frac{\partial{\bf G}({\boldsymbol{\xi}})}{\partial\xi_{i}}={\bf
C}\frac{\partial{\bf A}^{-1}({\boldsymbol{\xi}})}{\partial\xi_{i}}{\bf
F}({\boldsymbol{\xi}})+{\bf C}{\bf A}^{-1}({\boldsymbol{\xi}}){\bf
F}_{i}^{\prime}({\boldsymbol{\xi}}),$ (3.25)

since

${\bf A}({\boldsymbol{\xi}}){\bf A}^{-1}({\boldsymbol{\xi}})=I,$

we have

$\frac{\partial{\bf A}^{-1}({\boldsymbol{\xi}})}{\partial\xi_{i}}=-{\bf
A}^{-1}({\boldsymbol{\xi}}){\bf A}_{i}^{\prime}({\boldsymbol{\xi}}){\bf
A}^{-1}({\boldsymbol{\xi}}),$

substitute it into the equation $(\ref{GD})$, the corresponding expression of
$\frac{\partial{\bf G}({\boldsymbol{\xi}})}{\partial\xi_{i}}$ is resulted.

Let ${\bf v}$ and ${\bf z}$ satisfy the forward model

${\bf A}({\boldsymbol{\xi}}){\bf v}={\bf F}({\boldsymbol{\xi}}),$ (3.26)

and the costate equation

${\bf A}({\boldsymbol{\xi}})^{T}{\bf z}=-{\bf C}^{T}{\bf
r}(\boldsymbol{\xi}),$ (3.27)

respectively. The equation $(\ref{irow})$ can then be arranged as

${[{\bf J}(\boldsymbol{\xi})^{T}{\bf r}(\boldsymbol{\xi})]_{i}}=\left[{\bf
A}_{i}^{\prime}({\boldsymbol{\xi}}){\bf v}-{\bf
F}_{i}^{\prime}({\boldsymbol{\xi}})\right]\cdot{\bf z},$

it means that for any ${\boldsymbol{\xi}}$, we need only solve the forward
model and the costate equation to compute $\nabla\psi(\boldsymbol{\xi})$.
Meanwhile, the matrices ${\bf A}_{i}^{\prime}({\boldsymbol{\xi}})$ and ${\bf
F}_{i}^{\prime}({\boldsymbol{\xi}})$ can be calculated in parallel. We list
the discrete adjoint-based repCNLD method in Algorithm 1, the for-loop can be
proceeded in parallel to save the computational cost. For the case of
m-repCNLD, the forward model and costate equations are solved with the low-
fidelity method for the second chain, the term in the swapping rate should be
modified as equation $(\ref{modiS})$.

Input: The mean $m$ and covariance matrix $B$ of the prior, parameter $\beta$,
$\tau_{1}$, $\tau_{2}$, observed data ${\bf d}$, observation matrix $\bf C$,
noise of the observed data $\sigma_{o}$. The number of iterations $N_{max}$.
The intensity $c$ and time step $\delta$ are omitted in the swapping rate.

Output: Samples of Markov chains $\\{\Xi_{i}^{(1)}\\}_{i=1}^{N_{max}}$ and
$\\{\Xi_{i}^{(2)}\\}_{i=1}^{N_{max}}$

1 Generate the initial values $\xi_{0}^{(1)}$, $\xi_{0}^{(2)}$ for the two
chains randomly.

2Evaluate the forward model at $\xi_{0}^{(j)}$ by equation $(\ref{Forward})$,
obtain the solution ${\bf v}_{0}^{(j)}$ and errors ${\bf r}_{0}^{(j)}={\bf
C}{\bf v}_{0}^{(j)}-{\bf d}$, for $j=1,2$.

3while  _$i <N_{\max}$_ do

4 for  _$j=1:2$_ do

5

6 Solve the costate equation $(\ref{cos})$ and calculate the gradient term
$\nabla\psi(\xi_{0}^{(j)})$ by equation $(\ref{gradient})$.

7 Update the new particles by
$\xi^{(j)}=\sqrt{1-\beta^{2}}\xi_{0}^{(j)}+(1-\sqrt{1-\beta^{2}})(m-B\nabla\psi(\xi_{0}^{(j)}))+\beta\sqrt{\tau_{j}}w^{(j)},\
w^{(j)}\sim\mathcal{N}(0,B).$

8 Simulate the forward model at $\xi^{(j)}$ by equation $(\ref{Forward})$,
obtain the solution ${\bf v}^{(j)}$ and error ${\bf r}^{(j)}={\bf C}{\bf
v}^{(j)}-{\bf d}$.

9 Calculate the energy function $E^{(j)}=U(\xi^{(j)})$ according to equation
$(\ref{energy})$.

10

11 Draw $u$ from the uniform distribution $\mathcal{U}[0,1]$ and calculate the
exchange rate
$s=\min\bigg{\\{}1,\exp\left(\tau_{\delta}(E^{(1)}-E^{(2)})\right)\bigg{\\}},$

12 if  _$u <s$_ then

13 Update $\xi_{0}^{(1)}$ by $\xi^{(2)}$, $\xi_{0}^{(2)}$ by $\xi^{(1)}$.

14 Update ${\bf v}_{0}^{(1)}$ by ${\bf v}^{(2)}$, ${\bf v}_{0}^{(2)}$ by ${\bf
v}^{(1)}$.

15 Update ${\bf r}_{0}^{(1)}$ by ${\bf r}^{(2)}$, ${\bf r}_{0}^{(2)}$ by ${\bf
r}^{(1)}$.

16 else

17 Update $\xi_{0}^{(j)}$ by $\xi^{(j)}$, for $j=1,2$.

18 Update ${\bf v}_{0}^{(j)}$ by ${\bf v}^{(j)}$, for $j=1,2$.

19 Update ${\bf r}_{0}^{(j)}$ by ${\bf r}^{(j)}$, for $j=1,2$.

20 Save the samples as $\Xi_{i}^{(1)}=\xi_{0}^{(1)},\ \
\Xi_{i}^{(2)}=\xi_{0}^{(2)}.$

Algorithm 1 Discrete adjoint-based repCNLD method.

##  4 Numerical experiments

In this section, we present some numerical experiments results. In particular,
we apply the proposed repCNLD method to sample Gaussian mixture distribution
and solve Bayesian nonlinear inverse problems. The Gaussian mixture
distribution simulation shows that the method can capture multi-mode but not
get stuck in local modes. For the nonlinear PDE inverse problems, the adjoint
method is used to compute the gradient of the logarithmic likelihood function
effectively. The proposed methods are able to explore the multi-modal and
high-dimensional posterior density, meanwhile improving the efficiency of the
MCMC method.

###  4.1 Gaussian mixture density simulation

In this subsection, we will show the performance of the proposed method in
sampling the Gaussian mixture density, i.e., for
${\boldsymbol{\xi}}\in{\mathbb{R}}^{n}$, the target distribution is

$\pi({\boldsymbol{\xi}}|\theta)=\sum_{k=1}^{M}\gamma_{k}\phi({\boldsymbol{\xi}}|\theta_{k}),$

where $\gamma_{k}\geq 0$ is the weight coefficients and
$\sum_{k=1}^{M}\gamma_{k}=1$, $\phi({\boldsymbol{\xi}}|\theta_{k})$ is the
Gaussian density and has the form

$\phi({\boldsymbol{\xi}}|\theta_{k})=(2\pi)^{-\frac{n}{2}}|B_{k}|^{-\frac{n}{2}}\exp\bigg{\\{}-\frac{1}{2}({\boldsymbol{\xi}}-m_{k})^{T}B_{k}^{-1}({\boldsymbol{\xi}}-m_{k})\bigg{\\}},$

where $\theta_{k}=(m_{k},B_{k})$ is the given parameters of the density. We
use the scheme $(\ref{pCN-langevin})$ and $(\ref{repCNLDs})$ for the pCNLD and
repCNLD method, respectively. The potential function is set as

$\psi({\boldsymbol{\xi}})=-\log\pi({\boldsymbol{\xi}}|\theta)-\frac{1}{2}(\boldsymbol{\xi}-m)^{T}B^{-1}(\boldsymbol{\xi}-m),$
(4.28)

such that $\exp\\{-U({\boldsymbol{\xi}})\\}$ equals to the target distribution
$\pi({\boldsymbol{\xi}}|\theta)$. Substitute the equation
$(\ref{m-potential})$ into the equation $(\ref{pCN-langevin})$, the updating
scheme for a single chain can be rearranged as

${\boldsymbol{\xi}}_{k+1}={\boldsymbol{\xi}}_{k}+(1-\sqrt{1-\beta^{2}})B\nabla\log\pi({\boldsymbol{\xi}}_{k}|\theta)+\beta\sqrt{B\tau}{\boldsymbol{w}}_{k},$

the scheme for the repCNLD method can be correspondingly derived. It can be
seen that the schemes does not depend on the mean of the prior, the matrix $B$
works as a scalable parameter here.

####  4.1.1 1-dimensional example

In this experiment, we demonstrate the repCNLD method on a 1-dimensional
Gaussian mixture model with ${M}=2$, $\gamma_{1}=0.4$, $\gamma_{2}=0.6$,
$\theta_{1}=(-3,0.7^{2})$, $\theta_{2}=(2,0.5^{2})$, which is the first
example in paper [11]. The time steps are set as $\delta=0.001$, the
corresponding values of $\beta$ are calculated as $\beta=0.0447$, and the
temperatures $\tau_{1}$ and $\tau_{2}$ are set to be 1 and 15, respectively.
We also use the pCNLD method to simulate the target distribution as a
comparison; the time step and temperature are set to be $\delta_{1}$ and
$\tau_{1}$. The auxiliary prior is chosen as $\mathcal{N}(0,3)$ for both
methods. We simulate $10^{5}$ samples from each distribution and the kernel
estimations of the samples are displayed in Figure 4.1(Left). As can be seen,
the pCNLD gets stuck in a single mode, while repCNLD successfully yields a
close approximation to the ground truth distribution.

We also simulate the distribution with $\mu_{1}=-6$, $\mu_{2}=4$, In this
example, the distance between the mean of two distributions is significantly
larger than in the first example. As a result, it will be harder to capture
both modes. The temperature for the second chain of the repCNLD is set
$\tau_{2}=40$, the auxiliary prior is set $\mathcal{N}(0,3^{2})$, of which the
range is wider, and the other parameters are kept the same. The results are
shown in Figure 4.1(Right), the repCNLD can still yield a good approximation,
and the pCNLD still fails to escape from the local modes as expected.

Figure 4.1: Kernel density estimations for samples obtained by the pCNLD and
repCNLD methods with: (Left) $\mu_{1}=-3,\mu_{2}=2$ and (Right)
$\mu_{1}=-6,\mu_{2}=4$. The solid line is the true density.

Figure 4.2: Density contour plot of (Left) samples obtained by the repCNLD
method and (Right) samples obtained by the pCNLD method. The ground true
density is plotted in the solid line in both pictures, the density of samples
obtained by the repCNLD method is plotted in the dashed line, while the
density of samples obtained by the pCNLD method is shown in dot.

####  4.1.2 2-dimensional example

In this example, we demonstrate the repCNLD method on a 2-dimensional Gaussian
mixture model with $M=3$, $\gamma_{1}=0.3$, $\gamma_{2}=0.3$,
$\gamma_{3}=0.4$, $m_{1}=[4,2]^{T}$, $m_{2}=[-4,2]^{T}$, $m_{3}=[0,-3]^{T}$,
and

$\begin{gathered}B_{1}=\begin{bmatrix}1&0.6\\\ 0.6&1\end{bmatrix},\quad
B_{2}=\begin{bmatrix}1&-0.6\\\ -0.6&1\end{bmatrix},\quad
B_{3}=\begin{bmatrix}1&0\\\ 0&1\end{bmatrix}.\end{gathered}$

The time steps are set as $\delta=0.001$, the temperatures $\tau_{1}$ and
$\tau_{2}$ are set to be 1 and 20 for the repCNLD method, while the time step
and temperature are set as $\delta$ and $\tau_{1}$ for the single chain pCNLD
method. The mean and covariance matrix of the Gaussian prior is set as
$[0,0]^{T}$ and $B=\left[\begin{array}[]{cc}10&0\\\ 0&10\end{array}\right]$ .
We draw $10^{5}$ samplers for both the repCNLD and pCNLD methods, the contour
plot of the kernel density estimations are illustrated in Figure 4.2. The
samples drawn by the repCNLD method provide a good approximation of the target
distribution; they explore all modes of the target distribution. The samples
obtained by the pCNLD get trapped in local modes with mean $m_{1}$ and
$m_{3}$, in addition, few samples are exploring the mode with mean $m_{2}$.

###  4.2 Initial center identification

In this section, we consider apply the repCNLD method to nonlinear inverse
problems whose underlying forward model is the single-phase Darcy flow. The
flow is described in terms of concentration $u(x,t)$, which is governed by the
following parabolic equation

$\frac{\partial u(x,t)}{\partial t}=\text{div}(\kappa(x)\nabla
u(x,t))+f(x,t),\ x\in\Omega,t\in(0,T]\\\ $ (4.29)

subject to an appropriate boundary condition and initial condition. The
studied physical domain is $\Omega=[0,1]\times[0,1]$, the terminal time is set
as $T=0.1$.

The initial condition has the form

$u_{0}(x;{\boldsymbol{\xi}})=\frac{q}{2\pi
l_{x_{1}}l_{x_{2}}}\exp\left(-\frac{(x_{1}-\xi_{1})^{2}}{2l_{x_{1}}^{2}}-\frac{(x_{2}-\xi_{2})^{2}}{2l_{x_{2}}^{2}}\right),$

the unknown parameter ${\boldsymbol{\xi}}=(\xi_{1},\xi_{2})$ is the location
of the pollution, the other parameters are prescribed as

$l_{x_{1}}=0.1,\ l_{x_{2}}=0.2,\ q=1.$

The accurate solution of the forward problem is set to be

$u(x,t;{\boldsymbol{\xi}})=u_{0}(x;{\boldsymbol{\xi}})e^{-t},$

$\kappa(x)=1$, the inhomogeneous Dirichlet boundary condition and source term
are preset to match the solution. Our goal is to infer ${\boldsymbol{\xi}}$
from some measurements of $u(x,t)$ at time $T$. We observe the data at the
place $x_{o}=(0.5,0.3)$, the measurement error is set as $\sigma_{o}=0.1$. In
this experiment, $u(x_{o},T)=Me^{-r^{2}}e^{-T}$, where $r=1$. This setting
guarantees an infinite number of solutions to the inverse problem [24]. We set
the mean and covariance of the prior as $m=[1/2,1/2]^{T}$ and
$B=\left[\begin{array}[]{cc}1/4&0\\\ 0&1/4\end{array}\right]$ to guarantee the
domain $\Omega$ be sampled with a certain probability.

In order to show the efficiency and accuracy of the methods, we compare the
pCNLD, the repCNLD proceeded with the same forward solver, and the m-repCNLD.
We use the standard finite element method to solve the forward models, the
backward Euler method is used for the temporal discretization, and the time
step is set to be $\Delta t=0.001$ for all methods. For the m-repCNLD, the
solver with $\Delta x=1/40$ is assigned to the chain with low temperature,
while in the high-temperature chain, the forward model is solved with $\Delta
x=1/20$. $\Delta x=1/40$ is set for both the single chain pCNLD and the
repCNLD methods.

The time steps used in discretizing the Langevin diffusion chains are set as
$\delta=10^{-5}$, the temperatures are set as $\tau_{1}=1$, $\tau_{2}=15$ for
both the repCNLD and m-repCNLD methods. For the pCNLD method, the parameters
are prescribed as $\delta$ and $\tau_{1}$. Starting with the same initial
value, chains with $3\times 10^{5}$ iterations are run. We note that running
the chain for $10^{5}$ iterations, the pCNLD simulation time is about
$41462.9$s, while the repCNLD and m-repCNLD methods take $49549.6$s and
$47660.7$s to complete all iterations, respectively. Actually, when the number
of iterations is set the same, simulations for the repCNLD or m-repCNLD
methods double the simulations for the single chain pCNLD method. However, the
computational cost of the repCNLD method is much less than twice the
consumption of the single-chain pCNLD method; this is because the samples from
the two chains of the repCNLD or m-repCNLD are updated in parallel. Meanwhile,
the parallel technique makes the speeding-up effects of the m-repCNLD not
significant. If the parallel toolbox is unavailable, the multi-variance
strategy can significantly save the computational cost significantly [24].

Figure 4.3: Autocorrelation functions (ACF) of the Onsager-Machlup functional
(OMF) or the energy function for different methods. The values of ACF for the
replica exchange methods decrease dramatically to 0, and the ACF for the
single chain method decreases more slowly compared with the double chain
methods. Figure 4.4: Effective sample size of the posterior samples for
different methods. With the same length of Markov chain, the independent
samples obtained by the pCNLD method is less than the ones obtained by the
replica exchange methods for both $\xi_{1}$ and $\xi_{2}$.  Figure 4.5:
Scatter plot of posterior samples for (Left) pCNLD method, (Middle) repCNLD
method, and (Right) m-repCNLD method, respectively. The solid line is the true
solution to the inverse problem.

We use the autocorrelation functions (ACF) of the Onsager-Machlup functional
(OMF)[10, 22] to compare the efficiency of the three mentioned methods, the
OMF here is the corresponding energy function we defined in equation (2.3).
The ACF of the OMF for different methods are displayed as functions of the lag
in Figure 4.3. It can be seen from the figure that all of the methods result
in low ACF values. The ACF values of repCNLD and m-repCNLD methods decrease
much more dramatically than the ACF of the pCNLD method, the values are also
much lower. Both the repCNLD and m-repCNLD methods perform better than the
pCNLD method.

The effective sample size (ESS) is another scholar to measure the sampling
efficiency of MCMC [18]. The most common definition of ESS is

$\text{ESS}=\frac{N_{max}}{1+2\rho},$

where $\rho$ is the integrated autocorrelation time and $N_{max}$ is the total
sample size, and it is used to estimate the number of effectively independent
samples in the Markov chain. With the same length of the Markov chain, the
bigger the ESS is, the more effective the algorithm is. We compute the ESS of
the unknown $\xi$ and show them in Figure 4.4, the bar-plot indicates that the
replica exchange method generates much more effectively independent draws than
the single chain method, while the repCNLD method outperforms the m-repCNLD
method in both $\xi_{1}$ and $\xi_{2}$.

In Figure 4.5, the last $3\times 10^{4}$ samples of the posterior are plotted
for each method, and the solid line is the true solution. We can see that the
pCNLD method fails to capture all of the solutions, while the repCNLD and
m-repCNLD methods capture all of the solutions to the inverse problem. The
second chain of the replica exchange methods successfully helps the particles
in the first chain escape from the local modes as expected. In summary, in
this inverse problem with infinite solutions, the replica exchange methods are
significantly more efficient than the single-chain method.

###  4.3 Permeability identification

In this subsection, we consider the permeability identification of the flow.
The Gaussian process is imposed as the prior of the log permeability.
Moreover, we adopt the truncated Karhunen Loève expansion to parameterize the
Gaussian field $Y(x,\omega)=\log(\kappa(x;\boldsymbol{\xi}))$ and it can be
expressed as

$Y(x,\omega)=\bar{Y}(x)+\sum_{i=1}^{n}\sqrt{\lambda_{i}}\varphi(x)\xi_{i}(\omega),$

where $\bar{Y}(x)$ is the mean of the random field,
$\\{\xi_{i}(\omega)\\}_{i=1}^{n}$ are uncorrelated random variables, and
$(\lambda_{i},\varphi_{i}(x))$ are the eigenvalue and eigenfunctions of the
eigenvalue problem

$\int_{\Omega}\text{cov}[Y](x_{1},x_{2})\varphi(x_{2})dx_{2}=\lambda\varphi(x_{1}),$

where $\text{cov}[Y]$ is the covariance function of $Y(x,\omega)$. In this
experiment, the mean of the random field is set as $\bar{Y}(x)=0$, the
covariance function of $Y(x,\omega)$ is given by [34]

$\text{cov}[Y](x_{1},x_{2})=\sigma^{2}_{Y}\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}\frac{|x_{1}-x_{2}|}{l}\right)^{\nu}B_{\nu}\left(\sqrt{2\nu}\frac{|x_{1}-x_{2}|}{l}\right)$

where $\Gamma(\cdot)$ is the Gamma function, and $B_{\nu}(\cdot)$ is the
modified Bessel function. We choose $\sigma^{2}_{Y}=1$, $\nu=0.2$, the length
scale along the horizontal and vertical direction is set as 1 and 0.5,
respectively. $86\%$ of the total prior energy defined by

$e(n)=\frac{\sum_{i=1}^{n}\lambda_{i}}{\sum_{i=1}^{\infty}\lambda_{i}}$

is retained, where the eigenvalues are sorted in ascending order, resulting in
$n=15$.

The force term in the equation $(\ref{ex1})$ is set as

$f(x,t)=\sum_{i=1}^{N_{w}}\frac{q_{i}}{2\pi
l_{w}^{2}}\exp\bigg{\\{}-\frac{\|x-x_{i}\|^{2}}{2l_{w}^{2}}\bigg{\\}},$

where $\\{q_{i}\\}_{i=1}^{N_{w}}$ and $\\{x_{i}\\}_{i=1}^{N_{w}}$ are the
prescribed production rates and well locations [6]. The production rate is set
as 1.5 and $l_{w}=0.05$ in the force term. The true coefficient
$\kappa(x;\boldsymbol{\xi})$ is a realization from the prior, the true log
permeability and the well locations are shown in Figure 4.6. In addition, we
consider the initial condition $u(x,0)=4$ and no-flow boundary conditions for
this problem.

The synthetic data is generated by using FEM on a $120\times 120$ grid and
100-time steps. The observed data is the pressure collected at the well
locations, the measurement time used are $t_{1}=0.01$, $t_{n}=0.02n$
($n=2,\cdots,5$), the observed noise is $\sigma_{o}=0.1$. The forward model is
solved with $\Delta x=1/60$, and the time interval $[0,T]$ is discretized with
50 steps for the repCNLD and pCNLD methods. For the m-repCNLD method, the
forward model is solved with $\Delta t=0.002$ for the first chain, while the
forward model simulated in the second chain is run by $\Delta t=0.005$.

Figure 4.6: (Left) The true log permeability, (Right) well configuration.

The time steps involved in the Langevin dynamics are set as $\delta=2\times
10^{-3}$, the temperatures are set as $\tau_{1}=1$, $\tau_{2}=1.5$ for both
the repCNLD and m-repCNLD methods. The parameters are prescribed as
$\delta_{1}$ and $\tau_{1}$ for the pCNLD method. The starting point is set as
the mean of the prior. The computational cost of $10^{5}$ iterations for the
pCNLD, repCNLD and m-repCNLD methods is $57263.6$s, $58535.2$s, $57212.7$s,
respectively. This is also due to the usage of the adjoint method and parallel
technique. $5\times 10^{5}$ samples are drawn for each method, the ACF of the
OMF is plotted against the lag in Figure 4.7, the ACF values of the repCNLD
and m-repCNLD methods are lower than the pCNLD method, which means that they
have higher mixing rate than the pCNLD method.

Figure 4.7: Autocorrelation functions (ACF) of the Onsager-Machlup functional
(OMF) or the energy function for different methods. The values of ACF for the
replica exchange methods are lower than the ones of the pCNLD method, the
chains obtained by the replica exchange methods have a higher mixing rate than
the ones obtained by the pCNLD method.  Figure 4.8: The 1-D and 2-D marginal
posterior probability distributions for part of the parameters for the pCNLD
method. The uncertainty of the modes $[\xi_{1},\xi_{3},\xi_{5}]$ is reduced by
the observed data, while for the modes that are associated with small
eigenvalues, the data is not informative. Figure 4.9: The 1-D and 2-D marginal
posterior probability distributions for part of the parameters for the repCNLD
method. Figure 4.10: The 1-D and 2-D marginal posterior probability
distributions for part of the parameters for the m-repCNLD method.

We visualize the 1-dimensional (1-D) and 2-dimensional (2-D) marginal
posterior probability distributions for part of the parameters for the three
methods. Based on the shape of their 2-D marginals, some correlations are
apparent between several modes, such as $\xi_{1}$ and $\xi_{7}$, $\xi_{3}$ and
$\xi_{5}$, $\xi_{7}$ and $\xi_{13}$. The ranges of the 1-D marginal posterior
for parameters $\xi_{1}$, $\xi_{3}$ and $\xi_{5}$ are much narrower than the
ones of the high-indexed modes, which indicates that the infinite data is not
informative about the small scales of the unknown log permeability. It can be
observed from Figure 4.11, the posterior mean and standard deviation for all
the methods have common spatial features. There is more uncertainty around the
wells $P_{1}$ and $P_{4}$, and less uncertainty around the other wells. The
non-zero values of the posterior skewness indicate that the posterior is non-
Gaussian; this may be contributed by the non-linearity of the forward map and
the insufficient data information.

Figure 4.11: The estimated posterior (Upper) mean, (Middle) standard deviation
and (Bottom) skewness of the log $\kappa(x;\xi)$ for different methods. The
spatial distributions of estimates are basically the same. The skewness values
around the wells $P_{1}$ and $P_{7}$ are larger than the ones around the other
wells, since the data provided by the two wells is less informative.

##  5 Conclusions

In this work, we propose the repCNLD MC algorithm to explore target or
posterior distributions; two preconditioned Langevin diffusion chains with
different temperatures are simulated for the sampling process. The diffusion
with low-temperature aims to exploit the geometry locally, while the diffusion
with high temperature works to explore the domain globally. They are allowed
to swap with the specified probability. The swapping does not change the
invariant distribution of the diffusion chains but accelerates the convergence
of the target distribution. We use the Crank-Nicolson scheme to discretize the
preconditioned Langevin diffusion chains; the covariance of the Gaussian prior
is cooperated with the temperature parameters to scale the noise of the
diffusion together. We analyze the discretization error of using the pCN
scheme. The proposed method performs well for multimodal distribution taking
advantage of the replica exchange method. In order to save the computational
cost, the forward model applied in the second chain is solved with the low-
fidelity method; we derive an unbiased estimator for the swapping exchange
rate by tracing to the forward model errors. Meanwhile, we provide the
framework of a discrete adjoint method for the calculation of the gradient. We
demonstrate the efficiency of the proposed method through some numerical
examples. For the nonlinear inverse problems, we use the ACF values of OMF
(energy function) values to assess the MCMC chains’ efficiency. The replica
exchange methods outperform the single-chain methods with the same starting
points.

We use the adjoint method to calculate the derivative information of the log-
likelihood function, which saves the computational cost. However, it would be
unattainable for problems with black-box models. The covariance matrix of the
prior works is a preconditioned matrix here; it is preset and fixed; some
adaptive methods may be used to improve the efficiency of our proposed method
to deal with complex problems where the scales of the objective function are
different in different modes. Further investigation of these issues is worth
pursuing in the future.

## Acknowledgments

NO acknowledges the support of Hunan Provincial Natural Science Foundation of
China 2021JJ40557 and Chinese NSF 11901060. GL gratefully acknowledges the
support of the National Science Foundation (DMS-1555072, DMS-2053746, and
DMS-2134209), and Brookhaven National Laboratory Subcontract 382247, and U.S.
Department of Energy (DOE) Office of Science Advanced Scientific Computing
Research program DE-SC0021142 and DE-SC0023161.

## References

  * [1] Johnathan M Bardsley.  Computational uncertainty quantification for inverse problems, volume 19.  SIAM, 2018.
  * [2] Alexandros Beskos, Gareth Roberts, and Andrew Stuart.  Optimal scalings for local metropolis–hastings chains on nonproduct targets in high dimensions.  The Annals of Applied Probability, 19(3):863–898, 2009.
  * [3] Alexandros Beskos, Gareth Roberts, Andrew Stuart, and Jochen Voss.  Mcmc methods for diffusion bridges.  Stochastics and Dynamics, 8(03):319–350, 2008.
  * [4] Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng.  Handbook of markov chain monte carlo.  CRC press, 2011.
  * [5] Yi Chen, Jinglin Chen, Jing Dong, Jian Peng, and Zhaoran Wang.  Accelerating nonconvex learning via replica exchange langevin diffusion.  arXiv preprint arXiv:2007.01990, 2020.
  * [6] Zhangxin Chen, Guanren Huan, and Yuanle Ma.  Computational methods for multiphase flows in porous media.  SIAM, 2006.
  * [7] Eric Chung, Yalchin Efendiev, Wing Tat Leung, Sai-Mang Pun, and Zecheng Zhang.  Multi-agent reinforcement learning accelerated mcmc on multiscale inversion problem.  arXiv preprint arXiv:2011.08954, 2020.
  * [8] Simon L Cotter, Gareth O Roberts, Andrew M Stuart, and David White.  Mcmc methods for functions: modifying old algorithms to make them faster.  Statistical Science, 28(3):424–446, 2013.
  * [9] Arnak S Dalalyan and Lionel Riou-Durand.  On sampling from a log-concave density using kinetic langevin diffusions.  Bernoulli, 26(3):1956–1988, 2020.
  * [10] Masoumeh Dashti, Kody JH Law, Andrew M Stuart, and Jochen Voss.  Map estimators and their consistency in bayesian nonparametric inverse problems.  Inverse Problems, 29(9):095017, 2013.
  * [11] Wei Deng, Qi Feng, Liyao Gao, Faming Liang, and Guang Lin.  Non-convex learning via replica exchange stochastic gradient mcmc.  In International Conference on Machine Learning, pages 2474–2483. PMLR, 2020.
  * [12] Alain Durmus and Eric Moulines.  High-dimensional bayesian inference via the unadjusted langevin algorithm.  Bernoulli, 25(4A):2854–2882, 2019.
  * [13] Raaz Dwivedi, Yuansi Chen, Martin J Wainwright, and Bin Yu.  Log-concave sampling: Metropolis-hastings algorithms are fast!  In Conference on learning theory, pages 793–797. PMLR, 2018.
  * [14] Y. Efendiev, T. Hou, and W. Luo.  Preconditioning markov chain monte carlo simulations using coarse-scale models.  SIAM Journal on Scientific Computing, 28(2):776–803, 2006.
  * [15] Joel N Franklin.  Well-posed stochastic extensions of ill-posed linear problems.  Journal of mathematical analysis and applications, 31(3):682–716, 1970.
  * [16] Nils Lid Hjort, Chris Holmes, Peter Müller, and Stephen G Walker.  Bayesian nonparametrics, volume 28.  Cambridge University Press, 2010.
  * [17] M. Izzatullah, T. V. Leeuwen, and D. Peter.  Bayesian seismic inversion: A fast sampling langevin dynamics markov chain monte carlo method.  Geophysical Journal International, 2021.
  * [18] Robert E Kass, Bradley P Carlin, Andrew Gelman, and Radford M Neal.  Markov chain monte carlo in practice: a roundtable discussion.  The American Statistician, 52(2):93–100, 1998.
  * [19] PE Kloeden and S Shott.  Linear-implicit strong schemes for itô-galkerin approximations of stochastic pdes.  Journal of Applied Mathematics and Stochastic Analysis, 14(1):47–53, 2001.
  * [20] V. Leeuwen and P. Jan.  Nonlinear data assimilation.  Frontiers in Applied Dynamical Systems Reviews & Tutorials, 2, 2015\.
  * [21] Don S Lemons and Anthony Gythiel.  Paul langevin’s 1908 paper “on the theory of brownian motion”[“sur la théorie du mouvement brownien,” cr acad. sci.(paris) 146, 530–533 (1908)].  American Journal of Physics, 65(11):1079–1081, 1997.
  * [22] Jinglai Li.  A note on the karhunen–loève expansions for infinite-dimensional bayesian inverse problems.  Statistics & Probability Letters, 106:1–4, 2015.
  * [23] Guang Lin, Christian Moya, and Zecheng Zhang.  Accelerated replica exchange stochastic gradient langevin diffusion enhanced bayesian deeponet for solving noisy parametric pdes.  arXiv preprint arXiv:2111.02484, 2021.
  * [24] Guang Lin, Yating Wang, and Zecheng Zhang.  Multi-variance replica exchange stochastic gradient mcmc for inverse and forward bayesian physics-informed neural network.  arXiv preprint arXiv:2107.06330, 2021.
  * [25] Jun S Liu and Jun S Liu.  Monte Carlo strategies in scientific computing, volume 10.  Springer, 2001.
  * [26] Jonathan C Mattingly, Natesh S Pillai, and Andrew M Stuart.  Diffusion limits of the random walk metropolis algorithm in high dimensions.  The Annals of Applied Probability, 22(3):881–930, 2012.
  * [27] Bernt Oksendal.  Stochastic differential equations: an introduction with applications.  Springer Science & Business Media, 2013.
  * [28] K. Platen and E. Peter.  Numerical solution of stochastic differential equations.  quantitative finance, 1992.
  * [29] Christian Robert and George Casella.  Monte Carlo statistical methods.  Springer Science & Business Media, 2013.
  * [30] Christian P Robert, George Casella, and George Casella.  Monte Carlo statistical methods, volume 2.  Springer, 1999.
  * [31] Gareth O Roberts and Richard L Tweedie.  Exponential convergence of langevin distributions and their discrete approximations.  Bernoulli, pages 341–363, 1996.
  * [32] Michael L Stein.  Interpolation of spatial data: some theory for kriging.  Springer Science & Business Media, 1999.
  * [33] Andrew M Stuart, Jochen Voss, and Petter Wilberg.  Conditional path sampling of sdes and the langevin mcmc method.  Communications in Mathematical Sciences, 2(4):685–697, 2004.
  * [34] Christopher KI Williams and Carl Edward Rasmussen.  Gaussian processes for machine learning.  MIT press Cambridge, MA, 2006.

## Appendix

### Proof of theorem 2.1

###### Proof.

The dynamics in equations $(\ref{cns-U})$ and $(\ref{swap})$ defines a Markov
jump process, following [27] and [5], the infinitesimal generator has the form

$\begin{array}[]{l}\mathcal{L}^{c}\left(f(x,y)\right)=\underbrace{-[L\nabla_{x}U(x)]\cdot\nabla_{x}f(x,y)+\tau_{1}L:\nabla_{x}\nabla_{x}f(x,y)}_{\mathcal{L}^{c}_{1}\left(f(x,y)\right)}\\\
\underbrace{-[L\nabla_{y}U(y)]\cdot\nabla_{y}f(x,y)+\tau_{2}L:\nabla_{y}\nabla_{y}f(x,y)}_{\mathcal{L}^{c}_{2}\left(f(x,y)\right)}+\underbrace{cs(x,y)\left(f(y,x)-f(x,y)\right)}_{\mathcal{L}^{c}_{s}\left(f(x,y)\right)},\end{array}$

As demonstrated in the proof of the Lemma 3.2 in paper [5], we only need to
show that

$\int\int g(x,y)\mathcal{L}^{c}\left(f(x,y)\right)\pi(x,y){\rm d}x{\rm
d}y=\int\int f(x,y)\mathcal{L}^{c}\left(g(x,y)\right)\pi(x,y){\rm d}x{\rm d}y$
(5.30)

for all $f,g\in\mathcal{D}(\mathcal{L}^{c})$, where the density $\pi(x,y)$ is
defined as $(\ref{invariant})$, to prove the invariance and reversibility.
Since the term arises from the swap $\mathcal{L}_{s}^{c}$ is the same with the
one in paper [5], we directly have

$\int\int g(x,y)\mathcal{L}_{s}^{c}\left(f(x,y)\right)\pi(x,y){\rm d}x{\rm
d}y=\int\int f(x,y)\mathcal{L}_{s}^{c}\left(g(x,y)\right)\pi(x,y){\rm d}x{\rm
d}y.$

For the term $\mathcal{L}_{1}^{c}$, by integration by parts and the fact that
$f$ and $g$ have compact support, we have that for all fixed
$y\in\mathbb{R}^{n}$,

$\begin{array}[]{l}{\int}\tau_{1}g(x,y)\left[L:\nabla_{x}\nabla_{x}f(x,y)\right]\exp(\frac{U(x)}{\tau_{1}}){\rm
d}x\\\
=-\tau_{1}{\int}L:\left(\nabla_{x}f(x,y)\left(\nabla_{x}\left[g(x,y)\exp(\frac{U(x)}{\tau_{1}})\right]\right)^{T}\right){\rm
d}x\\\
=-{\int}\tau_{1}\exp(\frac{U(x)}{\tau_{1}})\left[L:\left(\nabla_{x}f(x,y)\left(\nabla_{x}g(x,y)\right)^{T}\right)\right]{\rm
d}x\\\
\quad+{\int}g(x,y)\exp(\frac{U(x)}{\tau_{1}})\left[L:\left(\nabla_{x}f(x,y)\left(\nabla_{x}U(x)\right)^{T}\right)\right]{\rm
d}x,\\\ \end{array}$

where the symbol $T$ means the transposition of a matrix. Note that the term
$[L\nabla_{x}U(x)]\cdot\nabla_{x}f(x,y)$ can be rewritten as
$L:\left(\nabla_{x}f(x,y)\left(\nabla_{x}U(x)\right)^{T}\right)$, we then have

$\int g(x,y)\mathcal{L}_{1}^{c}\left(f(x,y)\right)\pi(x,y){\rm
d}x=-\int\tau_{1}\left[L:\left(\nabla_{x}f(x,y)\left(\nabla_{x}g(x,y)\right)^{T}\right)\right]\pi(x,y){\rm
d}x,$

hence,

$\int\int g(x,y)\mathcal{L}_{1}^{c}\left(f(x,y)\right)\pi(x,y){\rm d}x{\rm
d}y=-\int\int\tau_{1}\left[L:\left(\nabla_{x}f(x,y)\left(\nabla_{x}g(x,y)\right)^{T}\right)\right]\pi(x,y){\rm
d}x{\rm d}y.$

By switching the positions of $f$ and $g$ in the above equation, we have

$\int\int f(x,y)\mathcal{L}_{1}^{c}\left(g(x,y)\right)\pi(x,y){\rm d}x{\rm
d}y=-\int\int\tau_{1}\left[L:\left(\nabla_{x}g(x,y)\left(\nabla_{x}f(x,y)\right)^{T}\right)\right]\pi(x,y){\rm
d}x{\rm d}y,$

since the matrix $L$ is symmetric, then

$L:\left(\nabla_{x}f(x,y)\left(\nabla_{x}g(x,y)\right)^{T}\right)=L:\left(\nabla_{x}g(x,y)\left(\nabla_{x}f(x,y)\right)^{T}\right),$

we finally have

$\int\int g(x,y)\mathcal{L}_{1}^{c}\left(f(x,y)\right)\pi(x,y){\rm d}x{\rm
d}y=\int\int f(x,y)\mathcal{L}_{1}^{c}\left(g(x,y)\right)\pi(x,y){\rm d}x{\rm
d}y.$ (5.31)

By the same derivation, the equation $(\ref{ig1})$ also holds for
$\mathcal{L}_{2}^{c}$. In summary, the equation $(\ref{detailbalance})$ is
proved to be established for the density $\pi$ defined in $(\ref{invariant})$.
Thus, $\\{\boldsymbol{\xi}_{t}\\}_{t\geq 0}$ is a reversible Markov process,
and its invariant distribution has the density
$\pi(\boldsymbol{\xi}^{(1)},\boldsymbol{\xi}^{(2)})\propto\exp\\{-\frac{U(\boldsymbol{\xi}^{(1)})}{\tau_{1}}-\frac{U(\boldsymbol{\xi}^{(2)})}{\tau_{2}}\\}$.
∎

### Proof of theorem 2.2

###### Proof.

According to the property of $\nabla$, we have

$\displaystyle\nabla_{z}f^{2}$ $\displaystyle=$ $\displaystyle 2f\nabla_{z}f,$
$\displaystyle\nabla_{z}\nabla_{z}f^{2}$ $\displaystyle=$ $\displaystyle
2\left(\nabla_{z}f\left(\nabla_{z}f\right)^{T}\right)+2f\left(\nabla_{z}\nabla_{z}f\right)$

for $z=x,y$. Then we have

$\begin{array}[]{l}\Gamma^{c}\left(f\left(x,y\right)\right)=\frac{1}{2}\mathcal{L}^{c}\left(f^{2}\left(x,y\right)\right)-f\left(x,y\right)\mathcal{L}^{c}\left(f\left(x,y\right)\right)\\\
=\tau_{1}\left[L\nabla_{x}f(x,y)\right]\cdot\nabla_{x}f(x,y)+\tau_{2}\left[L\nabla_{y}f(x,y)\right]\cdot\nabla_{y}f(x,y)\\\
\quad+\frac{c}{2}s(x,y)\left(f(y,x)-f(x,y)\right)^{2},\\\ \end{array}$

since the matrix $L$ is symmetric positive definite, the first two terms of
the above equation is consequently non-negative. Hence, the Dirichlet form
$\mathcal{E}^{c}$ can be given by

$\displaystyle\mathcal{E}^{c}(f)=\int$
$\displaystyle\left(\tau_{1}\left\|\sqrt{L}\nabla_{x}f(x,y)\right\|^{2}+\tau_{2}\left\|\sqrt{L}\nabla_{y}f(x,y)\right\|^{2}\right)\mathrm{d}\mu\left(x,y\right)$
$\displaystyle+\int\frac{c}{2}s\left(x,y\right)\left(f(y,x)-f(x,y)\right)^{2}\mathrm{d}\mu\left(x,y\right)$

∎

### Proof of theorem 2.3

###### Proof.

The integral form of equation $(\ref{continue})$ is

$\boldsymbol{\xi}_{t}=\boldsymbol{\xi}_{0}-\int_{0}^{t}H\nabla
D(\boldsymbol{\xi}_{s}){\rm
d}s+\int_{0}^{t}\sqrt{H}\Sigma_{s}\mathrm{d}\boldsymbol{W}_{s},$

compare it with the equation $(\ref{continue-interpolation})$, we have

$\displaystyle\boldsymbol{\xi}_{t}-\boldsymbol{\xi}_{t}^{\delta}$
$\displaystyle=-H\int_{0}^{t}\left(\nabla D(\boldsymbol{\xi}_{t})-\nabla
D\left({\boldsymbol{\xi}}^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}\right)\right)\mathrm{d}s-(1-\eta)H\int_{0}^{t}\nabla
D\left({\boldsymbol{\xi}}^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}\right){\rm d}s$
$\displaystyle+\sqrt{H}\int_{0}^{t}\left(\Sigma_{s}-\Sigma^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}\right){\rm
d}\boldsymbol{W}_{s}+(1-\eta)\sqrt{H}\int_{0}^{t}\Sigma^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}{\rm d}\boldsymbol{W}_{s}.$

According to the compatibility of matrix norms, we apply the Cauchy-Schwartz
inequality and take expectation, then

$\displaystyle{\mathbb{E}}\left[\|\boldsymbol{\xi}_{t}-\boldsymbol{\xi}_{t}^{\delta}\|^{2}\right]$
$\displaystyle\leq
4\|H\|_{F}^{2}\underbrace{{\mathbb{E}}\left[\bigg{\|}\int_{0}^{t}\left(\nabla
D(\boldsymbol{\xi}_{s})-\nabla D\left({\boldsymbol{\xi}}^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}\right)\right)\mathrm{d}s\bigg{\|}^{2}\right]}_{{\cal
I}_{1}}$ (5.32)
$\displaystyle+4tr(H)\underbrace{{\mathbb{E}}\left[\bigg{\|}\int_{0}^{t}\left(\Sigma_{s}-\Sigma^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}\right){\rm
d}\boldsymbol{W}_{s}\bigg{\|}^{2}\right]}_{{\cal I}_{2}}$
$\displaystyle+4(1-\eta)^{2}\|H\|_{F}^{2}\underbrace{{\mathbb{E}}\left[\bigg{\|}\int_{0}^{t}\nabla
D\left({\boldsymbol{\xi}}^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}\right){\rm d}s\bigg{\|}^{2}\right]}_{{\cal J}_{1}}$
$\displaystyle+4(1-\eta)^{2}tr(H)\underbrace{{\mathbb{E}}\left[\bigg{\|}\int_{0}^{t}\Sigma^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}{\rm d}\boldsymbol{W}_{s}\bigg{\|}^{2}\right]}_{{\cal
J}_{2}}.$

where $tr(H)$ is the trace of the matrix $H$, the equality
$\|\sqrt{H}\|_{F}^{2}=tr(H)$ establishes.

For the term ${\cal J}_{1}$, by Cauchy-Schwarz inequality and the smooth
assumption, we have

$\displaystyle{\cal J}_{1}$
$\displaystyle={\mathbb{E}}\left[\bigg{\|}\int_{0}^{t}\left(\nabla
D\left({\boldsymbol{\xi}}^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}\right)-\nabla
D(\boldsymbol{\xi}^{*})\right)\mathrm{d}s\bigg{\|}^{2}\right]\leq
t{\mathbb{E}}\left[\int_{0}^{t}\bigg{\|}\left(\nabla
D\left({\boldsymbol{\xi}}^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}\right)-\nabla
D(\boldsymbol{\xi}^{*})\right)\bigg{\|}^{2}\mathrm{d}s\right]$
$\displaystyle\leq
tK^{2}{\mathbb{E}}\left[\int_{0}^{t}\left(\|{\boldsymbol{\xi}}^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}-\boldsymbol{\xi}^{*}\|^{2}\right)\mathrm{d}s\right]\leq
2tK^{2}{\mathbb{E}}\left[\int_{0}^{t}\left(\|{\boldsymbol{\xi}}^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}\|^{2}+\|\boldsymbol{\xi}^{*}\|^{2}\right)\mathrm{d}s\right]$
$\displaystyle\leq 2tK^{2}{\mathbb{E}}\left[\sum_{k=0}^{\lfloor
t/{\delta}\rfloor}\int_{k\delta}^{(k+1)\delta}\left(\|{\boldsymbol{\xi}}^{\delta}_{k\delta}\|^{2}+\|\boldsymbol{\xi}^{*}\|^{2}\right)\mathrm{d}s\right]$
$\displaystyle\leq 2tK^{2}(1+t/\delta)\delta\left(\sup_{k\geq
0}{\mathbb{E}}\left[\|{\boldsymbol{\xi}}^{\delta}_{k\delta}\|^{2}\right]+\|\boldsymbol{\xi}^{*}\|^{2}\right)$

where $\boldsymbol{\xi}^{*}$ is a stationary point of $D(\cdot)$, i.e.,
$\nabla D(\boldsymbol{\xi}^{*})=0$. Analogue to the Lemma C.2 in paper [5],
when $0<\frac{2\delta}{2+\delta}<\frac{\alpha_{1}}{\|B\|_{F}^{2}K^{2}}$, we
have

$\sup_{k\geq
0}{\mathbb{E}}\left[\|{\boldsymbol{\xi}}^{\delta}_{k\delta}\|^{2}\right]\leq\epsilon_{1}(n,\tau_{2},K,\alpha_{1},\alpha_{2},B)$
(5.33)

where the constant $\epsilon_{1}(n,\tau_{2},K,\alpha_{1},\alpha_{2},B)$
depends on the dimension $n$, the temperature parameter $\tau_{2}$, the
smoothness constant $K$, the constants $(\alpha_{1},\alpha_{2})$ in equation
$(\ref{B-dissipative})$ and the covariance matrix $B$. Hence, we have

${\cal J}_{1}\leq
2tK^{2}(1+t/\delta)\delta\epsilon_{2}(n,\tau_{2},K,\alpha_{1},\alpha_{2},B).$
(5.34)

For the term ${\cal J}_{2}$, according to the It$\hat{o}$ isometry, we have

$\displaystyle{\cal J}_{2}$
$\displaystyle=\sum_{j=1}^{2n}\int_{0}^{t}{\mathbb{E}}\left[\left(\Sigma^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}(j)\right)^{2}\right]{\rm d}s$
$\displaystyle\leq\sum_{j=1}^{2n}\sum_{k=0}^{\lfloor
t/\delta\rfloor}\int_{k\delta}^{(k+1)\delta}{\mathbb{E}}\left[\left(\Sigma^{\delta}_{k\delta}(j)\right)^{2}\right]{\rm
d}s$

where $\Sigma^{\delta}_{k\delta}(j)$ is the $j$-th diagonal entry of
$\Sigma^{\delta}_{k\delta}$. Recall that the values of
$\left(\Sigma^{\delta}_{k\delta}(j)\right)^{2}$ would be $2\tau_{1}$ or
$2\tau_{2}$, suppose the probability of being $2\tau_{2}$ is $p\in[0,1]$, then
the probability of being $2\tau_{1}$ is $1-p$, then

${\mathbb{E}}\left[\left(\Sigma^{\delta}_{k\delta}(j)\right)^{2}\right]=2p(\tau_{2}-\tau_{1})+2\tau_{1}\leq
2\tau_{2},$

hence we have

${\cal J}_{2}\leq 4(1+t/\delta)\delta\tau_{2}n.$ (5.35)

Similar to the analysis in paper [5], we have

$\displaystyle{\cal I}_{1}$ $\displaystyle\leq
t{\mathbb{E}}\left[\int_{0}^{t}\bigg{\|}\nabla D(\boldsymbol{\xi}_{s})-\nabla
D\left({\boldsymbol{\xi}}^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}\right)\bigg{\|}^{2}\mathrm{d}s\right]\leq
tK^{2}{\mathbb{E}}\left[\int_{0}^{t}\|\boldsymbol{\xi}_{s}-{\boldsymbol{\xi}}^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}\|^{2}\mathrm{d}s\right]$ $\displaystyle\leq
2K^{2}t\left({\mathbb{E}}\left[\int_{0}^{t}\|\boldsymbol{\xi}_{s}-{\boldsymbol{\xi}}^{\delta}_{s}\|^{2}\mathrm{d}s\right]+{\mathbb{E}}\left[\int_{0}^{t}\|\boldsymbol{\xi}^{\delta}_{s}-{\boldsymbol{\xi}}^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}\|^{2}\mathrm{d}s\right]\right)$

where

${\mathbb{E}}\left[\int_{0}^{t}\|\boldsymbol{\xi}^{\delta}_{s}-{\boldsymbol{\xi}}^{\delta}_{\lfloor
s/{\delta}\rfloor\delta}\|^{2}\mathrm{d}s\right]\leq(1+t/\delta)\left(4K^{2}\|H\|_{F}^{2}\delta^{3}\left(\sup_{k\geq
0}{\mathbb{E}}\left[\|{\boldsymbol{\xi}}^{\delta}_{k\delta}\|^{2}\right]+\|\boldsymbol{\xi}^{*}\|^{2}\right)+4tr(H)n\tau_{2}\delta^{2}\right),$

hence we have

${\cal I}_{1}\leq
2K^{2}t\left({\mathbb{E}}\left[\int_{0}^{t}\|\boldsymbol{\xi}_{s}-{\boldsymbol{\xi}}^{\delta}_{s}\|^{2}\mathrm{d}s\right]+\epsilon_{3}(n,\tau_{2},K,\alpha_{1},\alpha_{2},B)\delta\right).$
(5.36)

Since the term ${\cal I}_{2}$ is the same with the one in paper [5], we use
the result directly and have

${\cal I}_{2}\leq\epsilon_{4}(n,\tau_{1},\tau_{2},c)t\delta.$ (5.37)

By substituting the equations $(\ref{J1})$, $(\ref{J2})$, $(\ref{I1})$ and
$(\ref{I2})$ into $(\ref{total})$, we have

$\displaystyle{\mathbb{E}}\left[\|\boldsymbol{\xi}_{t}-\boldsymbol{\xi}_{t}^{\delta}\|^{2}\right]$
$\displaystyle\leq 8\|B\|_{F}^{2}\cdot
2K^{2}t\left({\mathbb{E}}\left[\int_{0}^{t}\|\boldsymbol{\xi}_{s}-{\boldsymbol{\xi}}^{\delta}_{s}\|^{2}\mathrm{d}s\right]+\epsilon_{3}(n,\tau_{2},K,\alpha_{1},\alpha_{2},B)\delta\right)$
$\displaystyle+8(1-\eta)^{2}(1+t/\delta)\|B\|_{F}^{2}\cdot
2tK^{2}\epsilon_{2}(n,\tau_{2},K,\alpha_{1},\alpha_{2},B)\delta$
$\displaystyle+8(1-\eta)^{2}(1+t/\delta)tr(B)\cdot
4\tau_{2}n\delta+8tr(B)\cdot\epsilon_{4}(n,\tau_{1},\tau_{2},c)t\delta,$

recall that $\eta=2/(2+\delta)$, $\delta\in(0,2)$, then for $t\in[0,T]$, we
have

$(1-\eta)^{2}(1+t/\delta)=\frac{\delta^{2}+t\delta}{(\delta+2)^{2}}\leq(\frac{1}{4}+\frac{T}{8}).$

Finally, we have

$\displaystyle{\mathbb{E}}\left[\|\boldsymbol{\xi}_{t}-\boldsymbol{\xi}_{t}^{\delta}\|^{2}\right]$
$\displaystyle\leq 8\|B\|_{F}^{2}\cdot
2K^{2}t\left({\mathbb{E}}\left[\int_{0}^{t}\|\boldsymbol{\xi}_{s}-{\boldsymbol{\xi}}^{\delta}_{s}\|^{2}\mathrm{d}s\right]+\epsilon_{3}(n,\tau_{2},K,\alpha_{1},\alpha_{2},B)\delta\right)$
$\displaystyle+\epsilon_{6}(n,\tau_{1},\tau_{2},c,K,\alpha_{1},\alpha_{2},T,B)\delta.$

Hence, by applying the Gr$\ddot{o}$nwall’s equality, there exists a constant
$\epsilon(n,\tau_{1},\tau_{2},c,K,\alpha_{1},\alpha_{2},T,B)$ such that for
$t\in[0,T]$,

${\mathbb{E}}\left[\|\boldsymbol{\xi}_{t}-\boldsymbol{\xi}_{t}^{\delta}\|^{2}\right]\leq\epsilon(n,\tau_{1},\tau_{2},c,K,\alpha_{1},\alpha_{2},T,B)\delta.$

∎