# Simulating Diffusion Bridges with Score Matching

Jeremy Heng 111Corresponding author<EMAIL_ADDRESS>ESSEC Business School,
Singapore 139408, Singapore Valentin De Bortoli Center for Science of Data,
ENS Ulm, Paris 75005, France Arnaud Doucet Department of Statistics,
University of Oxford, Oxford OX1 3LB, U.K. James Thornton44footnotemark: 4

###### Abstract

We consider the problem of simulating diffusion bridges, which are diffusion
processes that are conditioned to initialize and terminate at two given
states. The simulation of diffusion bridges has applications in diverse
scientific fields and plays a crucial role in the statistical inference of
discretely-observed diffusions. This is known to be a challenging problem that
has received much attention in the last two decades. This article contributes
to this rich body of literature by presenting a new avenue to obtain diffusion
bridge approximations. Our approach is based on a backward time representation
of a diffusion bridge, which may be simulated if one can time-reverse the
unconditioned diffusion. We introduce a variational formulation to learn this
time-reversal with function approximation and rely on a score matching method
to circumvent intractability. Another iteration of our proposed methodology
approximates the Doob’s $h$-transform defining the forward time representation
of a diffusion bridge. We discuss algorithmic considerations and extensions,
and present numerical results on an Ornstein–Uhlenbeck process, a model from
financial econometrics for interest rates, and a model from genetics for cell
differentiation and development to illustrate the effectiveness of our
approach.

Keywords: Diffusion; Diffusion bridge; Score matching; Stochastic differential
equation; Time-reversal.

##  1 Introduction

Diffusion processes have been used extensively in mathematical and natural
sciences. A diffusion process $X=(X_{t})_{t\in[0,T]}$ in $\mathbb{R}^{d}$ is
defined by the stochastic differential equation

$\displaystyle\mathrm{d}X_{t}=f(t,X_{t})\mathrm{d}t+\sigma(t,X_{t})\mathrm{d}W_{t},$
(1)

where $f:[0,T]\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is a drift
function, $\sigma:[0,T]\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times d}$
is a diffusion coefficient, and $W=(W_{t})_{t\in[0,T]}$ is a $d$-dimensional
Brownian motion. We suppose $f$ and $\sigma$ are sufficiently regular to
induce a unique weak solution and
$\Sigma(t,x_{t})=(\sigma\sigma)^{\mathrm{\scriptscriptstyle T}}(t,x_{t})$ is
uniformly positive definite for all $(t,x_{t})\in[0,T]\times\mathbb{R}$. For
any $0\leq s<t\leq T$, we denote the transition density of (1) with respect to
the Lebesgue measure on $\mathbb{R}^{d}$ as $p(t,x_{t}\mid s,x_{s})$ and
assume that it is positive for ease of exposition. While the numerical
simulation of $X$ can be routinely handled by time-discretization schemes
(Kloeden and Platen, 1992), the task of simulating $X$ initialized at
$X_{0}=x_{0}$ and conditioned to terminate at $X_{T}=x_{T}$ is a challenging
problem that has received much attention in the last two decades.

Simulating the conditioned process $X^{\star}=(X_{t}^{\star})_{t\in[0,T]}$,
commonly referred to as a diffusion bridge, has applications in diverse fields
such as computational chemistry (Bolhuis et al., 2002; Wang et al., 2020),
financial econometrics (Elerian et al., 2001; Durham and Gallant, 2002),
genetics (Wang et al., 2011), and shape analysis (Arnaudon et al., 2022). When
performing statistical inference for parameters of $f$ and $\sigma$ in the
case where $X$ is observed at discrete time points, diffusion bridge
simulation is a crucial tool that allows one to impute missing paths between
observations within an expectation-maximization algorithm or a Gibbs sampler
(Pedersen, 1995; Roberts and Stramer, 2001; Eraker, 2001; Beskos et al., 2006;
Golightly and Wilkinson, 2008; van der Meulen and Schauer, 2017).

By Doob’s $h$-transform (Rogers and Williams, 2000, p. 83), it is well-known
that $X^{\star}$ satisfies

$\displaystyle\mathrm{d}X_{t}^{\star}=\\{f(t,X_{t}^{\star})+\Sigma(t,X_{t}^{\star})\nabla\log
h(t,X_{t}^{\star})\\}\mathrm{d}t+\sigma(t,X_{t}^{\star})\mathrm{d}W_{t},\quad
X_{0}^{\star}=x_{0},$ (2)

where $h(t,x_{t})=p(T,x_{T}\mid t,x_{t})$ and $\nabla$ denotes the gradient
operator. The term $\Sigma(t,x_{t})\nabla\log h(t,x_{t})$ forces the
conditioned process towards the terminal condition $X_{T}^{\star}=x_{T}$. As
the transition density and hence its logarithmic gradient is intractable for
most diffusions, exploiting this result to simulate diffusion bridges is
highly non-trivial. To this end, one can characterize $h$ as the solution of
the backward Kolmogorov equation

$\displaystyle\partial_{t}h(t,x_{t})+(\mathcal{L}h)(t,x_{t})=0,$ (3)

with terminal condition at time $T$ given by the Dirac measure at $x_{T}$,
where $\mathcal{L}$ denotes the generator of $X$. Equation (3) reveals that
$h$ propagates information about the terminal constraint backwards in time.
However, numerical resolution of this partial differential equation is
particularly challenging due to the singularity at time $T$, and
computationally demanding when the dimension $d$ is large (Wang et al., 2020).
Furthermore, one must run a solver for every pair of conditioned states
$(x_{0},x_{T})$ considered.

A common approach to address these difficulties is to simulate a proposal
bridge process $X^{\circ}=(X_{t}^{\circ})_{t\in[0,T]}$, satisfying
$\mathrm{d}X_{t}^{\circ}=f^{\circ}(t,X_{t}^{\circ})\mathrm{d}t+\sigma(t,X_{t}^{\circ})\mathrm{d}W_{t}$
with $X_{0}^{\circ}=x_{0}$. One constructs
$f^{\circ}:[0,T]\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ using a
tractable approximation of (2), and corrects for the discrepancy using
importance sampling or an independent Metropolis–Hastings algorithm
(Papaspiliopoulos and Roberts, 2012; Elerian et al., 2001). The simple choice
$f^{\circ}(t,x_{t})=f(t,x_{t})$ typically performs poorly as it does not take
the constraint $X_{T}=x_{T}$ into account (Pedersen, 1995). The drift of a
Brownian bridge $f^{\circ}(t,x_{t})=(x_{T}-x_{t})/(T-t)$ has been considered
in several works (Durham and Gallant, 2002; Delyon and Hu, 2006; Stramer and
Yan, 2007; Papaspiliopoulos et al., 2013), and improved by Whitaker et al.
(2017) using an innovative decomposition of the process into deterministic and
stochastic parts. Clark (1990) followed by Delyon and Hu (2006) studied the
choice $f^{\circ}(t,x_{t})=f(t,x_{t})+(x_{T}-x_{t})/(T-t)$ that incorporates
the dynamics of the original process $X$. To introduce more flexibility and
better mimic the structure of (2), Schauer et al. (2017) proposed setting
$f^{\circ}(t,x_{t})=f(t,x_{t})+\Sigma(t,x_{t})\nabla\log\tilde{h}(t,x_{t})$,
where $\tilde{h}(t,x_{t})=\tilde{p}(T,x_{T}\mid t,x_{t})$ is an analytically
tractable transition density of an auxiliary process. For tractability, the
latter is typically chosen from the class of linear processes and can be
optimized to get the best approximation within this class (van der Meulen and
Schauer, 2017). Other Markov chain Monte Carlo approaches include Gibbs
sampling (Eraker, 2001), Langevin-type stochastic partial differential
equations (Stuart et al., 2004; Beskos et al., 2008), and piecewise
deterministic Monte Carlo (Bierkens et al., 2021).

The exact simulation algorithms developed in Beskos and Roberts (2005) and
Beskos et al. (2006) can be employed to sample diffusion bridges without any
time-discretization error. However, these elegant methods are limited to the
class of diffusion processes that can be transformed to have unit diffusion
coefficient. Bladt and Sørensen (2014) and Bladt et al. (2016) devised an
ingenious methodology to simulate diffusion bridges based on coupling and
time-reversal of diffusions. Their proposed method is applicable to the class
of ergodic diffusions with an invariant density that is either explicitly
known or numerically approximated. Closely related approaches include
sequential Monte Carlo algorithms that resample using backward information
filter approximations (Guarniero, 2017), information from backward pilot paths
(Lin et al., 2010), or guided weight functions (Del Moral and Murray, 2015).
The main idea underlying these works is the representation of the diffusion
bridge in (2) and (3).

This work contributes to this rich body of literature by presenting a new
avenue to obtain diffusion bridge approximations. In Section 2.1, we begin by
examining the dynamics of the time-reversed bridge process $X^{\star}$, which
amounts to a time-reversal of (1). Section 2.2 introduces a variational
formulation to learn this time-reversal using function approximation and a
score matching method (Vincent, 2011). In Section 2.3, we explain how another
iteration of this methodology approximates the Doob’s $h$-transformed process
in (2). Our approach is widely applicable and can be used to construct
proposal bridge processes within existing frameworks. Lastly, in Section 3 we
illustrate computational improvements over common choices of proposal bridge
processes on three numerical examples. A Python package is available on
GitHub.

##  2 Diffusion bridges

###  2.1 Time-reversed bridge process

It can be shown that the time-reversed bridge process
$Z^{\star}=(Z_{t}^{\star})_{t\in[0,T]}=(X_{T-t}^{\star})_{t\in[0,T]}$
satisfies

$\displaystyle\mathrm{d}Z_{t}^{\star}=b(t,Z_{t}^{\star})\mathrm{d}t+\sigma(T-t,Z_{t}^{\star})\mathrm{d}B_{t},\quad
Z_{0}^{\star}=x_{T},$ (4)

with drift function
$b(t,z_{t})=-f(T-t,z_{t})+\Sigma(T-t,z_{t})s(T-t,z_{t})+\nabla\cdot\Sigma(T-t,z_{t})$,
another standard Brownian motion $B=(B_{t})_{t\in[0,T]}$,
$s(t,x_{t})=s^{\star}(t,x_{t})-\nabla\log h(t,x_{t})$, and
$\nabla\cdot\Sigma(t,x)=(\sum_{j=1}^{d}\partial_{x_{j}}\Sigma^{1,j}(t,x),\ldots,\sum_{j=1}^{d}\partial_{x_{j}}\Sigma^{d,j}(t,x))$
is the divergence of $\Sigma$. Here $s^{\star}(t,x_{t})=\nabla\log
p^{\star}(t,x_{t})$ denotes the score of the marginal density
$p^{\star}(t,x_{t})$ of the diffusion bridge process $X_{t}^{\star}$ at time
$t\in(0,T)$. We refer readers to Haussmann and Pardoux (1986) and Millet et
al. (1989) for conditions under which the representation in (4) holds. By the
Markov property, we have the relation

$\displaystyle
p^{\star}(t,x_{t})=p\\{t,x_{t}\mid(0,x_{0}),(T,x_{T})\\}=\frac{p(t,x_{t}\mid
0,x_{0})h(t,x_{t})}{p(T,x_{T}\mid 0,x_{0})},$ (5)

as $h(t,x_{t})=p(T,x_{T}\mid t,x_{t})$. This implies that
$s(t,x_{t})=\nabla\log p(t,x_{t}\mid 0,x_{0})$ is simply the score of the
transition density of $X$.

Exploiting this backward time representation to derive diffusion bridge
approximations is also highly non-trivial due to the intractability of the
transition density $\eta(t,x_{t})=p(t,x_{t}\mid 0,x_{0})$, which is now
characterized by the forward Kolmogorov equation
$\partial_{t}\eta(t,x_{t})+(\mathcal{F}\eta)(t,x_{t})=0$, with initial
condition at time $0$ given by the Dirac measure at $x_{0}$ and $\mathcal{F}$
denotes the Fokker–Planck operator of $X$. Numerical resolution of $\eta$
using partial differential equation solvers also suffers from the same
difficulties as (3). A key observation is that (4) can be understood as first
setting $Z_{0}^{\star}=x_{T}$ to satisfy the terminal constraint, and then
evolving $Z^{\star}$ using the time-reversal of (1). Due to the influence of
the score $s(t,x_{t})$, the process will end at the initial constraint
$Z_{T}^{\star}=x_{0}$ by construction. This connection between simulation of a
diffusion bridge and time-reversal of its original diffusion process will form
the basis of our score approximation.

###  2.2 Learning time-reversal with score matching

We introduce a variational formulation to learn the time-reversal of (1),
involving path measures on the space of continuous functions from $[0,T]$ to
$\mathbb{R}^{d}$, equipped with the cylinder $\sigma$-algebra. Let
$\mathbb{P}^{x_{0}}$ denote the path measure induced by (1) with
$X_{0}=x_{0}$, and $\mathbb{Q}_{\phi}^{x_{0}}$ be the path measure induced by
a time-reversed process $Z=(Z_{t})_{t\in[0,T]}$ satisfying

$\displaystyle\mathrm{d}Z_{t}=b_{\phi}(t,Z_{t})\mathrm{d}t+\sigma(T-t,Z_{t})\mathrm{d}B_{t},\quad
Z_{0}\sim p(T,dx_{T}\mid 0,x_{0}),$ (6)

with drift function
$b_{\phi}(t,z_{t})=-f(T-t,z_{t})+\Sigma(T-t,z_{t})s_{\phi}(T-t,z_{t})+\nabla\cdot\Sigma(T-t,z_{t})$
that mimics the form of $b(t,z_{t})$ in (4). Here
$s_{\phi}:[0,T]\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ represents a
function approximation of the score $s(t,x_{t})$ that depends on parameters
$\phi\in\Phi$ to be optimized. We shall measure the score approximation error
as

$\displaystyle
e(\phi)={E}^{x_{0}}\left\\{\int_{0}^{T}\left\|s_{\phi}(t,X_{t})-s(t,X_{t})\right\|_{\Sigma(t,X_{t})}^{2}\mathrm{d}t\right\\},$
(7)

where $E^{x_{0}}$ denotes expectation with respect to $\mathbb{P}^{x_{0}}$ and
$\|\cdot\|_{A}$ denotes the Euclidean norm weighted by a positive definite
$A\in\mathbb{R}^{d\times d}$.

Let $\mathrm{d}\mathbb{P}^{x_{0}}/\mathrm{d}\mathbb{Q}_{\phi}^{x_{0}}$ denote
the Radon–Nikodym derivative of $\mathbb{P}^{x_{0}}$ with respect
$\mathbb{Q}_{\phi}^{x_{0}}$. We consider minimizing the Kullback–Leibler
divergence
$\textsc{kl}(\mathbb{P}^{x_{0}}|\mathbb{Q}_{\phi}^{x_{0}})=E^{x_{0}}\\{\log\mathrm{d}\mathbb{P}^{x_{0}}/\mathrm{d}\mathbb{Q}_{\phi}^{x_{0}}(X)\\}$.
The following result gives an expression of this objective and shows that the
process $Z$ will end at the initial constraint $Z_{T}=x_{0}$ if the score
approximation error is finite.

######  Proposition 1.

Assuming $e(\phi)<\infty$, we have
$\textsc{kl}(\mathbb{P}^{x_{0}}|\mathbb{Q}_{\phi}^{x_{0}})=e(\phi)/2$ and
$Z_{T}=x_{0}$ holds $\mathbb{Q}_{\phi}^{x_{0}}$-almost surely.

Although this clearly relates learning the time-reversal of (1) with
approximation of the score $s(t,x_{t})$, its form is not amenable to
optimization. The following result gives an alternative and practical
expression by adapting the idea of denoising score matching (Vincent, 2011) to
our setting.

######  Proposition 2.

For any partition $(t_{m})_{m=0}^{M}$ of the interval $[0,T]$, we have
$\textsc{kl}(\mathbb{P}^{x_{0}}|\mathbb{Q}_{\phi}^{x_{0}})=L(\phi)+C$ if
$e(\phi)<\infty$, where $C$ is a constant independent of $\phi\in\Phi$, the
loss function is defined as

$\displaystyle
L(\phi)=\frac{1}{2}\sum_{m=1}^{M}\int_{t_{m-1}}^{t_{m}}{E}^{x_{0}}\left\\{\left\|s_{\phi}(t,X_{t})-g(t_{m-1},X_{t_{m-1}},t,X_{t})\right\|_{\Sigma(t,X_{t})}^{2}\right\\}\mathrm{d}t,$
(8)

and $g(s,x_{s},t,x_{t})=\nabla\log p(t,x_{t}\mid s,x_{s})$ for $0\leq s<t\leq
T$ and $x_{s},x_{t}\in\mathbb{R}^{d}$.

Therefore minimizing the Kullback–Leibler divergence
$\textsc{kl}(\mathbb{P}^{x_{0}}|\mathbb{Q}_{\phi}^{x_{0}})$ is equivalent to
minimizing the loss function $L(\phi)$. This allows us to circumvent the
intractable score $s(t,x_{t})$ by working with
$g(t_{m-1},x_{t_{m-1}},t,x_{t})$, the score of the transition density
$p(t,x_{t}\mid t_{m-1},x_{t_{m-1}})$. Although the latter is also intractable,
approximations can be made when the sub-interval $[t_{m-1},t_{m}]$ is
sufficiently small. For example, under the Euler–Maruyama scheme (Kloeden and
Platen, 1992, p. 340) with stepsize $\delta t=T/M$,

$\displaystyle g(t_{m-1},x_{t_{m-1}},t_{m},x_{t_{m}})\approx-(\delta
t)^{-1}\Sigma^{-1}(t_{m-1},x_{t_{m-1}})\\{x_{t_{m}}-x_{t_{m-1}}-\delta
tf(t_{m-1},x_{t_{m-1}})\\}.$

Hence the loss $\mathrm{L}(\phi)$ can be approximated and minimized using
stochastic gradient algorithms by simulating time-discretized paths under (1).
The minimal loss of $\mathrm{L}(\phi)=-C$, achieved when
$s_{\phi}(t,x_{t})=s(t,x_{t})$ $\mathbb{P}^{x_{0}}$-almost surely, is unknown
as the constant $C$ is typically incomputable. After obtaining the score
approximation $s_{\phi}$, we can simulate a proposal bridge $Z$ from (6) with
$Z_{0}=x_{T}$ and correct it using importance sampling or independent
Metropolis–Hastings. Time-discretization considerations and proposal
correction procedures are detailed in Appendix C.

In scenarios where one is interested in multiple pairs of conditioned states
$(x_{0},x_{T})$, we can extend the above methodology to avoid having to learn
multiple score approximations as follows. We let the score approximation
$s_{\phi}(t,x_{t})$ in (6) also depend on the initial condition $x_{0}$, and
average the Kullback–Leibler objective
$\textsc{kl}(\mathbb{P}^{x_{0}}|\mathbb{Q}_{\phi}^{x_{0}})$ with a
distribution $p_{0}(\mathrm{d}x_{0})$ on $\mathbb{R}^{d}$ that can be sampled
from. By applying the arguments of Proposition 2 conditionally on
$X_{0}=x_{0}$, we obtain a loss function given by averaging (8) over
$p_{0}(\mathrm{d}x_{0})$, which can be minimized using time-discretization and
stochastic gradient algorithms.

###  2.3 Learning Doob’s $h$-transform

It is instructive to consider the time-reversal of (4), which gives

$\displaystyle\mathrm{d}X_{t}^{\star}=f^{\star}(t,X_{t}^{\star})\mathrm{d}t+\sigma(t,X_{t}^{\star})\mathrm{d}W_{t},\quad
X_{0}^{\star}=x_{0},$ (9)

with drift function
$f^{\star}(t,x_{t})=-b(T-t,x_{t})+\Sigma(t,x_{t})s^{\star}(t,x_{t})+\nabla\cdot\Sigma(t,x_{t})$.
Using the form of $b(t,z_{t})$ and the relation
$s(t,x_{t})=s^{\star}(t,x_{t})-\nabla\log h(t,x_{t})$, we can rewrite
$f^{\star}$ as

$\displaystyle
f^{\star}(t,x_{t})=f(t,x_{t})+\Sigma(t,x_{t})\\{s^{\star}(t,x_{t})-s(t,x_{t})\\}=f(t,x_{t})+\Sigma(t,x_{t})\nabla\log
h(t,x_{t}),$ (10)

therefore (9) recovers the Doob’s $h$-transform in (2). Although this is to be
expected as the reversal of the time-reversed process should recover the
original process, it forms the basis of our approximation of (2).

After obtaining an approximation of $s(t,x_{t})$, another iteration of our
proposed methodology can be used to learn a function approximation of
$s^{\star}(t,x_{t})$. For brevity, this is detailed in Appendix B. Plugging in
both approximations in (10) then gives an approximation of $X^{\star}$, which
could be of interest in algorithms where one requires a forward time
representation of the proposal bridge process (Lin et al., 2010; Del Moral and
Murray, 2015).

##  3 Numerical examples

###  3.1 Preliminaries

As our methodology allows one to employ any function approximator, we harness
the flexibility of neural networks and the ease of implementation using modern
software to approximate score functions. The choice of neural network and
stochastic gradient algorithm is detailed in Appendix D. Optimization times
ranged between several seconds to a few minutes on a simple desktop machine
and can be reduced with hardware accelerators. As such computational overheads
are marginal when deploying proposal bridge processes within an importance
sampler with many samples or an independent Metropolis–Hastings algorithm with
many iterations, we focus on assessing the quality of our proposals in
settings where existing proposal methods are unsatisfactory. The performance
measures considered are the importance sampling effective sample size
proportion, the acceptance rate of independent Metropolis–Hastings, and the
variance of the importance sampling log-transition density estimator or its
mean squared error when the true transition density is known. These measures
were computed using $1024$ samples or iterations, and $100$ independent
repetitions of each method. We benchmark our approximations of the backward
and forward diffusion bridges in (4) and (9) against the forward diffusion
method of Pedersen (1995), the modified diffusion bridge of Durham and Gallant
(2002), and the Clark–Delyon–Hu proposal bridge process studied by Clark
(1990) and Delyon and Hu (2006). Given the difficulty of comparing the wide
range of methods for diffusion bridges in a completely fair manner, as their
strengths and weaknesses can depend on the specificities of the problem under
consideration, we note that our objective is merely to illustrate a new avenue
to improve the construction of proposal bridge processes.

###  3.2 Ornstein–Uhlenbeck process

Consider (1) with linear drift function $f(t,x_{t})=-2x_{t}$, identity
diffusion coefficient $\sigma(t,x_{t})=I_{d}$, and the constraints
$x_{0}=x_{T}=(1,\ldots,1)$. The transition density and score function are
explicitly known and used as ground truth. Fig. 1 illustrates the impact of
the time horizon $T$ and dimension $d$ on algorithmic performance. We observe
that our proposal bridge processes offer substantial improvements over
existing methods in all settings considered. When comparing our forward and
backward diffusion bridge approximations, we notice some accumulation of
error, which is to be expected as the forward process is constructed using an
additional score approximation. Given this observation, we will only consider
the backward process in the following.

Figure 1: Results for Ornstein–Uhlenbeck example with fixed $d=1$ (_first
row_) and fixed $T=1$ (_second row_). The symbols and lines correspond to our
approximations of the backward (_cross-dotted_) and forward (_plus-dotted_)
diffusion bridges, the forward diffusion (_circle-solid_), the modified
diffusion bridge (_triangle-dashed_), and the Clark–Delyon–Hu bridge (_square-
dash-dotted_).

###  3.3 Interest rates model

Next, we consider a special case of an interest rates model in Aït-Sahalia and
Lo (1998), defined by (1) with $f(t,x_{t})=4/x_{t}-x_{t}$ and
$\sigma(t,x_{t})=1$. This specification admits a tractable transition density
and score function as ground truth, given in terms of modified Bessel
functions and its derivative. For each $T$, we first learn a single score
approximation to handle multiple conditioned states $(x_{0},x_{T})$ by
minimizing the loss in (8) averaged over initial states $X_{0}=x_{0}$ from the
gamma distribution with shape $5$ and rate $2$. We visualize our score
approximation for $T=1$ in Fig. 2, and examine how algorithmic performance
depends on $T$ and $(x_{0},x_{T})$ in Fig. 3. While the performance of
modified diffusion bridge typically degrades with $T$, the forward diffusion
and Clark–Delyon–Hu bridge may yield more satisfactory results for specific
conditioned states (second column). In contrast, our backward diffusion bridge
approximation performs well for all $T$ considered and when conditioning
requires the process to move away from the stable stationary point of $f$ at
$x=2$ (first and third columns).

Figure 2: Comparing true score function $s(t,x)$ (_first row_) and our
function approximation $s_{\phi}(t,x)$ (_second row_) for various time $t$ and
initial condition $x_{0}$. Figure 3: Results for interest rates model for
three pairs of conditioned states $(x_{0},x_{T})$. The symbols and lines
correspond to our approximation of the backward diffusion bridge (_cross-
dotted_), the forward diffusion (_circle-solid_), the modified diffusion
bridge (_triangle-dashed_), and the Clark–Delyon–Hu bridge (_square-dash-
dotted_).

###  3.4 Cell model

We end with a cell differentiation and development model from Wang et al.
(2011). Cellular expression levels $X_{t}=(X_{t,1},X_{t,2})$ of two genes are
modelled by (1) with $f=(f_{1},f_{2})$, where
$f_{i}(t,x_{t})=x_{t,i}^{4}/(2^{-4}+x_{t,i}^{4})+2^{-4}/(2^{-4}+x_{t,j}^{4})-x_{t,i}$,
for $(i,j)\in\\{(1,2),(2,1)\\}$, describe self-activation, mutual inhibition
and inactivation respectively, and $\sigma(t,x_{t})=\sigma_{X}I_{2}$ captures
intrinsic or external fluctuations. We consider the cell development from the
undifferentiated state of $x_{0}=(1,1)$ to a differentiated state
$x_{T}=(x_{T,1},x_{T,2})$ defined by another stable fixed point of $f$
satisfying $x_{T,1}<x_{T,2}$. Fig. 4 displays the performance of all proposal
methods to approximate this diffusion bridge. For this application, we observe
significant improvement over other methods when the time horizon is long and
when the diffusion coefficient $\sigma_{X}$ is smaller.

Figure 4: Results for cell model with $\sigma_{X}^{2}=0.1$ (_first row_) and
$\sigma_{X}^{2}=1$ (_second row_). The symbols and lines correspond to our
approximation of the backward diffusion bridge (_cross-dotted_), the forward
diffusion (_circle-solid_), the modified diffusion bridge (_triangle-dashed_),
and the Clark–Delyon–Hu bridge (_square-dash-dotted_).

## Acknowledgements

Jeremy Heng was funded by CY Initiative of Excellence (grant “Investissements
d’Avenir” ANR-16-IDEX-0008). Valentin De Bortoli and Arnaud Doucet are partly
supported by the EPSRC grant CoSInES EP/R034710/1.

## References

  * Aït-Sahalia and Lo (1998) Aït-Sahalia, Y. and A. W. Lo (1998).  Nonparametric estimation of state-price densities implicit in financial asset prices.  The Journal of Finance 53(2), 499–547.
  * Andrieu et al. (2010) Andrieu, C., A. Doucet, and R. Holenstein (2010).  Particle Markov chain Monte Carlo methods (with discussion).  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72(3), 269–342.
  * Arnaudon et al. (2022) Arnaudon, A., F. van der Meulen, M. Schauer, and S. Sommer (2022).  Diffusion bridges for stochastic Hamiltonian systems and shape evolutions.  SIAM Journal on Imaging Sciences 15(1), 293–323.
  * Beskos et al. (2006) Beskos, A., O. Papaspiliopoulos, and G. O. Roberts (2006).  Retrospective exact simulation of diffusion sample paths with applications.  Bernoulli 12(6), 1077–1098.
  * Beskos et al. (2006) Beskos, A., O. Papaspiliopoulos, G. O. Roberts, and P. Fearnhead (2006).  Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes (with discussion).  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68(3), 333–382.
  * Beskos and Roberts (2005) Beskos, A. and G. O. Roberts (2005).  Exact simulation of diffusions.  The Annals of Applied Probability 15(4), 2422–2444.
  * Beskos et al. (2008) Beskos, A., G. O. Roberts, A. Stuart, and J. Voss (2008).  MCMC methods for diffusion bridges.  Stochastics and Dynamics 8(03), 319–350.
  * Bierkens et al. (2021) Bierkens, J., S. Grazzi, F. van der Meulen, and M. Schauer (2021).  A piecewise deterministic Monte Carlo method for diffusion bridges.  Statistics and Computing 31(3), 1–21.
  * Bladt et al. (2016) Bladt, M., S. Finch, and M. Sørensen (2016).  Simulation of multivariate diffusion bridges.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 343–369.
  * Bladt and Sørensen (2014) Bladt, M. and M. Sørensen (2014).  Simple simulation of diffusion bridges with application to likelihood inference for diffusions.  Bernoulli 20(2), 645–675.
  * Bolhuis et al. (2002) Bolhuis, P. G., D. Chandler, C. Dellago, and P. L. Geissler (2002).  Transition path sampling: Throwing ropes over rough mountain passes, in the dark.  Annual Review of Physical Chemistry 53(1), 291–318.
  * Clark (1990) Clark, J. M. C. (1990).  The simulation of pinned diffusions.  In 29th IEEE Conference on Decision and Control, pp. 1418–1420. IEEE.
  * Del Moral and Murray (2015) Del Moral, P. and L. M. Murray (2015).  Sequential Monte Carlo with highly informative observations.  SIAM/ASA Journal on Uncertainty Quantification 3(1), 969–997.
  * Delyon and Hu (2006) Delyon, B. and Y. Hu (2006).  Simulation of conditioned diffusion and application to parameter estimation.  Stochastic Processes and their Applications 116(11), 1660–1675.
  * Durham and Gallant (2002) Durham, G. B. and A. R. Gallant (2002).  Numerical techniques for maximum likelihood estimation of continuous-time diffusion processes.  Journal of Business & Economic Statistics 20(3), 297–338.
  * Elerian et al. (2001) Elerian, O., S. Chib, and N. Shephard (2001).  Likelihood inference for discretely observed nonlinear diffusions.  Econometrica 69(4), 959–993.
  * Eraker (2001) Eraker, B. (2001).  MCMC analysis of diffusion models with application to finance.  Journal of Business & Economic Statistics 19(2), 177–191.
  * Golightly and Wilkinson (2008) Golightly, A. and D. J. Wilkinson (2008).  Bayesian inference for nonlinear multivariate diffusion models observed with error.  Computational Statistics & Data Analysis 52(3), 1674–1693.
  * Guarniero (2017) Guarniero, P. (2017).  The Iterated Auxiliary Particle Filter and Applications to State Space Models and Diffusion Processes.  Ph. D. thesis, University of Warwick.
  * Haussmann and Pardoux (1986) Haussmann, U. G. and E. Pardoux (1986).  Time reversal of diffusions.  The Annals of Probability, 1188–1205.
  * Kingma and Ba (2014) Kingma, D. P. and J. Ba (2014).  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Kloeden and Platen (1992) Kloeden, P. E. and E. Platen (1992).  Stochastic differential equations.  In Numerical Solution of Stochastic Differential Equations, pp. 103–160. Springer.
  * Lin et al. (2010) Lin, M., R. Chen, and P. Mykland (2010).  On generating Monte Carlo samples of continuous diffusion bridges.  Journal of the American Statistical Association 105(490), 820–838.
  * Middleton et al. (2019) Middleton, L., G. Deligiannidis, A. Doucet, and P. E. Jacob (2019).  Unbiased smoothing using particle independent Metropolis–Hastings.  In International Conference on Artificial Intelligence and Statistics. PMLR.
  * Millet et al. (1989) Millet, A., D. Nualart, and M. Sanz (1989).  Integration by parts and time reversal for diffusion processes.  The Annals of Probability, 208–238.
  * Papaspiliopoulos and Roberts (2012) Papaspiliopoulos, O. and G. O. Roberts (2012).  Importance sampling techniques for estimation of diffusion models.  In M. Kessler, A. Lindner, and M. Sorensen (Eds.), Statistical Methods for Stochastic Differential Equations, pp. 311–340. Chapman and Hall/CRC.
  * Papaspiliopoulos et al. (2013) Papaspiliopoulos, O., G. O. Roberts, and O. Stramer (2013).  Data augmentation for diffusions.  Journal of Computational and Graphical Statistics 22(3), 665–688.
  * Pedersen (1995) Pedersen, A. R. (1995).  Consistency and asymptotic normality of an approximate maximum likelihood estimator for discretely observed diffusion processes.  Bernoulli, 257–279.
  * Roberts and Stramer (2001) Roberts, G. O. and O. Stramer (2001).  On inference for partially observed nonlinear diffusion models using the Metropolis–Hastings algorithm.  Biometrika 88(3), 603–621.
  * Rogers and Williams (2000) Rogers, L. C. G. and D. Williams (2000).  Diffusions, Markov processes and Martingales: Volume 2: Itô Calculus, Volume 2.  Cambridge university press.
  * Schauer et al. (2017) Schauer, M., F. Van Der Meulen, and H. Van Zanten (2017).  Guided proposals for simulating multi-dimensional diffusion bridges.  Bernoulli 23(4A), 2917–2950.
  * Stramer and Yan (2007) Stramer, O. and J. Yan (2007).  On simulated likelihood of discretely observed diffusion processes and comparison to closed-form approximation.  Journal of Computational and Graphical Statistics 16(3), 672–691.
  * Stuart et al. (2004) Stuart, A. M., J. Voss, and P. Wilberg (2004).  Conditional path sampling of SDEs and the Langevin MCMC method.  Communications in Mathematical Sciences 2(4), 685–697.
  * van der Meulen and Schauer (2017) van der Meulen, F. and M. Schauer (2017).  Bayesian estimation of discretely observed multi-dimensional diffusion processes using guided proposals.  Electronic Journal of Statistics 11(1), 2358–2396.
  * Vaswani et al. (2017) Vaswani, A., N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin (2017).  Attention is all you need.  In Advances in Neural Information Processing Systems, pp. 5998–6008.
  * Vincent (2011) Vincent, P. (2011).  A connection between score matching and denoising autoencoders.  Neural Computation 23(7), 1661–1674.
  * Wang et al. (2011) Wang, J., K. Zhang, L. Xu, and E. Wang (2011).  Quantifying the Waddington landscape and biological paths for development and differentiation.  Proceedings of the National Academy of Sciences 108(20), 8257–8262.
  * Wang et al. (2020) Wang, S., D. Ramkrishna, and V. Narsimhan (2020).  Exact sampling of polymer conformations using Brownian bridges.  The Journal of Chemical Physics 153(3), 034901.
  * Whitaker et al. (2017) Whitaker, G. A., A. Golightly, R. J. Boys, and C. Sherlock (2017).  Improved bridge constructs for stochastic differential equations.  Statistics and Computing 27(4), 885–900.

##  Appendix A Proofs of Propositions 1 and 2

In the following, we write $\langle
x,y\rangle_{A}=x^{\mathrm{\scriptscriptstyle T}}Ay$ for the inner product of
$x,y\in\mathbb{R}^{d}$ weighted by a positive definite matrix
$A\in\mathbb{R}^{d\times d}$ and $\|x\|_{A}=\sqrt{\langle x,x\rangle_{A}}$ for
the induced weighted Euclidean norm.

###### Proof of Proposition 1.

It follows from Section 2.1 that the path measure $\mathbb{P}^{x_{0}}$ is also
induced by the backward process $Z=(Z_{t})_{t\in[0,T]}=(X_{T-t})_{t\in[0,T]}$
satisfying

$\displaystyle\mathrm{d}Z_{t}=b(t,Z_{t})\mathrm{d}t+\sigma(T-t,Z_{t})\mathrm{d}B_{t},\quad
Z_{0}\sim p(T,dx_{T}\mid 0,x_{0}).$ (11)

Using Girsanov’s theorem (Rogers and Williams, 2000, p. 79) and Equations (11)
and (6), the Radon–Nikodym derivative of $\mathbb{P}^{x_{0}}$ with respect
$\mathbb{Q}_{\phi}^{x_{0}}$ is

$\displaystyle\frac{\mathrm{d}\mathbb{P}^{x_{0}}}{\mathrm{d}\mathbb{Q}_{\phi}}(Z)=\exp\left\\{-\frac{1}{2}\int_{0}^{T}\|(b-b_{\phi})(t,Z_{t})\|_{\Sigma^{-1}(T-t,Z_{t})}^{2}\mathrm{d}t+\int_{0}^{T}(b-b_{\phi})(t,Z_{t})^{\mathrm{\scriptscriptstyle
T}}(\Sigma^{-1}\sigma)(T-t,Z_{t})\mathrm{d}B_{t}\right\\}.$

By properties of Itô integrals, the form of $b(t,z_{t})$ and
$b_{\phi}(t,z_{t})$, and the identity $\|Ax\|_{A^{-1}}^{2}=\|x\|_{A}^{2}$ for
a positive definite matrix $A\in\mathbb{R}^{d\times d}$ and
$x\in\mathbb{R}^{d}$, we can write the Kullback–Leibler divergence as

$\displaystyle\textsc{kl}(\mathbb{P}^{x_{0}}|\mathbb{Q}_{\phi}^{x_{0}})=\frac{1}{2}{E}^{x_{0}}\left\\{\int_{0}^{T}\left\|s_{\phi}(t,X_{t})-s(t,X_{t})\right\|_{\Sigma(t,X_{t})}^{2}\mathrm{d}t\right\\}=\frac{1}{2}e(\phi).$
(12)

Hence $\textsc{kl}(\mathbb{P}^{x_{0}}|\mathbb{Q}_{\phi}^{x_{0}})<\infty$ under
the assumption that $e(\phi)<\infty$. As this implies $\mathbb{P}^{x_{0}}$ is
absolutely continuous with respect $\mathbb{Q}_{\phi}^{x_{0}}$, the event
$Z_{T}=x_{0}$ holds $\mathbb{Q}_{\phi}^{x_{0}}$-almost surely as it is the
case under $\mathbb{P}^{x_{0}}$. ∎

###### Proof of Proposition 2.

By expanding the square in (12), we can decompose the Kullback–Leibler
divergence as

$\displaystyle\textsc{kl}(\mathbb{P}^{x_{0}}|\mathbb{Q}_{\phi}^{x_{0}})=C_{1}+C_{2}-C_{3},$
(13)

where

$\displaystyle C_{1}$
$\displaystyle=\frac{1}{2}\int_{0}^{T}\left\|s(t,x_{t})\right\|_{\Sigma(t,x_{t})}^{2}p(t,x_{t}\mid
0,x_{0})\mathrm{d}x_{t}\mathrm{d}t,$ $\displaystyle C_{2}$
$\displaystyle=\frac{1}{2}\sum_{m=1}^{M}\int_{t_{m-1}}^{t_{m}}\int_{\mathbb{R}^{d}}\|s_{\phi}(t,x_{t})\|_{\Sigma(t,x_{t})}^{2}p(t,x_{t}\mid
0,x_{0})\mathrm{d}x_{t}\mathrm{d}t,$ (14) $\displaystyle C_{3}$
$\displaystyle=\sum_{m=1}^{M}\int_{t_{m-1}}^{t_{m}}\int_{\mathbb{R}^{d}}\left\langle
s_{\phi}(t,x_{t}),s(t,x_{t})\right\rangle_{\Sigma(t,x_{t})}p(t,x_{t}\mid
0,x_{0})\mathrm{d}x_{t}\mathrm{d}t.$

We examine the term $C_{3}$ that depends on the unknown score function
$s(t,x_{t})$. Firstly, we can write

$\displaystyle
C_{3}=\sum_{m=1}^{M}\int_{t_{m-1}}^{t_{m}}\int_{\mathbb{R}^{d}}\left\langle
s_{\phi}(t,x_{t}),\nabla p(t,x_{t}\mid
0,x_{0})\right\rangle_{\Sigma(t,x_{t})}\mathrm{d}x_{t}\mathrm{d}t.$

By differentiating the Chapman–Kolmogorov equation with respect to the
variable $x_{t}$

$\displaystyle\nabla p(t,x_{t}\mid 0,x_{0})$
$\displaystyle=\nabla\int_{\mathbb{R}^{d}}p(t,x_{t}\mid
t_{m-1},x_{t_{m-1}})p(t_{m-1},x_{t_{m-1}}\mid 0,x_{0})\mathrm{d}x_{t_{m-1}}$
$\displaystyle=\int_{\mathbb{R}^{d}}\nabla\log p(t,x_{t}\mid
t_{m-1},x_{t_{m-1}})p(t,x_{t}\mid
t_{m-1},x_{t_{m-1}})p(t_{m-1},x_{t_{m-1}}\mid 0,x_{0})\mathrm{d}x_{t_{m-1}},$

we obtain

$\displaystyle
C_{3}=\sum_{m=1}^{M}\int_{t_{m-1}}^{t_{m}}{E}^{x_{0}}\left\\{\left\langle
s_{\phi}(t,X_{t}),\nabla\log p(t,X_{t}\mid
t_{m-1},X_{t_{m-1}})\right\rangle_{\Sigma(t,X_{t})}\right\\}\mathrm{d}t.$ (15)

By expanding the square in (8) and using (14), (15), and

$\displaystyle C_{4}$
$\displaystyle=\frac{1}{2}\sum_{m=1}^{M}\int_{t_{m-1}}^{t_{m}}{E}^{x_{0}}\left\\{\|g(t_{m-1},X_{t_{m-1}},t,X_{t})\|_{\Sigma(t,X_{t})}^{2}\right\\}\mathrm{d}t,$

we have

$\displaystyle L(\phi)=C_{2}-C_{3}+C_{4}.$

The claim follows by noting the decomposition in (13) and taking
$C=C_{1}-C_{4}$. ∎

##  Appendix B Learning Doob’s $h$-transform

Suppose we have found a minimizer $\hat{\phi}\in\arg\min_{\phi\in\Phi}L(\phi)$
satisfying $e(\hat{\phi})<\infty$ and denote the corresponding score
approximation as $\hat{s}(t,x_{t})=s_{\hat{\phi}}(t,x_{t})$ and drift function
as $\hat{b}(t,z_{t})=b_{\hat{\phi}}(t,z_{t})$. Consider a time-reversed bridge
process $\hat{Z}=(\hat{Z}_{t})_{t\in[0,T]}$ satisfying

$\displaystyle\mathrm{d}\hat{Z}_{t}=\hat{b}(t,\hat{Z}_{t})\mathrm{d}t+\sigma(T-t,\hat{Z}_{t})\mathrm{d}B_{t},\quad\hat{Z}_{0}=x_{T},$
(16)

which should be seen as an approximation of (4). Let $\hat{q}(t,z_{t}\mid
s,z_{s})$ denote its transition density for any $0\leq s<t\leq T$,
$\hat{\mathbb{Q}}^{x_{0},x_{T}}$ be the induced path measure, and
$\hat{E}^{x_{0},x_{T}}$ to denote expectation with respect to
$\hat{\mathbb{Q}}^{x_{0},x_{T}}$. Note that
$\hat{p}^{\star}(t,x_{t})=\hat{q}(T-t,x_{t}\mid 0,x_{T})$ is an approximation
of the marginal density $p^{\star}(t,x_{t})$ in (5) for each $t\in(0,T)$.

Our discussion in (10) prompts having the time-reversal of (16) as an
approximation of the Doob’s $h$-transform process $X^{\star}$ in (2). The
bridge process $\hat{X}=(\hat{X}_{t})_{t\in[0,T]}=(\hat{Z}_{T-t})_{t\in[0,T]}$
satisfies

$\displaystyle\mathrm{d}\hat{X}_{t}=\hat{f}(t,\hat{X}_{t})\mathrm{d}t+\sigma(t,\hat{X}_{t})\mathrm{d}W_{t},\quad\hat{X}_{0}=x_{0},$

with drift function
$\hat{f}(t,x_{t})=-\hat{b}(T-t,x_{t})+\Sigma(t,x_{t})\hat{s}^{\star}(t,x_{t})+\nabla\cdot\Sigma(t,x_{t})$
which is to be seen as an approximation of $f^{\star}(t,x_{t})$. We can
approximate the score
$\hat{s}^{\star}(t,x_{t})=\nabla\log\hat{p}^{\star}(t,x_{t})$ of the marginal
density $\hat{p}^{\star}(t,x_{t})$ and hence the time-reversal of (16) using
the methodology described in Section 2.2. The following summarizes the key
elements involved.

Consider a path measure $\hat{\mathbb{Q}}_{\phi}^{x_{0},x_{T}}$ that is
induced by the bridge process $\hat{X}=(\hat{X}_{t})_{t\in[0,T]}$ satisfying

$\displaystyle\mathrm{d}\hat{X}_{t}=\hat{f}_{\phi}(t,\hat{X}_{t})\mathrm{d}t+\sigma(t,\hat{X}_{t})\mathrm{d}W_{t},\quad\hat{X}_{0}=x_{0},$
(17)

with drift function
$\hat{f}_{\phi}(t,x_{t})=-\hat{b}(T-t,x_{t})+\Sigma(t,x_{t})\hat{s}^{\star}_{\phi}(t,x_{t})+\nabla\cdot\Sigma(t,x_{t})$,
where
$\hat{s}_{\phi}^{\star}:[0,T]\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$
denotes a function approximation of the score $\hat{s}^{\star}(t,x_{t})$ with
parameters $\phi\in\Phi$ to be optimized. We now measure the score
approximation error as

$\displaystyle\hat{e}(\phi)=\hat{E}^{x_{0},x_{T}}\left\\{\int_{0}^{T}\left\|\hat{s}_{\phi}^{\star}(t,\hat{X}_{t})-\hat{s}^{\star}(t,\hat{X}_{t})\right\|_{\Sigma(t,\hat{X}_{t})}^{2}\mathrm{d}t\right\\}.$

######  Proposition 3.

Assuming $\hat{e}(\phi)<\infty$, we have
$\textsc{kl}(\hat{\mathbb{Q}}^{x_{0},x_{T}}|\hat{\mathbb{Q}}_{\phi}^{x_{0},x_{T}})=\hat{e}(\phi)/2$
and $\hat{X}_{T}=x_{T}$ holds $\hat{\mathbb{Q}}_{\phi}^{x_{0},x_{T}}$-almost
surely.

######  Proposition 4.

For any partition $(t_{m})_{m=0}^{M}$ of the interval $[0,T]$, we have
$\textsc{kl}(\hat{\mathbb{Q}}^{x_{0},x_{T}}|\hat{\mathbb{Q}}_{\phi}^{x_{0},x_{T}})=\hat{L}(\phi)+\hat{C}$
if $\hat{e}(\phi)<\infty$, where $\hat{C}$ is a constant independent of
$\phi\in\Phi$, the loss function is defined as

$\displaystyle\hat{L}(\phi)=\frac{1}{2}\sum_{m=1}^{M}\int_{t_{m-1}}^{t_{m}}\hat{E}^{x_{0},x_{T}}\left\\{\left\|\hat{s}_{\phi}^{\star}(T-t,\hat{Z}_{t})-\hat{g}(t_{m-1},\hat{Z}_{t_{m-1}},t,\hat{Z}_{t})\right\|_{\Sigma(t,\hat{Z}_{t})}^{2}\right\\}\mathrm{d}t,$
(18)

and $\hat{g}(s,z_{s},t,z_{t})=\nabla\log\hat{q}(t,z_{t}\mid s,z_{s})$ for
$0\leq s<t\leq T$ and $z_{s},z_{t}\in\mathbb{R}^{d}$.

The proof of these results is similar to Appendix A and is thus omitted. As
before, this allows us to circumvent intractability in the Kullback–Leibler
divergence
$\textsc{kl}(\hat{\mathbb{Q}}^{x_{0},x_{T}}|\hat{\mathbb{Q}}_{\phi}^{x_{0},x_{T}})$
by minimizing the loss function $\hat{L}(\phi)$. In the ideal case of
$\hat{s}_{\phi}^{\star}(t,x_{t})=\hat{s}^{\star}(t,x_{t})$
$\hat{\mathbb{Q}}^{x_{0},x_{T}}$-almost surely, the minimal loss of
$\hat{L}(\phi)=-\hat{C}$ is also unknown in practice, and
$\hat{\mathbb{Q}}_{\phi}^{x_{0},x_{T}}$ recovers the law of the diffusion
bridge process $X^{\star}$ only if the initial score approximation error
satisfies $e(\hat{\phi})=0$.

Given a minimizer $\phi^{\circ}\in\arg\min_{\phi\in\Phi}\hat{L}(\phi)$ and the
corresponding score approximation
$s^{\circ}(t,x_{t})=\hat{s}_{\phi^{\circ}}^{\star}(t,x_{t})$, by rewriting the
drift $f^{\circ}(t,x_{t})=\hat{f}_{\phi^{\circ}}(t,x_{t})$ as

$\displaystyle
f^{\circ}(t,x_{t})=f(t,x_{t})+\Sigma(t,x_{t})\\{s^{\circ}(t,x_{t})-\hat{s}(t,x_{t})\\}$

and comparing it with (10), we see that the last two terms on the right
provide an approximation of the term $\Sigma(t,x_{t})\nabla\log h(t,x_{t})$ in
Doob’s $h$-transform.

##  Appendix C Numerical implementation

In this section, we detail various numerical considerations to implement our
proposed methodology. For simplicity, we employ the Euler–Maruyama scheme
(Kloeden and Platen, 1992, p. 340) on a uniform discretization of the interval
$[0,T]$, denoted by $0=t_{0}<t_{1}<\cdots<t_{M}=T$, with stepsize $\delta
t=T/M$ for $M>1$. Non-uniform discretizations involve only minor
modifications; some higher-order schemes could also be considered. In the
following, we denote a multivariate normal distribution with mean vector
$\mu\in\mathbb{R}^{d}$ and covariance matrix $\Sigma\in\mathbb{R}^{d\times d}$
as $\mathcal{N}(\mu,\Sigma)$, and its density as
$x\mapsto\mathcal{N}(x;\mu,\Sigma)$. We write the zero vector as
$0_{d}\in\mathbb{R}^{d}$ and the identity matrix as
$I_{d}\in\mathbb{R}^{d\times d}$.

The time-discretization of the stochastic differential equation defining $X$
satisfies the following recursion

$\displaystyle X_{t_{m}}=X_{t_{m-1}}+\delta
tf(t_{m-1},X_{t_{m-1}})+\sigma(t_{m-1},X_{t_{m-1}})(W_{t_{m}}-W_{t_{m-1}}),\quad
X_{0}=x_{0},$ (19)

for $m\in\\{1,\ldots,M\\}$, with independent Brownian increments
$W_{t_{m}}-W_{t_{m-1}}\sim\mathcal{N}(0_{d},\delta tI_{d})$. Equation (19)
induces a normal approximation of the transition density
$p(t_{m},x_{t_{m}}\mid t_{m-1},x_{t_{m-1}})$ of the form

$\displaystyle p_{M}(t_{m},x_{t_{m}}\mid
t_{m-1},x_{t_{m-1}})=\mathcal{N}\\{x_{t_{m}};x_{t_{m-1}}+\delta
tf(t_{m-1},x_{t_{m-1}}),\delta t\Sigma(t_{m-1},x_{t_{m-1}})\\}.$

By replacing the score of $p(t_{m},x_{t_{m}}\mid t_{m-1},x_{t_{m-1}})$ with
that of $p_{M}(t_{m},x_{t_{m}}\mid t_{m-1},x_{t_{m-1}})$, the gradient
$g(t_{m-1},x_{t_{m-1}},t_{m},x_{t_{m}})$ in Proposition 2 can be approximated
by

$\displaystyle g_{M}(t_{m-1},x_{t_{m-1}},t_{m},x_{t_{m}})=-(\delta
t)^{-1}\Sigma^{-1}(t_{m-1},x_{t_{m-1}})\\{x_{t_{m}}-x_{t_{m-1}}-\delta
tf(t_{m-1},x_{t_{m-1}})\\}.$

We then define the following approximation of the loss function $L(\phi)$ in
(8)

$L_{M}(\phi)=\frac{1}{2}\delta
t\sum_{m=1}^{M}{E}_{M}^{x_{0}}\left\\{\left\|s_{\phi}(t_{m},X_{t_{m}})-g_{M}(t_{m-1},X_{t_{m-1}},t_{m},X_{t_{m}})\right\|_{\Sigma(t_{m},X_{t_{m}})}^{2}\right\\},$

where ${E}_{M}^{x_{0}}$ denotes expectation with respect to the law of the
time-discretized process under (19). To obtain a minimizer
$\hat{\phi}\in\arg\min_{\phi\in\Phi}L_{M}(\phi)$ using stochastic gradient
algorithms, the gradient with respect to parameters $\phi\in\Phi$

$\displaystyle\nabla_{\phi}L_{M}(\phi)=\delta
t\sum_{m=1}^{M}E_{M}^{x_{0}}\left[(\nabla_{\phi}s_{\phi}^{\mathrm{\scriptscriptstyle
T}}\Sigma)(t_{m},X_{t_{m}})\\{s_{\phi}(t_{m},X_{t_{m}})-g_{M}(t_{m-1},X_{t_{m-1}},t_{m},X_{t_{m}})\\}\right]$
(20)

can be unbiasedly estimated using independent sample paths from (19). The
above notation $\nabla_{\phi}s_{\phi}$ refers to the Jacobian of $s_{\phi}$.
Equation (20) can be seen as an approximation of the gradient
$\nabla_{\phi}{L}(\phi)$.

After obtaining $\hat{\phi}$ with optimization, we can simulate a proposal
bridge $\hat{Z}=(\hat{Z}_{t})_{t\in[0,T]}$ satisfying (16) with drift function
$\hat{b}(t,z_{t})=b_{\hat{\phi}}(t,z_{t})$. We employ the following modified
Euler–Maruyama scheme

$\displaystyle\hat{Z}_{t_{m}}=\hat{Z}_{t_{m-1}}+\delta
t\hat{b}(t_{m-1},\hat{Z}_{t_{m-1}})+\left(\frac{T-t_{m}}{T-t_{m-1}}\right)^{1/2}\sigma(T-t_{m-1},\hat{Z}_{t_{m-1}})(B_{t_{m}}-B_{t_{m-1}}),$
(21)

for $m\in\\{1,\ldots,M-1\\}$, with independent Brownian increments
$B_{t_{m}}-B_{t_{m-1}}\sim\mathcal{N}(0_{d},\delta tI_{d})$, initial condition
$\hat{Z}_{0}=x_{T}$, and terminal constraint $\hat{Z}_{T}=x_{0}$. This changes
the variance of the usual Euler–Maruyama transitions with a multiplier of
$(T-t_{m})/(T-t_{m-1})$ at time step $m$. We found that this modification can
improve practical performance for times near the endpoint by lowering the
transition variances. Such behaviour is consistent with findings in earlier
works by Durham and Gallant (2002) and Papaspiliopoulos et al. (2013) when
constructing proposal bridge processes with the drift of a Brownian bridge.
This gives a normal approximation of the transition density
$\hat{q}(t_{m},z_{t_{m}}\mid t_{m-1},z_{t_{m-1}})$

$\displaystyle\hat{q}_{M}(t_{m},z_{t_{m}}\mid t_{m-1},z_{t_{m-1}})$
$\displaystyle=\mathcal{N}\left\\{z_{t_{m}};z_{t_{m-1}}+\delta
t\hat{b}(t_{m-1},z_{t_{m-1}}),\delta
t\left(\frac{T-t_{m}}{T-t_{m-1}}\right)\Sigma(T-t_{m-1},z_{t_{m-1}})\right\\}.$
(22)

We can perform importance sampling on $\mathbb{E}_{M}=(\mathbb{R}^{d})^{M-1}$
to correct for the discrepancy between the law of our proposal bridge process

$\displaystyle\hat{Q}_{M}^{x_{0},x_{T}}(z_{t_{1}:t_{M-1}})=\prod_{m=1}^{M-1}\hat{q}_{M}(t_{m},z_{t_{m}}\mid
t_{m-1},z_{t_{m-1}}),\quad
z_{t_{1}:t_{M-1}}=(z_{t_{m}})_{m=1}^{M-1}\in\mathbb{E}_{M},$ (23)

and the law of the time-discretized diffusion bridge process

$\displaystyle
P_{M}^{x_{0},x_{T}}(x_{t_{1}:t_{M-1}})=\frac{\gamma_{M}^{x_{0},x_{T}}(x_{t_{1}:t_{M-1}})}{p_{M}(T,x_{T}\mid
0,x_{0})},\quad x_{t_{1}:t_{M-1}}=(x_{t_{m}})_{m=1}^{M-1}\in\mathbb{E}_{M},$
(24)

with
$\gamma_{M}^{x_{0},x_{T}}(x_{t_{1}:t_{M-1}})=\prod_{m=1}^{M}p_{M}(t_{m},x_{t_{m}}\mid
t_{m-1},x_{t_{m-1}})$, and also estimate

$\displaystyle p_{M}(T,x_{T}\mid
0,x_{0})=\int_{\mathbb{E}_{M}}\gamma_{M}^{x_{0},x_{T}}(x_{t_{1}:t_{M-1}})\mathrm{d}x_{t_{1}:t_{M-1}},$
(25)

which is an approximation of the transition density $p(T,x_{T}\mid 0,x_{0})$
under the Euler–Maruyama scheme. The corresponding unnormalized importance
weight is

$\omega(z_{t_{1}:t_{M-1}})=\frac{\gamma_{M}^{x_{0},x_{T}}(x_{t_{1}:t_{M-1}})}{\hat{Q}_{M}^{x_{0},x_{T}}(z_{t_{1}:t_{M-1}})}$

with $x_{t_{1}:t_{M-1}}=(z_{T-t_{m}})_{m=1}^{M-1}$, and an unbiased importance
sampling estimator of the transition density $p_{M}(T,x_{T}\mid 0,x_{0})$ is
$N^{-1}\sum_{n=1}^{N}\omega(z_{t_{1}:t_{M-1}}^{n})$ where
$(z_{t_{1}:t_{M-1}}^{n})_{n=1}^{N}$ denote $N\in\mathbb{N}$ independent sample
paths from $\hat{Q}_{M}^{x_{0},x_{T}}$. As noted by Lin et al. (2010), the
root mean squared error of this transition density estimator is approximately
equals to the $\chi^{2}$-divergence of $\hat{Q}_{M}^{x_{0},x_{T}}$ from
$P_{M}^{x_{0},x_{T}}$ divided by the sample size $N$. One can also employ
proposals from (23) within an independent Metropolis–Hastings algorithm that
has (24) as its invariant law (Elerian et al., 2001). At each iteration of the
algorithm, a sample path
$z_{t_{1}:t_{M-1}}^{\circ}\sim\hat{Q}_{M}^{x_{0},x_{T}}$ is accepted with
probability
$\min\\{1,\omega(z_{t_{1}:t_{M-1}}^{\circ})/\omega(z_{t_{1}:t_{M-1}})\\}$,
where $z_{t_{1}:t_{M-1}}$ denotes the current state of the Markov chain. The
efficiency of this Markov chain Monte Carlo algorithm can be assessed by
monitoring its acceptance probability. To improve the acceptance probability,
we can also combine independent Metropolis–Hastings with importance sampling
within a particle independent Metropolis–Hastings algorithm (Andrieu et al.,
2010) that has invariant law (24). Each iteration of this algorithm involves
selecting a proposed sample path $z_{t_{1}:t_{M-1}}^{\circ}$ among $N$
candidates $(z_{t_{1}:t_{M-1}}^{n})_{n=1}^{N}\sim\hat{Q}_{M}^{x_{0},x_{T}}$
according to probabilities proportional to their weights
$(\omega(z_{t_{1}:t_{M-1}}^{n}))_{n=1}^{N}$, and accepting it with probability
$\min\\{1,\hat{p}_{M}^{\circ}(T,x_{T}\mid 0,x_{0})/\hat{p}_{M}(T,x_{T}\mid
0,x_{0})\\}$ that depends on the ratio of the new and current transition
density estimators $\hat{p}_{M}^{\circ}(T,x_{T}\mid
0,x_{0})=N^{-1}\sum_{n=1}^{N}\omega(z_{t_{1}:t_{M-1}}^{n})$ and
$\hat{p}_{M}(T,x_{T}\mid 0,x_{0})$, respectively. Under mild assumptions,
consistency of importance sampling estimators as $N\rightarrow\infty$ implies
that the acceptance probability of particle independent Metropolis–Hastings
algorithm converges to one. This algorithm can also be combined with unbiased
Markov chain Monte Carlo methods to provide unbiased estimates of expectations
with respect to the law of the time-discretized diffusion bridge (Middleton et
al., 2019).

Lastly, we sketch the key steps to learn the Doob’s $h$-transform process for
the sake of brevity. Using the score of the normal transition density in (C),
we may approximate $\hat{g}(t_{m-1},z_{t_{m-1}},t_{m},z_{t_{m}})$ and hence
the loss function $\hat{L}(\phi)$ in (18). The approximate loss can be
minimized using stochastic gradient algorithms and sample paths from (21). By
time-discretizing the resulting proposal bridge process in (17), we may then
employ it as an importance proposal to approximate the law in (24) and the
transition density in (25), or to generate proposal distributions within
independent Metropolis–Hastings algorithms.

##  Appendix D Implementation details

###  D.1 Neural network and stochastic optimization

The architecture of the neural networks we employed is illustrated in Fig. 5.
For all numerical experiments, optimization was performed using the stochastic
gradient algorithm of Kingma and Ba (2014) with a momentum of $0.99$ and
learning rate of $0.01$.

Figure 5: Neural network architectures involve multi-layer perceptron (MLP)
blocks and an “Encoding” block which applies the sine transform described in
Vaswani et al. (2017). MLP (1a) and (1b) have one hidden layer and MLP (2) has
two hidden layers. All neurons use the Leaky ReLU activation function.

###  D.2 Ornstein–Uhlenbeck process

Let $X$ be an Ornstein–Uhlenbeck process, defined by (1) with linear drift
function $f(t,x_{t})=\alpha-\beta x_{t}$ and identity diffusion coefficient
$\sigma(t,x_{t})=I_{d}$. In this analytically tractable example, for any
$0\leq s<t\leq T$, the transition density of $X$ is a normal density
$p(t,x_{t}\mid s,x_{s})=\mathcal{N}\\{x_{t};m(t-s,x_{s}),v(t-s)I_{d}\\}$, with
the following mean and variance

$\displaystyle
m(t-s,x_{s})=\frac{\alpha}{\beta}+\left(x_{s}-\frac{\alpha}{\beta}\right)\exp\\{-\beta(t-s)\\},\quad
v(t-s)=\frac{1-\exp\\{-2\beta(t-s)\\}}{2\beta}.$

Hence the logarithmic gradient term in the Doob’s $h$-transform of (2) is

$\displaystyle\nabla\log
h(t,x_{t})=\frac{\exp\\{-\beta(T-t)\\}}{v(T-t)}\\{x_{T}-m(T-t,x_{t})\\},$

and the score of the transition density in (4) is

$\displaystyle s(t,x_{t})=v(t)^{-1}\\{m(t,x_{0})-x_{t}\\}.$

For dimension $d=1$ and varying time horizon $T\in\\{1,2,4,8\\}$, we employed
a time-discretization stepsize of $\delta t=0.02$, $500$ optimization
iterations, and $100$ sample paths per iteration. For the case of $T=1$ and
varying $d\in\\{1,2,4,8\\}$, we decreased the stepsize and increased the
number of optimization iterations and the capacity of the neural network with
dimension.

###  D.3 Interest rates model

We consider a special case of an interest rates model in Aït-Sahalia and Lo
(1998), defined by (1) with drift function $f(t,x_{t})=\theta/x_{t}-x_{t}$ and
diffusion coefficient $\sigma(t,x_{t})=1$. The transition density of this
diffusion is known and given by

$\displaystyle\log p(t,x_{t}\mid s,x_{s})$
$\displaystyle=\theta\log(x_{t}/x_{s})+\frac{1}{2}\log(x_{t}x_{s})-x_{t}^{2}+\left(\theta+\frac{1}{2}\right)(t-s)$
$\displaystyle-\log\mathrm{sinh}(t-s)-\frac{x_{t}^{2}+x_{s}^{2}}{\exp\\{2(t-s)\\}-1}+\log
I_{\theta-1/2}\left\\{\frac{x_{t}x_{s}}{\mathrm{sinh}(t-s)}\right\\},$

for $0\leq s<t\leq T$, where $I_{\nu}$ denotes the modified Bessel function of
order $\nu$. The logarithmic gradient term in the Doob’s $h$-transform of (2)
is

$\displaystyle\nabla\log h(t,x_{t})$
$\displaystyle=-\frac{\theta}{x_{t}}+\frac{1}{2x_{t}}-\frac{2x_{t}}{\exp\\{2(T-t)\\}-1}$
$\displaystyle+I_{\theta-1/2}\left\\{\frac{x_{T}x_{t}}{\mathrm{sinh}(T-t)}\right\\}^{-1}J_{\theta-1/2}\left\\{\frac{x_{T}x_{t}}{\mathrm{sinh}(T-t)}\right\\}\frac{x_{T}}{\mathrm{sinh}(T-t)},$

and the score of the transition density in (4) is

$\displaystyle
s(t,x_{t})=\frac{\theta}{x_{t}}+\frac{1}{2x_{t}}-2x_{t}-\frac{2x_{t}}{\exp(2t)-1}+I_{\theta-1/2}\left\\{\frac{x_{t}x_{0}}{\mathrm{sinh}(t)}\right\\}^{-1}J_{\theta-1/2}\left\\{\frac{x_{t}x_{0}}{\mathrm{sinh}(t)}\right\\}\frac{x_{0}}{\mathrm{sinh}(t)},$

where $J_{\nu}$ denotes the derivative of $I_{\nu}$.

Numerical experiments for all $T\in\\{1,2,4,8\\}$ employed a time-
discretization stepsize of $\delta t=0.02$, $1000$ optimization iterations,
$1000$ sample paths per iteration with $10$ unique initial conditions
$X_{0}=x_{0}$ sampled from the gamma distribution with shape $5$ and rate $2$.

###  D.4 Cell model

Our numerical experiments for all $T\in\\{2,4,8,16\\}$ and
$\sigma_{X}^{2}\in\\{0.1,1\\}$ employed a time-discretization stepsize of
$\delta t=0.02$, $2000$ optimization iterations, and $100$ sample paths per
iteration.


# RamanNet: A generalized neural network architecture for Raman Spectrum
Analysis

Nabil Ibtehaz  Department of Computer Science, Perdue University, West
Lafayette, IN 47907, United States<EMAIL_ADDRESS>(NI)  Muhammad E. H.
Chowdhury  Department of Electrical Engineering, Qatar University, Doha-2713,
Qatar<EMAIL_ADDRESS>(MEHC<EMAIL_ADDRESS>(AK<EMAIL_ADDRESS>(SK)  Corresponding author  Amith Khandakar  Department of Electrical
Engineering, Qatar University, Doha-2713, Qatar<EMAIL_ADDRESS>(MEHC),
<EMAIL_ADDRESS>(AK<EMAIL_ADDRESS>(SK)  Susu M. Zughaier  Department
of Basic Medical Sciences, College of Medicine, QU Health, Qatar University,
Doha, Qatar, PO Box 2713<EMAIL_ADDRESS>(SMZ)  Serkan Kiranyaz
Department of Electrical Engineering, Qatar University, Doha-2713, Qatar;
<EMAIL_ADDRESS>(MEHC<EMAIL_ADDRESS>(AK<EMAIL_ADDRESS>(SK)
M. Sohel Rahman  Department of Computer Science and Engineering, Bangladesh
University of Engineering and Technology, Dhaka-1205, Bangladesh,
<EMAIL_ADDRESS>(MSR)

###### Abstract

Raman spectroscopy provides a vibrational profile of the molecules and thus
can be used to uniquely identify different kind of materials. This sort of
fingerprinting molecules has thus led to widespread application of Raman
spectrum in various fields like medical dignostics, forensics, mineralogy,
bacteriology and virology etc. Despite the recent rise in Raman spectra data
volume, there has not been any significant effort in developing generalized
machine learning methods for Raman spectra analysis. We examine, experiment
and evaluate existing methods and conjecture that neither current sequential
models nor traditional machine learning models are satisfactorily sufficient
to analyze Raman spectra. Both has their perks and pitfalls, therefore we
attempt to mix the best of both worlds and propose a novel network
architecture RamanNet. RamanNet is immune to invariance property in CNN and at
the same time better than traditional machine learning models for the
inclusion of sparse connectivity. Our experiments on 4 public datasets
demonstrate superior performance over the much complex state-of-the-art
methods and thus RamanNet has the potential to become the defacto standard in
Raman spectra data analysis.

##  1 Introduction

Raman scattering is one the various light-matter interaction, comprising
absorption and subsequent emission of light by matter [1]. Spectroscopy being
the study of interaction of light or broader electromagnetic radiation with
matter, thus has also focused on Raman scattering, ever since it’s discovery
in 1928 by Raman and Krishnan [2]. Unlike elastic scatering, e.g. Raleigh
scattering, the wavelength of incident light changes in Raman scattering [3].
Following the typical norm of inelastic light scattering, when a photon
excites the sample the electrons are raised to a higher virtual energy state
[4]. This excitation event is usually short-lived and the molecule soon
reaches a new stable energy state, either lower (Stokes shift) or higher
(anti-Stokes) [5]. Based on the difference in energy, the sample achieves a
different vibrational and rotational state. Therefore, Raman spectroscopy can
be used to analyze the vibrational modes of various molecules, extracting the
structural fingerprint of such materials in the process [1].

Being capable of uniquely fingerprinting materials, Raman spectra has been
used in a wide variety of applications, covering medical diagnoistics,
forensics, mineralogy, bacteriology and virology [6]. **ADDING CITATIONS**

Conventionally Raman spectra are analyzed in terms of wavenumber $\tilde{v}$
$(cm^{-1})$. The standard practice is to present them with the wavenumber
shift linearly increasing along the horizontal axis. On the contrary, the
vertical axis ordinate is proportional to intensity [5]. Therefore, this
setting is not apparently too different from traditional spectograms. However,
the issue in treating Raman Spectra as typical spectrograms is that for Raman
spectrum we don’t have time along the horizontal axis, thus we cannot apply
models that are capable of spectrum analysis, e.g., CNN. On the contrary,
although the traditional machine learning methods are more capable, they
suffer from curse of dimensionality as Raman spectra data is usually long.
Principal Component Analysis (PCA) has been widely used for feature reduction,
as 34 out of recent 52 papers used that [6], still it may not always be able
to compress Raman spectra as PCA is more suitable for tabular data source,
where the features are ideally correlated, but in Raman spectra the
intensities at nearby Raman shifts are expected to demonstrate correlation.

With the reduction of cost and complexity related to Raman data extraction
pipeline, there has been an unprecedented expansion in Raman datasets in
recent years. These sudden growth of available Raman data requires suitable
methods to analyze them properly. However, to the best of our knowledge, till
date there has not been any machine learning methodology devised solely
focused on Raman spectra data, considering the pattern and properties of this
unique data source. To this end, we carefully contemplate the properties of
Raman spectra data and corresponding attributes of contemporary machine
learning methods. We argue why application of convolutional neural networks
(CNN) may not be appropriate for invariance properties, but at the same time
acknowledge that CNN is more capable than traditional machine learning methods
for the sparse connection and weight sharing properties. We make an attempt to
fuse the best of both worlds and propose RamanNet, a deep learning model that
also follows the paradigm of sparse connection without the limitation of
temporal or spatial invariance. We achieve this by using shifted densely
connected layers and emulate sparse connectivity found in CNN, without the
concern of invariance. Furthermore, dimensionality reduction has been involved
in most Raman spectra applications, thus we employ Triplet loss in our hidden
layers and make more separable and informative embeddings in lower dimensional
space. Our experiments on 4 public datasets demonstrate superior performance
over the much complex state-of-the-art methods and thus RamanNet has the
potential to become the defacto standard in Raman spectra data analysis.

##  2 Motivations and High Level Considerations

From a visual perception, Raman spectrum resembles a signal-like waveform.
Which has led to the application of 1D Convolutional Neural Networks to
analyze Raman spectrums [7, 8, 9]. However, Raman spectrum is not a typical
signal rather it’s an energy distribution plot, which may prevent appropriate
utilization of CNNs. In this section, we briefly discuss the properties of
CNNs and present our motivations and rationale.

Convolutional Neural Networks are one of the most successful and widely used
neural network architectures, particularly in computer vision [10] and signal
processing [11] domains. The success of CNNs can be largely attributed to
three primary properties of CNNs, namely sparse interaction, parameter sharing
and equivariant representations [12]. The nature of image or signal data,
i.e., spatial or temporal properties, works in perfect harmony with the
properties of CNN. Unlike traditional feed forward neural networks, where all
the input features are connected with all the hidden or output layers, CNNs
employ sparse connectivity by leveraging kernels. Only a small portion of
input is analyzed at a specific step using kernels, thus this prevents the
computation from being overwhelmed with analyzing the entire input at once.
Although focus is put on local information, global information is also
considered by appropriate use of pooling operations along with broader field
of vision in deeper networks. Furthermore, the use of kernel based computation
also faciliatates parameter sharing, reducing number of paramaters and risk of
overfitting simultaneously. Another vital property of CNN is the equivariance
in representation. Mathematically for a convolutional operation $f$ on an
image $I$, $f(I(x,y))=f((I(x-h,y-k))$ , i.e., a point of interest whether it
resides at location $(x,y)$ or at a shifted location $(x-h,y-k)$, the output
of the convoluation operation remains the same.

This equivarience to translation plays a pivotal role in working image or
signal data. When processing signal data, this property implies that CNN
generates a timeline of the emergence of different keypoints in the signal. As
a result, regardless of the time of occurence, all the features in a signal is
captured, unlike traditional neural networks which would have only sought for
the features at the exact fixed timestamps. Furthermore, the application of
global pooling operations makes sure that all the feature signatures are
preserved in the final representation.

In Fig. 1, we have presented a simplified example of how CNN works with a
signal as input. In Fig. 1a, we have three shifted ECG signals, since we
cannot always be certain of the occurrence of a particular feature when
dealing with signals, it is imperative that the model can identify the feature
irrespective of the timestamp. In Fig. 1b, we can observe the featuremaps from
a convoluational layer, it is apparent that the equivalent features are
detected, albeit shifted in accordance with the shifted nature the signals,
complying to the property $f(I(x,y))=f(I(x-h,y-k))$. Finally, a global pooling
operation summarizes the featuremaps and provides identical representations of
all the three signals in Fig. 1c.

(a)

(b)

(c)

Figure 1: Analysis of CNN for a signal input.

From the example, it is obvious how applicable CNNs are for signal data. Since
the pattern of Raman spectrum resembles signal for quite a bit, it is trivial
to use CNNs to analyze raman spectra. The sparse connectivity and parameter
sharing truely helps in this regard, as it makes the computation simpler and
less prone to overfitting. However, the equivariance to translation which
proved vital for analyzing signals acts troublesome when working with Raman
spectrum. The Raman spectrum is plotted as intensity vs Raman shifts. Thus the
concept of equivariant translation is not applicable in this case, because
similar intensity at different Raman shifts implies a completely different
meaning. For example, in Fig. 2, we present a similar scenario but with Raman
spectrum as input. It can be seen that, we have three completely different
Raman spectrums, having similar patterns in different Raman shifts. But the
CNN model, due to translational equivariance, treats them as shifted inputs
and generates the same output for all of them (Fig. 2c).

Therefore, we hypothesize that the use of CNN in Raman spectrum analysis is
inappropriate. This is has been addressed in prior works using CNNs by using
deeper networks (broadening the field of vision thereby) along with discarding
pooling layers (attempting to preserve the locations of spectral peaks) [8].
On the contrary we can use traditional neural networks or machine learning
models which is capable of keeping track of the exact locations. However, they
succumb to the curse of dimensionality [13], due to the long length of Raman
spectrum. Although PCA has been mostly used to reduce that [6], still it feels
less intuitive to model spectrum inputs using hyperplane optimization, which
is more suitable for tabular data. As the spectra are apparently correlated in
neighboring Raman shifts, which motivates the use of sparse connectivity
through kernel-like operations in CNN instead. This conjecture is also
supported by the reduced accuracy in classical models [7, 8].

(a)

(b)

(c)

Figure 2: Analysis of CNN for a Raman spectrum input.

To sum up, we are faced with a dilemma between CNN and classical machine
learning pipelines. On one hand, CNNs help in analyzing Raman spectra for the
benefits of sparse connectivity (considering the correlation in neighboring
Raman shifts) and parameter sharing (reducing the number of parameters and
overfitting), but loses its applicability due to translational equivariance.
On the other hand, traditional models are free from translational
equivariance, but they can’t comprehend the correlation in the neighboring
raman shifts and suffers from the curse of dimensionality as trying to
optimize too many parameters at once.

Therefore, we propose a middle ground between the two approaches, fusing the
best of the both worlds. We propose to use shifted multi-layer perceptrons to
analyze shifted windows of Raman spectra. This facilitates sparse connectivity
as the shifted MLP layers mimic kernels and only analyze a part of the input.
Moreover, parameter sharing is redundant in this regard, as the kernel-like
operations of a particular kernel are only performed in one part and not
elsewhere. Finally, this also eliminates the issues with translation
equivarience, as for differnt locations or different shifted windows we have
different MLP layers. Thus,

Mathematically, a typical 1D CNN operation can be simplified as,

$y(i)=\sigma(\sum_{h}x(i+h)k(h)+b)$ (1)

Here, $x$ is the 1D input, $y$ is the output, $k$ is a learned kernel, $b$ is
the bias term and $\sigma$ is a nonlinear operation. The same kernel $k$ is
applied everywhere thus the translational equivarience is achieved.

On the contrary, our proposed modification from using a MLP,

$y(i)=\sigma(W_{f(i)}^{T}x+b)\equiv\sigma(\sum_{h}x(i+h)k_{f(i)}(h)+b)$ (2)

Here, the dot product $W^{T}x$ is mathematically equivalent to an 1D
convolutional operation with proper relation between the weight matrix $W$ and
kernel $k$. In addition since we are using sliding windows, the weight matrix
$W_{f(i)}$ and kernel $k_{f(i)}$ depends on the location, i.e., value of $i$.

##  3 Proposed Architecture

Hence, following our proposed modifications, we present RamanNet. As mentioned
in the previous section, we mimic the convolution operation using multi-layer
perceptrons or so called densely connected neural network layers. The input
Raman spectra is broken into overlaping sliding windows of length $w$ and step
size $dw$, and each of them is passed to a different dense block with $n_{1}$
neurons. This ensures sparse connectivity and reduces risk of overfitting.
Furthermore, this configuration somewhat resembles a convolution operation
without translation, thus we can consider the features extracted at the
neurons the same way as we would consider features learnt from kernels.

The features from all the dense blocks are concatenated together and a dropout
$dp_{1}$ is applied. The concatenated features are again summarized using
another dense layer with $n_{2}$ neurons. The outputs from the summarization
is regularized with a dropout of $dp_{2}$.

Finally, we compute $n_{f}$ features from the regulized outputs using another
dense layer with $n_{f}$ neurons. We name this layer Embedding Layer. Raman
spectra are known to be noisy with low SNR [8], which often leads to less
separation between the classes. Therefore, in order to improve this we
introduce Triplet Loss [14] as an auxiliary loss. Triplet Loss is defined
using eucleidian distance, $f$, between an anchor $A$, positive example $P$
and negative example $N$ as,

$L(A,P,N)=max(||f(A)-f(P)||^{2}-||f(A)-f(N)||^{2}+\alpha,0)$ (3)

Here, $\alpha$ is a margin term to ensure that the model learns non-trivial
information.

Using triplet loss we obtain a well separated embedding space. The embeddings
are finally used to predict the classes of input using Softmax activation
function and a dropout of $dp_{3}$.

Other than the output layer, all the layers uses LeakyRELU activation function
[15] and are batch normalized [16]. In our experiments we used the
hyperparameters as $w=50$, $dw=25$, $n_{1}=25$, $n_{2}=512$, $n_{f}=256$,
$dp_{1}=0.5$, $dp_{2}=0.4$, $dp_{3}=0.25$.

A simplified diagram of the RamanNet has been presented in Fig. 3.

Figure 3: RamanNet Architecture

##  4 Datasets

One particular limitation when working with deep learning for Raman spectrum
analysis is the lack of sufficient public benchmark datasets [6]. Although
there has been several of recent works utilizing deep learning models for
Raman spectrum analysis, the datasets are mostly proprietary or private [17,
18]. From an elaborate study of the recent works, investigating Raman Spectrum
analysis, we selected 4 publicly available datasets for our analysis and
benchmarking.

###  4.1 Covid Dataset

We used the publicly available data [19] from the recent work [20], which acts
as a pilot study of primary screening of COVID-19 by Raman spectroscopy. This
dataset contains a total of 177 serum samples collected from 63 COVID-19
patients, 59 suspected ones and 55 healthy people (i.e., control group). The
COVID-19 group was recruited at the Chengdu Public Health Clinical Medical
Center, and it includes 58 symptomatic and 5 asymptomatic patients. The
suspected group demonstrated ‘flu’ like symptoms but were tested negative but
at least to RT-PCR tests. For all the subjects, 1 hour repose of blood
sampling the serum was extracted by centrifuging at 3000 $rpm$ for 10 minutes,
and were stored at 4 ℃. Later, a single-mode diode laser with 785 $nm$
wavelength and 100 $mW$ power was used for Raman excitation. The laser power
on the sample was measured around 70 $mW$ and the spectra were recorded in the
range of 600-1800 $cm^{-1}$.

###  4.2 Melanoma Dataset

Erzina et al. [9] used Surface Enhanced Raman Spectroscopy to detect Melanoma.
Healthy residual skin and skin melanoma metastasis were collected. From cell
line and primary culture 12 categories of samples were considered, each having
8-9 samples. The SERS spectra were measured using ProRaman-L spectrometer at
785 $nm$ excitation wavelength and 33 $mW$ power. AuMs, functionalized by ADT-
NH2, ADT-COOH or ADT-(COOH)2 was used to measure the SERS spectra, i.e., 3
different spectra were obtained for each sample. Finally, the background was
removed using smooting algorithms and the recorded spectra were normalized to
0-1 intensity. The range of the recorded spectra is 100-4278 $cm^{-1}$.

###  4.3 Mineral Dataset

Mineral substance identification is another popular application of Raman
spectroscopy. Among the various mineral databases, we selected the RRUFF
database [21]. The RRUFF project aims at curating the most comprehensive set
of high quality spectral data from well characterized minerals, comprising
Raman spectra, X-ray diffraction and information from chemistry. The database
contains Raman spectra of various minerals at different configurations, i.e.
varied wavelengths and orientations, accompanied by various kinds of
processing. Furthermore, the Raman spectra computed from the different
materials are hardly consistent, e.g., the ranges of Raman shifts are also
different. In order to alleviate such irregularities, we have only considered
spectra computed at 432 $nm$, which is the majority. Moreover, we have only
collected the raw spectra, i.e., which were not processed any way. Since the
recorded spectra have different ranges of Raman shift, we have cropped the
region that is common in all the spectra and used cubic spline to interpolate
the locations if necessary.The resultant spectra cover the range 280 - 4237
$cm^{-1}$ and min-max normalization was performed. Finally, we take the
mineral classes with at least 10 samples, which reduces the database to twenty
mineral classes and use them to evaluate the proposed model.

###  4.4 Bacteria Dataset

We collected the Bacteria-ID dataset from [8], where the potential of Raman
spectrosocpy in label-free bacteria detection has been investigated. This
dataset consists of 30 bacterial and yeast isolates, including multiple
isolates of Gram-negative and Gram-positive bacteria. The dataset is organized
into a reference training dataset, reference fine-tuning set and test set. The
fine-tuning dataset is used to account for the changes in measurement caused
by optical system efficiency degradation. The training dataset contains 2000
spectra for each of the 30 isolates, whereas the fine-tuning and test set
contains 100 spectra for each isolate. The isolates were cultured on blood
agar plates sealed with Parafilm and stored at 4 ${}^{\circ}C$. The Raman
spectra was generated using Horiba LabRAM HR Evolution Raman microscope, with
633 $nm$ illumination at 13.17 $mW$ along with a 300 $l/mm$ grating. The
spectra was computed at 1.2 $cm^{-1}$ dispersion to simultaneously maximize
signal strength and minimize background signal. The recorded spectra was
normalized to $0-1$ intensity, covering the spectral range between
381.98-1792.4 $cm^{-1}$.

##  5 Experimental Setup

The experiments have been conducted in a server computer with Intel Xeon
@2.2GHz CPU, 24 GB RAM, and NVIDIA TESLA P100 (16 GB) GPU. We implemented the
RamanNet architecture using Tensorflow [22]. The codes are available in the
following github repository.

https://github.com/nibtehaz/RamanNet

In the following subsections, we briefly explain the experimental protocols
and evaluation procedures adopted for the different tasks. We have tried to
mimic the corresponding baseline papers to ensure a fair comparison. As
mentioned earlier, we were only compelled to follow a different evaluation
scheme for the Mineral dataset, thus we reproduced that baseline model’s
output following our exact protocol.

###  5.1 Covid Dataset

We followed the identical evaluation method as presented in [20]. In this
work, the authors conducted a “blind” validation. They randomly divided the
whole dataset into 70% for training and 30% as hold-out test set. In order to
further assess the independence of the data over model performance, this
process was repeated 50 times and the average values of the metrics were
recorded. We followed the same protocol to evaluate RamanNet. In order to
avoid overfitting, we used 10% data from the training set as validation data,
but the hold-out test data was left completely independent from the training
process and was only used for evaluation. The models were trained for 1000
epochs.

###  5.2 Melanoma Dataset

In the original work Erzina et al. [9], performed 75%:25% train-validation
split. Therefore, we attempted to follow a similar splitting criterion. Since,
the splitting information was not provided, we opted to perform a 4 fold
cross-validation instead, as it would also split the data in the same ratio.
The authors demonstrated 100% accuracy using 3 different spectra (AuMs
functionalized by ADT-NH2, ADT-COOH or ADT-(COOH)2 respectively) together.
However, to make the task difficult, we experimented with using only one type
of spectra as input. In order to achieve that, we reproduced the model as
described in [9] and evaluated the model with one particular spectra as input
at a time.

###  5.3 Mineral Dataset

Jichao et al. [7] also used the RRUFF database for mineral classification
task. However, they employed the leave one out cross validation scheme, which
is computationally too expensive. In order to reduce the computational
requirements, we thus opted for a 5 fold cross-validation instead. In order to
compare RamanNet with the model presented by Jichao et al. [7], we implemented
their proposed model and used it in our analysis.

###  5.4 Bacteria Dataset

In order to assess RamanNet on the Bacteria-ID dataset, we followed the same
training and evaluation procedure as presented in [8]. Similar to their
approach, we first pre-trained the model using the reference training dataset,
through a 5 fold cross-validation scheme. The five models obtained in this
process was then fine-tuned on the fine-tuning dateset, which was spilt into
90% training and 10% validation. The model with the highest accuracy on this
validation set was considered and was evaluated on the independent test
dataset. The models were trained for 100 and 250 epochs respectively on the
reference training and fine-tuning set.

##  6 Results

###  6.1 RamanNet Consistently Outperforms Existing Models

####  6.1.1 Covid Dataset

In [20], a SVM model was developed to distinguish the different categories,
namely healthy, suspected and COVID-19 patients. Instead of working with the
entire spectrum, wave points with significant differnces in ANOVA test were
selected. Thus, a statistically sound feature reduction was performed and the
reduced featureset was used as the input to the SVM model. On the contrary,
RamanNet takes the entire spectrum as input and adaptively finds the
significant region therein.

The average performance over 50 random trials for the different tasks are
presented in Table 1, here we present the accuracy, sensitivity and
specificity values. Among the 3 tasks, differentiating between suspected and
healthy subjects is the most challenging one, as evident from the inferior
performance of SVM (all metrics $\leq$ 70%). RamanNet, on the other hand,
performed comparatively better in this task. Notably, RamanNet improved
accuracy and specificity by 13% and 21% respectively. RamanNet also achieved
improved sensitivity score.

On the other two tasks the SVM model performed comparatively better. However,
our proposed RamanNet consistently outperforms SVM in all these two tasks as
well. Most promisingly, RamanNet improves sensitivity greatly in both the two
tasks. For this problem sensitivity is crucial as we need to correctly detect
the Covid-19 patients. This improvement in sensitivity did not come at any
cost of specificity, rather the specificity has also been improved compared to
the SVM model.

Table 1: Results on Covid Dataset Covid vs Suspected
---
Method | Accuracy | Sensitivity | Specificity
SVM | $87\pm 5$ | $89\pm 8$ | $86\pm 9$
RamanNet | $93\pm 3$ | $97\pm 4$ | $90\pm 6$
Covid vs Healthy
Method | Accuracy | Sensitivity | Specificity
SVM | $91\pm 4$ | $89\pm 7$ | $93\pm 6$
RamanNet | $95$ | $95\pm 4$ | $96\pm 3$
Suspected vs Healthy
Method | Accuracy | Sensitivity | Specificity
SVM | $69\pm 5$ | $70\pm 9$ | $66\pm 9$
RamanNet | $82\pm 6$ | $77\pm 15$ | $87\pm 11$

####  6.1.2 Melanoma Dataset

As described previously, we perform a 4 fold cross-validation on the melanoma
dataset. Following the original evaluation as performed by Erzina et al. [9],
we consider all the three different types of spectra simultaneously as input.
In addition, we perform a difficult version of the problem by taking only one
type of spectra as input at a time. The results are presented in Table 2.

Table 2: Results on Melanoma Dataset | -NH2 | -(COOH)2 | -COOH | All | #Parameters
---|---|---|---|---|---
Fold | RamanNet | CNN | RamanNet | CNN | RamanNet | CNN | RamanNet | CNN | RamanNet | CNN
1 | 100 | 97.42 | 100 | 100 | 99.35 | 94.19 | 100 | 100 | 1.3 M | 25.7M
2 | 99.35 | 96.13 | 99.35 | 98.71 | 98.71 | 96.12 | 100 | 98.71
3 | 100 | 95.45 | 100 | 98.05 | 99.35 | 87.01 | 100 | 100
4 | 100 | 97.40 | 100 | 98.70 | 96.75 | 96.10 | 100 | 99.35

From the results, it is evident that although the CNN model manages to achieve
perfect 100% accuracy (in 2 folds out of 4) when given 3 spectra as input, the
accuracy falls when a single spectra is given as input. This drop in
performance can be explained by the loss of information, when working with a
single spectra. For different functionalizations of AuMs, the sample is
observed from a different point of view and different information is obtained.
Thus, when working with reduced number of spectra, insightful information are
likely get lost and that negatively affects the performance. Although for the
-(COOH)2 as input the accuracy stays above 98%, it falls below 97% for the
other two input spectra.

RamanNet on the other hand seems to consistently outperform the CNN model for
both when all the 3 spectra are considered together or separately. RamanNet
not only consistently achieved 100% accuracy with all the 3 spectra as input,
but also with -NH2 and -(COOH)2 individual spectra as input. Only for -COOH
spectra as input, the performance was not up to the mark, but still the
performance was superior to the CNN. All these improvements become more
significant when we compare the number of the parameters of the two models.
The CNN model consists of 25.7 M parameters whereas RamanNet has only 1.3 M
parameters ($\sim 5\%$). Thus, RamanNet is not only more accurate, it is also
computationally efficient at the same time.

####  6.1.3 Mineral Dataset

As described in the previous sections, the mineral dataset consists of 20
mineral classes and in order to compare with the model proposed by [7], we
implement their model and perform a five fold cross-validation. The results
are presented in Table 3.

Table 3: Results on Mineral Dataset | Top 1 Accuracy | Top 5 Accuracy | Top 10 Accuracy | #Parameters
---|---|---|---|---
| RamanNet | CNN | RamanNet | CNN | RamanNet | CNN | RamanNet | CNN
Fold 1 | 87.5 | 81.25 | 95.83 | 87.5 | 100 | 93.75 | 1.3M | 6.6M
Fold 2 | 93.75 | 85.42 | 97.92 | 97.92 | 100 | 97.92
Fold 3 | 89.58 | 79.17 | 93.75 | 89.58 | 93.75 | 91.67
Fold 4 | 97.92 | 85.42 | 100 | 91.67 | 100 | 93.75
Fold 5 | 85.12 | 63.83 | 93.62 | 78.72 | 95.74 | 85.10

Here, we have presented the Top 1, Top 5 and Top 10 accuracy metrics, which
are popularly used for multiclass classification problems. It is evident that
RamanNet consistently outperforms the previous state of the art CNN model. For
Top 1 Accuracy, the improvement is more prominent, nevertheless RamanNet
performs better in regards to the other metrics as well. Another point worth
mentioning is that, even in this case the CNN model is heavier comparatively
(6.6 M parameters vs. 1.3 M parameters for RamanNet).

####  6.1.4 Bacteria Dataset

In the original work [8], the authors used a 25 layer deep residual
convolutional neural network to classify the bacteria isolate raman spectrum
to one of the 30 classes. The proposed model achieved an average accuracy of
82.2%, and the resulting confusion matrix is presented in Fig. 4(a). The
authors also experimented with common and popular classifiers like Logistic
Regression (LR) and Support Vector Machine (SVM) as baselines, but those
models could only reacg 75.7% and 74.9% accuracy respectively.

RamanNet on the other hand manages to achieve an average accuracy of 85.5% on
this dataset, despite only having 3 hidden layers. As shown in the confusion
matrix (Fig. 4(b)), the number of misclassifications has reduced. Although
erroneous predictions still exist, most Gram-positive and Gram-negative
bacteria have been misclassified as Gram-positive and Gram-negative bacteria
respectively, and the errors are also mostly confined within same genus, as
analyzed in [8]. Although this 3.3% improvement may seem minor, this is
achieve using a much shallower network (3 layers vs 25 layers). Furthermore,
the improvements become more apparent when we compare the individual
classifications next to each other. As presented in Fig. 4(c), RamanNet
achieves either equal or better accuracy for 25 out of 30 bacteria isolate
classes. For other than two isolate classes, ResNet is only better with a
small margin. On the contrary, the cases where RamanNet performed better it
surpassed ResNet with higher margin (e.g., MSSA 3, P. mirabilis, K. aerogenes
etc.).

(a) 25 Layer ResNet

(b) RamanNet

(c) Comparison

Figure 4:

###  6.2 Model Interpretation

Interpretability has been one of the focus in deep learning research in recent
years [23]. Deep learning models are competent function approximator and given
sufficient amount of data, they are capable of modelling almost any complex
functions. With this potential, also comes the concern of what the model is
actually learning from the data. The model can learn actual significant
information and perform prediction accordingly, or it can merely learn from
the noises and get confused from various confounding factors instead.
Therefore, it is imperative to investigate the model interpretability and
examine what the model is learning from.

For convolutional neural networks, we can use various methods like saliency
maps [24] or score-CAM [25] methods to study model interpretation. However,
for multi-layer perceptrons, it is non-trivial to do so. The various
visualization methods have been designed based on CNNs and they cannot be
directly translated to MLPs.

SHAP (SHapley Additive exPlanations) [26] is a game theoretic approach to
explain the output of machine learning models. The biggest advantage of using
SHAP is that it is model agnostic, thus it can be used to analyze any machine
learning model. We extracted the features extracted from different moving
windows of RamanNet and trained the identical top layer of RamanNet using
Scikit-Learn MLP implementation, for compatibility reasons of the SHAP
implementation [27]. We then computed the SHAP scores of all the -NH2 samples
in the melanoma dataset, the results are presented as a violin plot in Fig. 5.

Figure 5: SHAP values for RamanNet

From the violin plot it is evident that certain regions of the spectrum
contributes most to the prediction. Moreover, those regions correspond to
actual significant properties of the sample. For example, the region in
400-500 cm-1 corresponds to Fe containing proteins and contributes most to the
prediction. After this region the lipids (300-400 cm-1) contributes the most,
and it goes on. The unlabelled regions in the figure correspond to less
attributed regions, and it is apparent that the model was also aware enough to
put less focus on them. All these findings are consistent with the one
presented in [9]. Therefore, we can expect that RamanNet is learning
significant information from the data and ignoring the noises instead of
falling into conundrum with confounding factors.

##  7 Conclusion

Raman spectroscopy has slowly started to gain more attention with the advances
in SERS technolgy. The gradual decrease of cost and complexity in computing
raman spectrum, is paving the way to large scale raman spectrum data
collection for diverse tasks. Therefore, we need suitable machine learning
methods to analyze these large scale Raman spectrum data. However, there has
not yet been any model developed for the sole purpose of Raman spectroscopy
analysis, motivated and designed based on the properties of Raman spectrum.
Recently, existing methods like CNN or SVM has shown success in Raman spectrum
analysis, but in this work we presented reasonings that such methods may not
always be suitable for Raman spectra analysis.

In this work, we present RamanNet, a generalized neural network architecture
for Raman spectrum analysis. We take intuitions from the nature of data and
design our model accordingly. We propose modifications to the convolutional
network behaviors and emulate such operations using multi-layer perceptrons.
This adjustment brings the best out of both worlds and it is reflected in
evaluating the model on 4 public datasets. Not only, RamanNet outperforms all
the state-of-the-art approaches in Raman spectroscopy analysis, it achieves it
adopting less complexity. Moreover, RamanNet generates embeddings from
spectrum which is much better than what is obtained from PCA, the defacto
standard in Raman spectrum analysis. Furthermore, we examined the
interpretibility of RamanNet and it was reveled that the model is capable of
focusing on the actual important information.

The future direction of this research can be manifold. Firstly, we wish to
evaluate RamanNet on more large-scale datasets, as they become public.
Secondly, we have not performed any preprocessing or background removal of the
spectra, we wish to investigate this further with more dedicated study to
infer the denoising capabilities of RamanNet. Last but not least, we also wish
to experiment with multiple particle data to assess if RamanNet is capable of
identifying, segmenting and extracting the signatures of the different
particles.

## References

  * [1] Derek J Gardiner.  Introduction to raman scattering.  In Practical Raman Spectroscopy, pages 1–12. Springer, 1989.
  * [2] Chandrasekhara Venkata Raman and Kariamanikkam Srinivasa Krishnan.  A new type of secondary radiation.  Nature, 121(3048):501–502, 1928.
  * [3] Jeroen Cornel, Christian Lindenberg, Jochen Schöll, and Marco Mazzotti.  Raman spectroscopy.  Industrial Crystallization Process Monitoring and Control, pages 93–103, 2012.
  * [4] Gordon G Hammes.  Spectroscopy for the biological sciences.  John Wiley & Sons, 2005.
  * [5] Robin R Jones, David C Hooper, Liwu Zhang, Daniel Wolverson, and Ventsislav K Valev.  Raman techniques: fundamentals and frontiers.  Nanoscale research letters, 14(1):1–34, 2019.
  * [6] Félix Lussier, Vincent Thibault, Benjamin Charron, Gregory Q Wallace, and Jean-Francois Masson.  Deep learning and artificial intelligence methods for raman and surface-enhanced raman scattering.  TrAC Trends in Analytical Chemistry, 124:115796, 2020.
  * [7] Jinchao Liu, Margarita Osadchy, Lorna Ashton, Michael Foster, Christopher J Solomon, and Stuart J Gibson.  Deep convolutional neural networks for raman spectrum recognition: a unified solution.  Analyst, 142(21):4067–4074, 2017.
  * [8] Chi-Sing Ho, Neal Jean, Catherine A Hogan, Lena Blackmon, Stefanie S Jeffrey, Mark Holodniy, Niaz Banaei, Amr AE Saleh, Stefano Ermon, and Jennifer Dionne.  Rapid identification of pathogenic bacteria using raman spectroscopy and deep learning.  Nature communications, 10(1):1–8, 2019.
  * [9] M Erzina, A Trelin, O Guselnikova, B Dvorankova, K Strnadova, A Perminova, P Ulbrich, D Mares, V Jerabek, R Elashnikov, et al.  Precise cancer detection via the combination of functionalized sers surfaces and convolutional neural network with independent inputs.  Sensors and Actuators B: Chemical, 308:127660, 2020.
  * [10] Hyeon-Joong Yoo.  Deep convolution neural networks in computer vision: a review.  IEIE Transactions on Smart Processing and Computing, 4(1):35–43, 2015.
  * [11] Serkan Kiranyaz, Turker Ince, Osama Abdeljaber, Onur Avci, and Moncef Gabbouj.  1-d convolutional neural networks for signal processing applications.  In ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 8360–8364. IEEE, 2019.
  * [12] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep Learning.  MIT Press, 2016.  http://www.deeplearningbook.org.
  * [13] Michel Verleysen and Damien François.  The curse of dimensionality in data mining and time series prediction.  In International work-conference on artificial neural networks, pages 758–770. Springer, 2005.
  * [14] Florian Schroff, Dmitry Kalenichenko, and James Philbin.  Facenet: A unified embedding for face recognition and clustering.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 815–823, 2015.
  * [15] Andrew L Maas, Awni Y Hannun, Andrew Y Ng, et al.  Rectifier nonlinearities improve neural network acoustic models.  In Proc. icml, volume 30, page 3. Citeseer, 2013.
  * [16] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In International conference on machine learning, pages 448–456. PMLR, 2015.
  * [17] Hao Yan, Mingxin Yu, Jiabin Xia, Lianqing Zhu, Tao Zhang, and Zhihui Zhu.  Tongue squamous cell carcinoma discrimination with raman spectroscopy and convolutional neural networks.  Vibrational Spectroscopy, 103:102938, 2019.
  * [18] William John Thrift, Antony Cabuslay, Andrew Benjamin Laird, Saba Ranjbar, Allon I Hochbaum, and Regina Ragan.  Surface-enhanced raman scattering-based odor compass: Locating multiple chemical sources and pathogens.  ACS sensors, 4(9):2311–2319, 2019.
  * [19] Data and code on serum raman spectroscopy as an efficient primary screening of coronavirus disease in 2019 (covid-19).  https://springernature.figshare.com/articles/dataset/Data_and_code_on_serum_Raman_spectroscopy_as_an_efficient_primary_screening_of_coronavirus_disease_in_2019_COVID-19_/12159924/1, (Last accessed on July 2021).
  * [20] Gang Yin, Lintao Li, Shun Lu, Yu Yin, Yuanzhang Su, Yilan Zeng, Mei Luo, Maohua Ma, Hongyan Zhou, Lucia Orlandini, et al.  An efficient primary screening of covid-19 by serum raman spectroscopy.  Journal of Raman Spectroscopy, 52(5):949–958, 2021.
  * [21] Barbara Lafuente, Robert T Downs, Hexiong Yang, and Nate Stone.  1\. the power of databases: The rruff project.  In Highlights in mineralogical crystallography, pages 1–30. De Gruyter (O), 2015.
  * [22] Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et al.  Tensorflow: A system for large-scale machine learning.  In 12th $\\{$USENIX$\\}$ symposium on operating systems design and implementation ($\\{$OSDI$\\}$ 16), pages 265–283, 2016.
  * [23] Quanshi Zhang and Song-Chun Zhu.  Visual interpretability for deep learning: a survey.  arXiv preprint arXiv:1802.00614, 2018.
  * [24] Nabil Ibtehaz, Muhammad EH Chowdhury, Amith Khandakar, Serkan Kiranyaz, M Sohel Rahman, Anas Tahir, Yazan Qiblawey, and Tawsifur Rahman.  Edith: Ecg biometrics aided by deep learning for reliable individual authentication.  arXiv preprint arXiv:2102.08026, 2021.
  * [25] Tawsifur Rahman, Amith Khandakar, Muhammad Abdul Kadir, Khandaker Rejaul Islam, Khandakar F Islam, Rashid Mazhar, Tahir Hamid, Mohammad Tariqul Islam, Saad Kashem, Zaid Bin Mahbub, et al.  Reliable tuberculosis detection using chest x-ray with deep learning, segmentation and visualization.  IEEE Access, 8:191586–191601, 2020.
  * [26] Scott M Lundberg and Su-In Lee.  A unified approach to interpreting model predictions.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems 30, pages 4765–4774. Curran Associates, Inc., 2017.
  * [27] Slundberg/shap: A game theoretic approach to explain the output of any machine learning model.  https://github.com/slundberg/shap.  Accessed: November 7, 2021.