# A time-stepping deep gradient flow method for option pricing∑in (rough)
diffusion models

Antonis Papapantoleon  and  Jasper Rou

###### Abstract.

We develop a novel deep learning approach for pricing European options in
diffusion models, that can efficiently handle high-dimensional problems
resulting from Markovian approximations of rough volatility models. The option
pricing partial differential equation is reformulated as an energy
minimization problem, which is approximated in a time-stepping fashion by deep
artificial neural networks. The proposed scheme respects the asymptotic
behavior of option prices for large levels of moneyness, and adheres to a
priori known bounds for option prices. The accuracy and efficiency of the
proposed method is assessed in a series of numerical examples, with particular
focus in the lifted Heston model.

###### Key words and phrases:

Option pricing, diffusion models, rough volatility models, lifted volatility
models, PDE, time-stepping, gradient flow, artificial neural network.

###### 2020 Mathematics Subject Classification:

91G20, 91G60, 68T07.

The authors are grateful to Emmanuil Georgoulis and Costas Smaragdakis for
several helpful discussions and advice during the work on this project. AP
gratefully acknowledges the financial support from the Hellenic Foundation for
Research and Innovation (H.F.R.I.) under the “First Call for H.F.R.I. Research
Projects to support Faculty members and Researchers and the procurement of
high-cost research equipment grant” (Project Number: 2152). The authors also
acknowledge the use of computational resources of the DelftBlue supercomputer,
provided by the Delft High-Performance Computing Centre [19].

##  1\. Introduction

Stochastic volatility models have been popular in the mathematical finance
literature because they allow to accurately model and reproduce the shape of
implied volatility smiles for a single maturity. They require though certain
modifications, such as making the parameters time- or maturity-dependent, in
order to reproduce a whole volatility surface; see e.g. the comprehensive
books by Gatheral [25] or Bergomi [15]. The class of rough volatility models,
in which the volatility process is driven by a fractional Brownian motion,
offers an attractive alternative to classical volatility models, since they
allow to reproduce many stylized facts of asset and option prices with only a
few (constant) parameters; see e.g. the seminal articles by Gatheral et al.
[27] and Bayer et al. [9], and the recent volume by Bayer et al. [13]. The
presence of fractional Brownian motion in the dynamics of rough volatility
models yields that the volatility process is not a semimartingale, and the
model is not Markovian, which means that the simulation of the dynamics, or
the pricing and hedging of derivatives, become challenging tasks. This has
attracted increasing attention from the mathematical finance community and has
led to the creation of several innovative methods to tackle these challenges.

Markovian approximations of fractional processes is a popular method for
tackling the challenges arising in these models, since they allow to relate
rough volatility models to their classical counterparts. Abi Jaber and El Euch
[3, 2] have designed multifactor stochastic volatility models approximating
rough volatility models and enjoying a Markovian structure. Zhu et al. [48]
have showed that the rough Bergomi model can be well-approximated by the
forward-variance Bergomi model with appropriately chosen weights and mean-
reversion speed parameters, which has the Markovian property. Bayer and
Breneis [7, 8] have studied Markovian approximations of stochastic Volterra
equations using an $N$-dimensional diffusion process defined as the solution
to a system of (ordinary) stochastic differential equations. Cuchiero and
Teichmann [18] have provided existence, uniqueness and approximation results
for Markovian lifts of affine rough volatility models of general (jump)
diffusion type. Moreover, Bonesini et al. [17] have proposed a theoretical
framework that exploits convolution kernels to transform a Volterra path-
dependent (non-Markovian) stochastic process into a standard (Markovian)
diffusion process.

Novel simulation schemes have also been developed in order to speed up Monte
Carlo methods for rough volatility models. Bennedsen et al. [14] have
developed a hybrid simulation scheme for Brownian semistationary processes,
and applied this to the rough Bergomi model. This method was further refined
by McCrickerd and Pakkanen [41] who developed variance reduction techniques
for conditionally log-normal models, and by Fukasawa and Hirano [24] who
optimized the use of random numbers through orthogonal projections. Moreover,
Gatheral [26] has developed hybrid quadratic-exponential schemes for the
simulation of rough affine forward variance models, while Bayer and Breneis
[6] have developed efficient and accurate simulation schemes for the rough
Heston model based on low-dimensional Markovian representations.

Deep learning and neural network methods have also been applied in the field
of rough stochastic volatility modeling. Horvath et al. [35] presented a
neural network-based method that performs the calibration task for the full
implied volatility surface in (rough) volatility models. Jacquier and Oumgari
[36] introduced a deep learning–based algorithm to evaluate options in affine
rough stochastic volatility models. Jacquier and Zuric [38] constructed a deep
learning-based numerical algorithm to solve path-dependent partial
differential equations (PDEs) arising in the context of rough volatility.

Many other methods have been utilized and applied to rough stochastic
volatility models. Let us mention here, among others, that El Euch and
Rosenbaum [22] have computed the characteristic function of the log‐price in
the rough Heston model, thus opening the door for the application of Fourier
transform methods for option pricing, Bayer et al. [11] proposed the
application of adaptive sparse grids and quasi Monte Carlo methods for option
pricing in the rough Bergomi model, while Bayer et al. [12] have applied
Hairer’s theory of regularity structures in order to analyze rough volatility
models. In addition, Bayer et al. [10], Horvath et al. [34] and Jacquier and
Pannier [37], among others, have studied asymptotic properties of rough
volatility models in various regimes.

The connection between PDEs and option pricing is deep and well-studied in the
literature. Several methods have been developed in order to price options in
diffusion or stochastic volatility models, especially in high-dimensional
settings, using sparse grids, asymptotic expansions, finite difference and
finite elements methods; see e.g. Reisinger and Wittum [44], Düring and
Fournié [20], Hilber et al. [32] and the comprehensive book by Hilber et al.
[33].

The recent advances in deep learning methods have also created an intense
interest in the solution of PDEs by deep learning and neural network methods.
Sirignano and Spiliopoulos [46] proposed to solve high-dimensional PDEs by
approximating the solution with a deep neural network which is trained to
satisfy the differential operator, initial condition, and boundary conditions.
Raissi et al. [43] introduced physics informed neural networks, i.e. neural
networks that are trained to solve supervised learning tasks while respecting
any given law of physics described by general nonlinear PDEs. Van der Meer et
al. [47] discussed the choice of a norm in the loss function for the training
of neural networks. Han et al. [30] introduced a deep learning-based approach
that can handle general high-dimensional parabolic PDEs by reformulating the
PDE using backward stochastic differential equations and approximating the
gradient of the unknown solution by neural networks. E and Yu [21] proposed a
deep learning-based method for numerically solving variational problems,
particularly the ones that arise from PDEs. Liao and Ming [40] proposed a
method to deal with the essential boundary conditions encountered in the deep
learning-based numerical solvers for PDEs. Georgoulis et al. [28] considered
the approximation of initial and boundary value problems involving high-
dimensional, dissipative evolution PDEs using a deep gradient flow neural
network framework. In similar spirit, Park et al. [42] proposed a deep
learning-based minimizing movement scheme for approximating the solutions of
PDEs. The work of [28] was followed up by Georgoulis et al. [29], who
developed an implicit-explicit deep minimizing movement method for pricing
European basket options written on assets that follow jump-diffusion dynamics.
Let us refer the reader to the forthcoming book by Jentzen et al. [39], who
provide an introduction to the topic of deep learning algorithms and a survey
of the (vast) literature.

The starting point for our research is the existence of a Markovian
approximation for a rough volatility model or, more generally, another
Markovian model that describes empirical data well; see e.g. Rømer [45] and
Abi Jaber and Li [4]. Inspired by Georgoulis et al. [28], we consider the
option pricing PDE, introduce a time-stepping discretization of the time
derivative, and reformulate the spatial derivatives as an energy minimization
problem; the main difference with Georgoulis et al. [28] is the presence of
the drift (first order) term, which we treat in an explicit fashion in the
discretization. This gradient flow formulation of the PDE gives rise to a
“natural” cost function for the training, via stochastic gradient descent, of
deep neural networks that approximate the solution of the PDE. The main
advantages of this method are twofold: on the one hand, the neural network
from the previous time step is a good initial guess for the neural network
optimization in the current time step, thus reducing the number of training
stages necessary. On the other hand, the energy formulation allows to consider
only first-order derivatives, thus reducing the training time in each time
step / training stage significantly. In addition, we utilize available
information about the qualitative behavior of option prices in order to assist
the training of the neural network. More specifically, the architecture of the
neural network respects the asymptotic behavior of option prices for large
levels of moneyness, and adheres to a priori known bounds for option prices.
The empirical results show that the proposed method is accurate, with a small
overall error in the order of $10^{-3}$ and comparable with other deep
learning methods, while the training times are significantly small compared to
other methods and do not grow significantly with the dimension of the
underlying space.

This article is organized as follows: in Section 2, we revisit the option
pricing PDE for Markovian diffusion models. In Section 3, we explain how we
develop our method for solving option pricing PDEs, with focus on the time
discretization and energy formulation. In Section 4, we revisit classical
models for asset prices, with particular focus in the lifted Heston model. In
Section 5, we explain certain details for the design and implementation of the
method. In Section 6, we present and discuss the outcome of our numerical
experiments, while Section 7 concludes this article.

##  2\. Problem formulation

Let $S$ denote the price process of a financial asset that evolves according
to a (Markovian) diffusion model, and consider a European vanilla derivative
on $S$ with payoff $\Psi(S_{T})$ at maturity time $T>0$. Using the fundamental
theorem of asset pricing and the Feynman–Kac formula, the price of this
derivative can be written as the solution to a PDE in these models. Indeed,
let $u:[0,T]\times\Omega\to\mathbb{R}$ denote the price of this derivative,
where $\Omega\subseteq\mathbb{R}^{n+1}$ and $t$ denotes the time to maturity.
Then, $u$ solves the general PDE

$\displaystyle\begin{cases}u_{t}(t,x)+\mathcal{A}u(t,x)+ru(t,x)=0,\quad&(t,x)\in[0,T]\times\Omega,\\\
u(0,x)=\Psi(x),\quad&x\in\Omega,\end{cases}$ (2.1)

where $\mathcal{A}$ is a second order differential operator of the form

$\mathcal{A}u=-\sum_{i,j=0}^{n}a^{ij}\frac{\partial^{2}u}{\partial
x_{i}\partial x_{j}}+\sum_{i=0}^{n}\beta^{i}\frac{\partial u}{\partial
x_{i}}.$ (2.2)

The coefficients $a^{ij},\beta^{i}$ of the generator $\mathcal{A}$ are
directly related to the dynamics of the stochastic process $S$ and can, in
general, depend on the time and the spatial variables. Their particular form
will be clarified in the applications later; see Section 4. Moreover, for the
time being, let us also prescribe homogeneous Dirichlet boundary conditions on
$\partial\Omega$.

We will later rewrite the PDE (2.1) as an energy minimization problem,
therefore we need to split the operator in a symmetric and an asymmetric part.
The symmetric part will then be associated to a Dirichlet energy. Using the
product rule, we can rewrite $\mathcal{A}$ as follows

$\displaystyle\mathcal{A}u$
$\displaystyle=\sum_{i=0}^{n}\beta^{i}\frac{\partial u}{\partial
x_{i}}-\sum_{i,j=0}^{n}\frac{\partial}{\partial
x_{j}}\left(a^{ij}\frac{\partial u}{\partial
x_{i}}\right)+\sum_{i,j=0}^{n}\frac{\partial a^{ij}}{\partial
x_{j}}\frac{\partial u}{\partial x_{i}}$
$\displaystyle=\sum_{i=0}^{n}\left(\beta^{i}+\sum_{j=0}^{n}\frac{\partial
a^{ij}}{\partial x_{j}}\right)\frac{\partial u}{\partial
x_{i}}-\sum_{i,j=0}^{n}\frac{\partial}{\partial
x_{j}}\left(a^{ij}\frac{\partial u}{\partial x_{i}}\right).$

Define $b^{i}=\beta^{i}+\sum_{j=0}^{n}\frac{\partial a^{ij}}{\partial x_{j}}$.
Then, $\mathcal{A}$ takes the form

$\mathcal{A}u=-\nabla\cdot\left(A\nabla u\right)+\mathbf{b}\cdot\nabla u,$
(2.3)

where

$A=\begin{bmatrix}a^{00}&a^{10}&\dots&a^{n0}\\\ a^{01}&a^{11}&\dots&a^{n1}\\\
\vdots&\vdots&\ddots&\vdots\\\
a^{0n}&a^{1n}&\dots&a^{nn}\end{bmatrix}\quad\textrm{and}\quad\mathbf{b}=\begin{bmatrix}b^{0}\\\
b^{1}\\\ \vdots\\\ b^{n}\end{bmatrix}.$ (2.4)

##  3\. Energy minimization and a deep gradient flow method

Let us now explain how we will develop our deep neural network method for
solving the PDE (2.1). We proceed in three steps. First, we discretize the
time derivative in Subsection 3.1 and introduce a time-stepping scheme.
Second, we rewrite the solution of the PDE as an energy minimization problem
in Subsection 3.2, which gives rise to a “natural” cost function. Third, we
approximate the minimizer by a neural network, which is trained using
stochastic gradient descent, in Subsection 3.3

###  3.1. Time discretization

A key ingredient of our method is that we do not use time as an input variable
for the neural network, thus training a global space-time network for solving
the PDE, as in Sirignano and Spiliopoulos [46] and Raissi et al. [43].
Instead, following Georgoulis et al. [28], we wish to train a neural network
for each time step. The main insight is that the neural network from the
previous time step is a good initial guess for the neural network optimization
in the current time step, thus reducing the number of training stages
necessary.

Let us divide the time interval $[0,T]$ into $N$ equally spaced intervals
$(t_{k-1},t_{k}]$, with $h=t_{k}-t_{k-1}=\frac{1}{N}$ for $k=0,1,\dots,N$. Let
$U^{k}$ denote the approximation to the solution of the PDE $u(t_{k})$ at time
step $t_{k}$, using the following time discretization scheme

$\displaystyle\frac{U^{k}-U^{k-1}}{h}-\nabla\cdot\left(A\nabla
U^{k}\right)+\mathbf{b}\cdot\nabla U^{k}+rU^{k}$ $\displaystyle=0,$
$\displaystyle U^{0}$ $\displaystyle=\Psi.$

In order to be able to rewrite the PDE as an energy minimization problem, we
wish to treat the asymmetric part as a constant function. Therefore, we
substitute the function $U^{k}$ in the asymmetric part with its value from the
previous time point, and thus consider the backward differentiation scheme

$\displaystyle\frac{U^{k}-U^{k-1}}{h}-\nabla\cdot\left(A\nabla
U^{k}\right)+F\left(U^{k-1}\right)+rU^{k}$ $\displaystyle=0,$ (3.1)
$\displaystyle U^{0}$ $\displaystyle=\Psi,$

where $F(u)=\mathbf{b}\cdot\nabla u$. Convergence analysis of this scheme
appears in Akrivis and Smyrlis [5].

###  3.2. Variational formulation

We would like now to associate the discretized PDE in (3.1) with an energy
functional, in order to define later a “natural” cost function for the
training of the neural network. In other words, we would like to find an
energy functional $I(u)$ such that $U^{k}$ is a critical point of $I$. Let us
rewrite (3.1) as follows:

$\left(U^{k}-U^{k-1}\right)+h\left(-\nabla\cdot\left(A\nabla
U^{k}\right)+rU^{k}+F\left(U^{k-1}\right)\right)=0,\quad U^{0}=\Psi.$ (3.2)

Then, we have the following result.

######  Lemma 3.1.

$U^{k}$ solves (3.2) if and only if $U^{k}$ minimizes

$\displaystyle I^{k}(u)=$
$\displaystyle\frac{1}{2}\left\|u-U^{k-1}\right\|_{L^{2}(\Omega)}^{2}+h\left(\int_{\Omega}\frac{1}{2}\left(\left(\nabla
u\right)^{\mathsf{T}}A\nabla
u+ru^{2}\right)+F\left(U^{k-1}\right)u\mathrm{d}x\right).$

###### Proof.

Let $v$ be a smooth function. Consider the function

$\displaystyle i^{k}(\tau)$ $\displaystyle=I^{k}\left(U^{k}+\tau
v\right)=\int_{\Omega}\frac{1}{2}\left(U^{k}+\tau
v-U^{k-1}\right)^{2}\mathrm{d}x$
$\displaystyle\quad+h\int_{\Omega}\frac{1}{2}\left(\left(\nabla\left(U^{k}+\tau
v\right)\right)^{\mathsf{T}}A\nabla\left(U^{k}+\tau v\right)+r\left(U^{k}+\tau
v\right)^{2}\right)+F\left(U^{k-1}\right)\left(U^{k}+\tau
v\right)\mathrm{d}x,$

for $\tau\in\mathbb{R}$. $U^{k}$ minimizes $I^{k}$ if and only if $\tau=0$
minimizes $i^{k}$. Therefore, we get

$\displaystyle 0$
$\displaystyle=(i^{k})^{\prime}(0)=\left[\int_{\Omega}\left(U^{k}+\tau
v-U^{k-1}\right)v\mathrm{d}x\right.$
$\displaystyle\quad+\left.h\int_{\Omega}\frac{1}{2}\left(\left(\nabla
v\right)^{\mathsf{T}}A\nabla\left(U^{k}+\tau v\right)+\left(\nabla U^{k}+\tau
v\right)^{\mathsf{T}}A\nabla v\right)+r\left(U^{k}+\tau
v\right)v+F\left(U^{k-1}\right)v\mathrm{d}x\right]_{\tau=0}$
$\displaystyle=\int_{\Omega}\left(U^{k}-U^{k-1}\right)v+\frac{h}{2}\left(\left(\nabla
v\right)^{\mathsf{T}}A\nabla U^{k}+\left(\nabla
U^{k}\right)^{\mathsf{T}}A\nabla
v\right)+hrU^{k}v+hF\left(U^{k-1}\right)v\mathrm{d}x$
$\displaystyle=\int_{\Omega}\left(U^{k}-U^{k-1}\right)v+h\left(A\nabla
U^{k}\cdot\nabla v+rU^{k}+F\left(U^{k-1}\right)\right)v\mathrm{d}x$
$\displaystyle=\int_{\Omega}\left(\left(U^{k}-U^{k-1}\right)+h\left(-\nabla\cdot\left(A\nabla
U^{k}\right)+rU^{k}+F\left(U^{k-1}\right)\right)\right)v\mathrm{d}x,$

where in the second to last equality we used that $A$ is symmetric, and in the
last equality we used the divergence theorem. This equality holds for all $v$
if and only if (3.2) holds. ∎

###  3.3. Neural network approximation

Using the results of the previous subsection, we define the following cost
(loss) function for training the neural network:

$\displaystyle
L^{k}(u)=\frac{1}{2}\left\|u-U^{k-1}\right\|_{L^{2}(\Omega)}^{2}+h\int_{\Omega}\mathcal{E}[u]\mathrm{d}x+h\int_{\Omega}F\left(U^{k-1}\right)u\mathrm{d}x,$

where

$\displaystyle\mathcal{E}[u]=\frac{1}{2}\left(\left(\nabla
u\right)^{\mathsf{T}}A\nabla u+ru^{2}\right),$

denotes the respective Dirichlet energy.

Let $f^{k}(\mathbf{x};\theta)$ denote a neural network approximation of
$U^{k}$ with trainable parameters $\theta$. Applying a Monte Carlo
approximation to the integrals, the discretized cost functional takes the form

$\displaystyle L^{k}\left(\theta;\mathbf{x}\right)$
$\displaystyle=\frac{|\Omega|}{2M}\sum_{i=1}^{M}\left(f^{k}(\mathbf{x}_{i};\theta)-f^{k-1}(\mathbf{x}_{i})\right)^{2}+hN^{k}\left(\theta;\mathbf{x}\right),$

where

$\displaystyle N^{k}\left(\theta;\mathbf{x}\right)$
$\displaystyle=\frac{|\Omega|}{M}\sum_{i=1}^{M}\left[\frac{1}{2}\left(\left(\nabla
f^{k}(\mathbf{x}_{i};\theta)\right)^{\mathsf{T}}A\nabla
f^{k}(\mathbf{x}_{i};\theta)+r\left(f^{k}({\mathbf{x}_{i}};\theta)\right)^{2}\right)+\left(\mathbf{b}\cdot\nabla
f^{k-1}(\mathbf{x}_{i})\right)f^{k}(\mathbf{x}_{i};\theta)\right].$

Here, $M$ denotes the number of samples $\mathbf{x}_{i}$.

In order to minimize this cost function, we use the stochastic gradient
descent algorithm, i.e. the iterative scheme:

$\theta_{n+1}=\theta_{n}-\alpha_{n}\nabla_{\theta}L^{k}(\theta_{n};\mathbf{x}).$

The hyper-parameter $\alpha_{n}$ is the step-size of our update, called the
learning rate. An overview of the Time Deep Gradient Flow (TDGF) method
appears in Algorithm 1.

Algorithm 1 Time Deep Gradient Flow method

1:Initialize $\theta_{0}^{0}$.

2:Set $f^{0}(\mathbf{x};\theta)=\Psi(\mathbf{x})$.

3:for each time step $k=1,\dots,N$ do

4: Initialize $\theta_{0}^{k}=\theta^{k-1}$.

5: for each sampling stage $n$ do

6: Generate random points $\mathbf{x}_{i}$ for training.

7: Calculate the cost functional $L^{k}(\theta_{n}^{k};\mathbf{x}_{i})$ for
the sampled points.

8: Take a descent step
$\theta_{n+1}^{k}=\theta_{n}^{k}-\alpha_{n}\nabla_{\theta}L^{k}(\theta_{n}^{k};\mathbf{x}_{i})$.

9: end for

10:end for

##  4\. Examples

Let us now outline the different models to which we will apply our method, and
derive the appropriate form of the generator in the option pricing PDE for
each of these models. In particular, we will treat the Black–Scholes model in
Subsection 4.1, the Heston model in Subsection 4.2, and the lifted Heston
model in Subsection 4.3.

###  4.1. Black–Scholes model

In the Black and Scholes [16] model, the dynamics of the stock price $S$
follow a geometric Brownian motion:

$\mathrm{d}S_{t}=rS_{t}\mathrm{d}t+\sigma S_{t}\mathrm{d}W_{t},\quad S_{0}>0,$

with $r,\sigma\in\mathbb{R}_{+}$ the risk free rate and the volatility
respectively. The generator corresponding to these dynamics, in the form
(2.2), equals

$\mathcal{A}u=-\frac{1}{2}\sigma^{2}x^{2}\frac{\partial^{2}u}{\partial
x^{2}}-rx\frac{\partial u}{\partial x},$

and, applying the chain rule, we have

$\mathcal{A}u=-\frac{\partial}{\partial
x}\left(\frac{1}{2}\sigma^{2}x^{2}\frac{\partial u}{\partial
x}\right)+\sigma^{2}x\frac{\partial u}{\partial x}-rx\frac{\partial
u}{\partial x}.$

Therefore, the operator $\mathcal{A}$ takes the form (2.3) with the
coefficients in (2.4) provided by

$\displaystyle a$ $\displaystyle=\frac{1}{2}\sigma^{2}x^{2},$ $\displaystyle
b$ $\displaystyle=(\sigma^{2}-r)x.$

###  4.2. Heston model

The Heston [31] model is a popular stochastic volatility model with dynamics

$\displaystyle\mathrm{d}S_{t}$
$\displaystyle=rS_{t}\mathrm{d}t+\sqrt{V_{t}}S_{t}\mathrm{d}W_{t},\quad$
$\displaystyle S_{0}>0,$ $\displaystyle\mathrm{d}V_{t}$
$\displaystyle=\lambda(\kappa-
V_{t})\mathrm{d}t+\eta\sqrt{V_{t}}\mathrm{d}B_{t},\quad$ $\displaystyle
V_{0}>0.$

Here $V$ is the variance process, $W,B$ are correlated (standard) Brownian
motions, with correlation coefficient $\rho$, and
$\lambda,\kappa,\eta\in\mathbb{R}_{+}$. The option pricing PDE (2.1) in this
model has two spatial variables $(x,v)$, corresponding to the asset price and
the variance, and the generator corresponding to these dynamics, in the form
(2.2), equals

$\mathcal{A}u=-rx\frac{\partial u}{\partial x}-\lambda(\kappa-v)\frac{\partial
u}{\partial v}-\frac{1}{2}x^{2}v\frac{\partial^{2}u}{\partial
x^{2}}-\frac{1}{2}\eta^{2}v\frac{\partial^{2}u}{\partial v^{2}}-\rho\eta
xv\frac{\partial^{2}u}{\partial x\partial v}.$

Applying the chain rule again, we get

$\displaystyle\mathcal{A}u=$ $\displaystyle-rx\frac{\partial u}{\partial
x}-\lambda(\kappa-v)\frac{\partial u}{\partial v}-\frac{\partial}{\partial
x}\left(\frac{1}{2}x^{2}v\frac{\partial u}{\partial x}\right)+xv\frac{\partial
u}{\partial x}-\frac{\partial}{\partial
v}\left(\frac{1}{2}\eta^{2}v\frac{\partial u}{\partial
v}\right)+\frac{1}{2}\eta^{2}\frac{\partial u}{\partial v}$
$\displaystyle-\frac{\partial}{\partial x}\left(\frac{1}{2}\rho\eta
xv\frac{\partial u}{\partial v}\right)+\frac{1}{2}\rho\eta v\frac{\partial
u}{\partial v}-\frac{\partial}{\partial v}\left(\frac{1}{2}\rho\eta
xv\frac{\partial u}{\partial x}\right)+\frac{1}{2}\rho\eta x\frac{\partial
u}{\partial x}.$

Therefore, the operator $\mathcal{A}$ takes the form (2.3) with the
coefficients in (2.4) provided by

$\displaystyle a^{00}$ $\displaystyle=\frac{1}{2}x^{2}v,$ $\displaystyle
a^{10}=a^{01}$ $\displaystyle=\frac{1}{2}\rho\eta xv,$ $\displaystyle a^{11}$
$\displaystyle=\frac{1}{2}\eta^{2}v,$ $\displaystyle b^{0}$
$\displaystyle=\left(-r+v+\frac{1}{2}\rho\eta\right)x,$ $\displaystyle b^{1}$
$\displaystyle=\lambda(v-\kappa)+\frac{1}{2}\eta^{2}+\frac{1}{2}\rho\eta v.$

###  4.3. Lifted Heston model

The Heston model requires difficult modifications, such as making the
parameters time dependent, in order to adequately calibrate the model to
financial market data. The rough Heston model, in which the volatility process
is driven by a fractional Brownian motion, is an alternative where only a few
(constant) parameters are required in order to adequately calibrate the model;
see Gatheral et al. [27] and Bayer et al. [9]. The dynamics of the variance
process in the rough Heston model are described by

$V_{t}=V_{0}+\int_{0}^{t}K(t-s)\left(\lambda\left(\kappa-
V_{s}\right)\mathrm{d}s+\eta\sqrt{V_{s}}\mathrm{d}B_{s}\right),$

where $K(t)=\frac{t^{H-\frac{1}{2}}}{\Gamma\left(H+\frac{1}{2}\right)}$ is the
fractional kernel and $H\in\left(0,\frac{1}{2}\right)$ is the Hurst index.
This model is not Markovian, due to the presence of the fractional Brownian
motion in the dynamics, hence the pricing and hedging of derivatives becomes a
challenging task. In particular, we do not have access to a PDE of the form
(2.1) for the rough Heston model.

The lifted Heston model is a Markovian approximation of the rough Heston
model, see Abi Jaber and El Euch [3], that is interesting as a model in its
own right; see also Abi Jaber [1]. The dynamics of the asset price in the
lifted Heston model are the following:

$\displaystyle\mathrm{d}S_{t}$
$\displaystyle=rS_{t}\mathrm{d}t+\sqrt{V_{t}^{n}}S_{t}\mathrm{d}W_{t},\quad
S_{0}>0,$ $\displaystyle\mathrm{d}V_{t}^{n,i}$
$\displaystyle=-\left(\gamma_{i}^{n}V_{t}^{n,i}+\lambda
V_{t}^{n}\right)\mathrm{d}t+\eta\sqrt{V_{t}^{n}}\mathrm{d}B_{t},\quad
V_{0}^{n,i}=0,$ (4.1) $\displaystyle V_{t}^{n}$
$\displaystyle=g^{n}(t)+\sum_{i=1}^{n}c_{i}^{n}V_{t}^{n,i},$ (4.2)
$\displaystyle g^{n}(t)$
$\displaystyle=V_{0}+\lambda\kappa\sum_{i=1}^{n}c_{i}^{n}\int_{0}^{t}\mathrm{e}^{-\gamma_{i}^{n}(t-s)}\mathrm{d}s,$

where $\lambda,\eta,c_{i}^{n},\gamma_{i}^{n},V_{0},\kappa\in\mathbb{R}_{+}$,
and $W,B$ are two correlated (standard) Brownian motions with correlation
coefficient $\rho$. The variance factors $V^{n,i}$ start from zero and share
the same one-dimensional Brownian motion $B$.

The option pricing PDE (2.1) has $n+1$ spatial variables
$(x,v_{1},\dots,v_{n})$ in the lifted Heston model, corresponding to the asset
price $S$ and the $n$ variance processes $V^{n,i}$ in (4.1). The next section
explains the sampling of the process $V^{n}$ in (4.2).

Let us now derive the generator for this system of SDEs. We will denote by
$V_{i}=V^{n,i}$, $f_{S}=\frac{\partial f}{\partial S}$ and
$f_{V_{i}}=\frac{\partial f}{\partial V^{n,i}}$, and suppress the dependence
on time for simplicity. Let
$f\left(S,V_{1},\dots,V_{n}\right):\mathbb{R}^{n+1}\to\mathbb{R}$ be a
$\mathcal{C}^{2}$-function. Using Itô’s formula, we have that

$\displaystyle\mathrm{d}f\left(S,V_{1},\dots,V_{n}\right)$
$\displaystyle=f_{S}\mathrm{d}S+\sum_{i=1}^{n}f_{V_{i}}\mathrm{d}V_{i}+\frac{1}{2}f_{SS}\mathrm{d}\left\langle
S,S\right\rangle+\sum_{i=1}^{n}f_{SV_{i}}\mathrm{d}\left\langle
S,V_{i}\right\rangle+\frac{1}{2}\sum_{i,j=1}^{n}f_{V_{i}V_{j}}\mathrm{d}\left\langle
V_{i},V_{j}\right\rangle$
$\displaystyle=f_{S}rS\mathrm{d}t-\sum_{i=1}^{n}f_{V_{i}}\left(\gamma_{i}^{n}V_{i}+\lambda
V_{t}^{n}\right)\mathrm{d}t+\frac{1}{2}f_{SS}V_{t}^{n}\left(S\right)^{2}\mathrm{d}t+\sum_{i=1}^{n}f_{SV_{i}}\eta\rho
V_{t}^{n}S\mathrm{d}t$
$\displaystyle\quad+\frac{1}{2}\sum_{i,j=1}^{n}f_{V_{i}V_{j}}\eta^{2}V_{t}^{n}\mathrm{d}t+\textrm{(local)
martingale}.$

Hence, the generator $\tilde{\mathcal{A}}$ of this system takes the form:

$\tilde{\mathcal{A}}f=-rxf_{x}+\sum_{i=1}^{n}\left(\gamma_{i}^{n}v_{i}+\lambda
V_{t}^{n}\right)f_{v_{i}}-\frac{1}{2}V_{t}^{n}x^{2}f_{xx}-\eta\rho
V_{t}^{n}x\sum_{i=1}^{n}f_{xv_{i}}-\frac{\eta^{2}}{2}V_{t}^{n}\sum_{i,j=1}^{n}f_{v_{i}v_{j}}.$

In order to bring this generator in the form (2.2), we need to switch to the
time to maturity, denoted now $\tau=T-t$. We have to be particularly careful
in this case, because the function $g^{n}$ in $V^{n}$ depends on the time $t$.
Let us denote

$\tilde{V}_{\tau}^{n}:=g^{n}(T-\tau)+\sum_{i=1}^{n}c_{i}^{n}V_{\tau}^{n,i}.$

Then, applying the chain rule again, we have that

$\displaystyle\mathcal{A}u$
$\displaystyle=-rxu_{x}+\sum_{i=1}^{n}\left(\gamma_{i}^{n}v_{i}+\lambda\tilde{V}_{\tau}^{n}\right)u_{v_{i}}-\left(\frac{1}{2}\tilde{V}_{\tau}^{n}x^{2}u_{x}\right)_{x}+\tilde{V}_{\tau}^{n}xu_{x}-\frac{1}{2}\sum_{i=1}^{n}\left(\eta\rho\tilde{V}_{\tau}^{n}xu_{v_{i}}\right)_{x}$
$\displaystyle\quad+\frac{1}{2}\eta\rho\tilde{V}_{\tau}^{n}\sum_{i=1}^{n}u_{v_{i}}-\frac{1}{2}\sum_{i=1}^{n}\left(\eta\rho\tilde{V}_{\tau}^{n}xu_{x}\right)_{v_{i}}+\frac{1}{2}\eta\rho
x\sum_{i=1}^{n}c_{i}^{n}u_{x}-\sum_{i,j=1}^{n}\left(\frac{\eta^{2}}{2}\tilde{V}_{\tau}^{n}u_{v_{i}}\right)_{v_{j}}$
$\displaystyle\quad+\frac{\eta^{2}}{2}\sum_{i,j=1}^{n}c_{j}^{n}u_{v_{i}}$
$\displaystyle=\left(\tilde{V}_{\tau}^{n}x-rx+\frac{1}{2}\eta\rho
x\sum_{i=1}^{n}c_{i}^{n}\right)u_{x}+\sum_{i=1}^{n}\left(\gamma_{i}^{n}v_{i}+\lambda\tilde{V}_{\tau}^{n}+\frac{1}{2}\eta\rho\tilde{V}_{\tau}^{n}+\frac{\eta^{2}}{2}\sum_{j=1}^{n}c_{j}^{n}\right)u_{v_{i}}$
$\displaystyle\quad-\left(\frac{1}{2}\tilde{V}_{\tau}^{n}x^{2}u_{x}\right)_{x}-\frac{1}{2}\sum_{i=1}^{n}\left(\eta\rho\tilde{V}_{\tau}^{n}xu_{v_{i}}\right)_{x}-\frac{1}{2}\sum_{i=1}^{n}\left(\eta\rho\tilde{V}_{\tau}^{n}xu_{x}\right)_{v_{i}}-\sum_{i,j=1}^{n}\left(\frac{\eta^{2}}{2}\tilde{V}_{\tau}^{n}u_{v_{i}}\right)_{v_{j}}.$

Therefore, the operator $\mathcal{A}$ takes the form (2.3) with the
coefficients in (2.4) provided by

$\displaystyle a^{00}$ $\displaystyle=\frac{1}{2}\tilde{V}_{\tau}^{n}x^{2},$
$\displaystyle a^{i0}=a^{0i}$
$\displaystyle=\frac{1}{2}\eta\rho\tilde{V}_{\tau}^{n}x,\quad$ $\displaystyle
i=1,\dots,n,$ $\displaystyle a^{ij}$
$\displaystyle=\frac{\eta^{2}}{2}\tilde{V}_{\tau}^{n},\quad$ $\displaystyle
i,j=1,\dots,n,$ $\displaystyle b^{0}$
$\displaystyle=\left(\tilde{V}_{\tau}^{n}-r+\frac{1}{2}\eta\rho\sum_{i=1}^{n}c_{i}^{n}\right)x,$
$\displaystyle b^{i}$
$\displaystyle=\left(\gamma_{i}^{n}v_{i}+\lambda\tilde{V}_{\tau}^{n}+\frac{1}{2}\eta\rho\tilde{V}_{\tau}^{n}+\frac{\eta^{2}}{2}\sum_{j=1}^{n}c_{j}^{n}\right),\quad$
$\displaystyle i=1,\dots,n.$

##  5\. Implementation details

Let us now describe some details about the design of the neural network
architecture and the implementation of the numerical method. Motivated by
Georgoulis et al. [29], we would like to use information about the option
price in order to facilitate the training of the neural network. On the one
hand, using the lower no-arbitrage bound of a call option, i.e. that
$u(t,x)\geq\left(x-K\mathrm{e}^{-rt}\right)^{+}$, the neural network will
learn the difference between the option price and this bound, instead of the
option price itself. This difference turns out to be an easier function to
approximate than the option price.

On the other hand, after some large level of asset price / moneyness, the
option price will grow linearly with the asset price / moneyness. Let us
denote this level by $x_{p}$. Therefore, instead of letting the neural network
learn the option price at points larger than $x_{p}$, we will approximate this
price by adding the difference between $x_{p}$ and this point to the option
value at $x_{p}$. In other words,

$u(x_{p}+y;\theta)=u(x_{p};\theta)+y-x_{p},\quad y>0.$

This is consistent with the imposition of Dirichlet boundary conditions on
$\Omega$.

Moreover, considering the moneyness $\frac{S}{K}$ instead of the stock price
$S$, will return the relative option price $\frac{u}{K}$. Thus, we can
approximate the option prices for all combinations of $S$ and $K$ with only
one parameter.

The architecture of the neural network for the TDGF method is inspired by the
Deep Galerkin method (DGM) of Sirignano and Spiliopoulos [46], hence we call
the hidden layers “DGM layers”. Overall, we set the following:

$\displaystyle\text{input}\qquad\quad X^{0}$
$\displaystyle=\sigma_{1}\left(W^{1}\mathbf{x}+b^{1}\right),$ DGM layer
$\displaystyle|\quad Z^{l}$
$\displaystyle=\sigma_{1}\left(U^{z,l}\mathbf{x}+W^{z,l}X^{l}+b^{z,l}\right),\quad$
$\displaystyle l=1,\dots,L,$ $\displaystyle|\quad G^{l}$
$\displaystyle=\sigma_{1}\left(U^{g,l}\mathbf{x}+W^{g,l}X^{l}+b^{g,l}\right),\quad$
$\displaystyle l=1,\dots,L,$ $\displaystyle|\quad R^{l}$
$\displaystyle=\sigma_{1}\left(U^{r,l}\mathbf{x}+W^{r,l}X^{l}+b^{r,l}\right),\quad$
$\displaystyle l=1,\dots,L,$ $\displaystyle|\quad\\!H^{l}$
$\displaystyle=\sigma_{1}\left(U^{h,l}\mathbf{x}+W^{h,l}\left(X^{l}\odot
R^{l}\right)+b^{h,l}\right),\quad$ $\displaystyle l=1,\dots,L,$
$\displaystyle\lfloor\quad\\!X^{l}$ $\displaystyle=\left(1-G^{l}\right)\odot
H^{l}+Z^{l}\odot X^{l-1},\quad$ $\displaystyle l=1,\dots,L,$
$\displaystyle\text{output}\quad f(x;\theta)$
$\displaystyle=\left(x-K\mathrm{e}^{-rt}\right)^{+}+\sigma_{2}\left(WX^{L}+b\right).$

Here, $L$ is the number of hidden layers, $\sigma_{i}$ is the activation
function for $i=1,2$, and $\odot$ denotes the element-wise multiplication.
Moreover, let us point here that x denotes the vector of inputs, which are the
asset price and the variances (in the stochastic volatility models), and then
$\mathbf{x}_{1}=x$. In the numerical experiments, we have used 3 layers and 50
neurons per layer. The activation functions we have selected are the
hyperbolic tangent function, $\sigma_{1}(x)=\tanh(x)$, and the softplus
function, $\sigma_{2}(x)=\log\left(\mathrm{e}^{x}+1\right)$, which guarantees
that the option price remains above the no-arbitrage bound.

We consider a maturity of $T=1.0$ year, and take the number of time steps
equal to $N=100$. We use $2000$ sampling stages in each time step, and in each
sampling stage we take $600(n+1)$ samples, following the recommendations in
Georgoulis et al. [28]. As for the sampling domain, we consider the moneyness
$x\in[0.01,3.0]$, the Heston volatility $V\in[0.001,0.1]$ and the lifted
Heston volatilities
$V^{n,i}\in\left[-V_{i}^{\text{high}},V_{i}^{\text{high}}\right]$, with

$V_{i}^{\text{high}}=3\sqrt{\frac{\eta^{2}V_{0}}{2\gamma_{i}^{n}}\left(1-\mathrm{e}^{-2\gamma_{i}^{n}T}\right)}.$

After sampling $V^{n,i}$ uniformly, we calculate $V^{n}$ and discard all
combinations that result in a negative $V^{n}$. We set $x_{p}=2.0$ and use a
learning rate $\alpha=3\times 10^{-4}$. The training was performed on the
DelftBlue supercomputer [19], using a single NVidia Tesla V100S GPU .

##  6\. Numerical results

Let us now present and discuss the outcome of certain numerical experiments
involving all models outlined in Section 4. We are interested in both the
accuracy, in Subsection 6.1, and the speed, in Subsection 6.2, of the
numerical method developed in this article.

We will compare the TDGF method with another popular deep learning method for
solving PDEs, the DGM of Sirignano and Spiliopoulos [46]. In the DGM approach,
the solution of the PDE is translated into a quadratic minimization problem,
which is minimized using stochastic gradient descent. In our setting, this
minimization takes the form:

$\displaystyle\left\|u_{t}-\nabla\cdot\left(A\nabla
u\right)+\mathbf{b}\cdot\nabla
u+ru\right\|_{L^{2}([0,T]\times\Omega)}^{2}+\left\|u(0,\mathbf{x})-\Psi(x)\right\|_{L^{2}(\Omega)}^{2}\to\min.$

Let us recall that, for large levels of asset price / moneyness, the option
price grows linearly with the asset price / moneyness. Therefore, we add a
loss term punishing the derivative of the option price with respect to the
stock deviating from 1 at $x=3.0$.

In order to have a fair comparison between the two methods, we use 200,000
sampling stages for the DGM. In each stage, we use $600(n+2)$ samples for the
PDE in the interior domain and $600(n+1)$ samples for the initial and boundary
conditions. The learning rate follows the schedule from Georgoulis et al. [28,
eq. (24)].

In order to compare the accuracy of the two methods, we need a reference
value. In the Black–Scholes model, in the case of a European call option, the
option pricing PDE has an exact solution:

$u(t,x)=x\Phi(d_{1})-K\mathrm{e}^{-r\tau}\Phi(d_{2}),$

with $\Phi$ the cumulative standard normal distribution function, $\tau=T-t$
the time to maturity,

$d_{1}=\frac{\log\left(\frac{x}{K}\right)+\left(r+\frac{\sigma^{2}}{2}\right)\tau}{\sigma\sqrt{\tau}}\quad\textrm{and}\quad
d_{2}=d_{1}-\sigma\sqrt{\tau}.$

In the Heston model, the option pricing PDE does not admit an analytical
solution. However, the model is affine and the characteristic function does
have an analytical expression; see Heston [31]. Using the characteristic
function, we compute the reference price with the COS method of Fang and
Oosterlee [23], with $N=512$ terms and truncation domain
$\left[-10\tau^{\frac{1}{4}},10\tau^{\frac{1}{4}}\right]$.

Finally, in the lifted Heston model, the characteristic function does not have
an analytical expression but, as this model is again affine, the
characteristic function is known up to the solution of a system of Riccati
equations; see Abi Jaber [1]. We solve these equations using an implicit-
explicit discretization scheme with 500 time steps, after which we apply the
COS method again to compute the reference price.

###  6.1. Accuracy

In the Black–Scholes model, we only have one input variable, the asset price
$S$. Figure 1 presents the relative $L^{2}$\- and the absolute maximum error
for the TDGF and the DGM methods in the Black–Scholes model. These errors are
calculated over 41 evenly spaced grid points in the interval $[0.01,3]$. Both
methods have a comparable (small) error.

(a) Relative $L^{2}$-error.

(b) Maximal absolute error.

Figure 1. Errors of the two methods in the Black–Scholes model against time,
with interest rate $r=0.05$ and volatility $\sigma=0.25$.

In the Heston model, we have two input variables, the asset price $S$ and the
variance $V$. Figure 2 presents the option price against the moneyness level
for different time instances, and Figure 3 presents the relative $L^{2}$\- and
the absolute maximum error for the TDGF and the DGM methods in the Heston
model. These figures are presented for fixed $V_{0}=0.03$, however other
values of $V_{0}$ between 0.01 and 0.09 provide similar results. The DGM is
more accurate than the TDGF method, although the TDGF is just as accurate as
in the Black–Scholes model, and the error is any case small, in the order of
$10^{-3}$.

Figure 2. Relative European call option prices in the Heston model against the
moneyness of the stock, compared to the COS and the DGM methods, for four
different times to maturity with $r=0.0$ and $\eta=0.1$, $\rho=0.0$,
$\kappa=0.01$, $V_{0}=0.03$ and $\lambda=2.0$.

(a) Relative $L^{2}$-error

(b) Maximal absolute error

Figure 3. Error of the two methods in the Heston model against time, with
$r=0.0$ and $\eta=0.1$, $\rho=0.0$, $\kappa=0.01$, $V_{0}=0.03$ and
$\lambda=2.0$.

In the lifted Heston model, we have $n+1$ input variables, the asset price $S$
and the variances $V^{n,i}$, with $i=1,\dots,n$. Figure 4 presents the
relative $L^{2}$\- and the maximum absolute error for the TDGF and the DGM
methods in this model, in the case $n=1$, using the same parameters as in the
Heston model before. The numerical results are comparable for the two methods,
and both methods even seem to perform slightly better than in the previous
example.

(a) Relative $L^{2}$-error

(b) Maximal absolute error

Figure 4. Errors of the two methods in the lifted Heston model with $n=1$
variance process against time, with $r=0.0$, $\eta=0.1$, $\rho=0.0$,
$\kappa=0.01$, $V_{0}=0.01$ and $\lambda=2.0$.

Figure 5 presents the option price against the moneyness level, and Figure 6
presents the relative $L^{2}$\- and the maximum absolute error for the TDGF
and the DGM methods again in the lifted Heston model, for the same parameter
set as previously, but now with $n=20$ variance factors. The errors of both
methods are comparable and, although they have increased compared to the $n=1$
case, they are still low, in the order of $10^{-3}$.

Figure 5. Relative European call option prices in the lifted Heston model with
$n=20$ variance processes against the moneyness of the stock, compared to the
COS and the DGM methods, for four different times to maturity with $r=0.0$,
$\eta=0.1$, $\rho=0.0$, $\kappa=0.01$, $V_{0}=0.01$ and $\lambda=2.0$.

(a) Relative $L^{2}$-error

(b) Maximal absolute error

Figure 6. Errors of the two methods in the lifted Heston model with $n=20$
variance process against time, with $r=0.0$, $\eta=0.1$, $\rho=0.0$,
$\kappa=0.01$, $V_{0}=0.01$ and $\lambda=2.0$.

Figure 7 presents the option price against the moneyness level, and Figure 8
presents the relative $L^{2}$\- and the absolute maximum error for the two
methods in the lifted Heston model, now for the parameter set suggested by Abi
Jaber [1] with $n=1$. The main difference here is that the correlation between
the asset price and the variance processes is no longer zero. The error in
both methods has slightly increased with respect to the previous example, with
the DGM performing slightly better than the TDGF. However, the overall error
is still small, in the order of $10^{-3}$.

Figure 7. Relative European call option prices in the lifted Heston model with
$n=1$ variance process against the moneyness of the stock, compared to the COS
and the DGM methods, for four different times to maturity with $r=0.0$,
$\eta=0.3$, $\rho=-0.7$, $\kappa=0.02$, $V_{0}=0.02$ and $\lambda=0.3$.

(a) Relative $L^{2}$-error

(b) Maximal absolute error

Figure 8. Errors of the two methods in the lifted Heston model with $n=1$
variance process against time, with $r=0.0$, $\eta=0.3$, $\rho=-0.7$,
$\kappa=0.02$, $V_{0}=0.02$ and $\lambda=0.3$.

Finally, Figure 9 presents the option price against the moneyness level, and
Figure 10 presents the relative $L^{2}$\- and the absolute maximum error for
the two methods in the lifted Heston model, for the same parameter set as in
the previous two figures, but now with $n=5$ and $\rho=-0.2$. The combination
of increased dimension (i.e. 5 variance processes compared to one before) and
decreased correlation results in the errors of both methods being similar to
the previous case, and still small, in the order of $10^{-3}$. However, as the
number of variance processes increase while the correlation is non-zero, then
the errors in both methods are increasing rapidly, and we can no longer
accommodate e.g. 20 variance factors as we did in the uncorrelated case.

Figure 9. Relative European call option prices in the lifted Heston model with
$n=5$ variance processes against the moneyness of the stock, compared to the
COS and the DGM methods, for four different times to maturity with $r=0.0$,
$\eta=0.3$, $\rho=-0.2$, $\kappa=0.02$, $V_{0}=0.02$ and $\lambda=0.3$.

(a) Relative $L^{2}$-error

(b) Maximal absolute error

Figure 10. Errors of the two methods in the lifted Heston model with $n=5$
variance processes against time, with $r=0.0$, $\eta=0.3$, $\rho=-0.2$,
$\kappa=0.02$, $V_{0}=0.02$ and $\lambda=0.3$.

###  6.2. Training and computational times

Table 1 summarizes the training times for the TDGF and the DGM methods in the
different models. As we have expected, due to the time stepping and the
absence of a second derivative in the cost function, the training of the TDGF
method is faster than for the DGM method. Moreover, the training time for the
TDGF method does not increase as the dimension of the problem (i.e. the number
of variance processes) increases. On the contrary, the training time for the
DGM method increases significantly with the number of dimensions, and becomes
more than 4 times slower in the lifted Heston model when the dimension of the
variance processes increases from 1 to 20.

The computational times for the TDGF and the DGM methods in all models are
presented in Table 2. The computational times are taken as the average over 34
time points of computing the option prices on a grid of 41 moneyness points.
Even though we take an average, there is still a high variance in the
computational times. Both methods are significantly faster than the COS method
in the lifted Heston model, in which we do not know the characteristic
function explicitly and a system of Riccati equations needs to be solved for
every evaluation of the characteristic function. The DGM method tends to be
slightly faster than the TDGF method, since the latter needs to load the model
corresponding to the time at which we want to compute the option price,
whereas the DGM method only needs to load the weights once for all times.

Model | Black–Scholes | Heston | lifted Heston, $n=1$ | lifted Heston, $n=5$ | lifted Heston, $n=20$
---|---|---|---|---|---
DGM | $7.5\times 10^{3}$ | $12.5\times 10^{3}$ | $13.3\times 10^{3}$ | $22.3\times 10^{3}$ | $56.1\times 10^{3}$
TDGF | $4.1\times 10^{3}$ | $6.0\times 10^{3}$ | $6.4\times 10^{3}$ | $6.7\times 10^{3}$ | $7.6\times 10^{3}$
Table 1. Training time in seconds of the different methods for a European call option in the different models. Model | Black–Scholes | Heston | lifted Heston, $n=1$ | lifted Heston, $n=5$ | lifted Heston, $n=20$
---|---|---|---|---|---
Exact/COS | 0.00025 | 0.018 | 8.9 | 9.2 | 10.4
DGM | 0.0043 | 0.0016 | 0.0034 | 0.027 | 0.0053
TDGF | 0.024 | 0.0060 | 0.020 | 0.010 | 0.025
Table 2. Computational time in seconds of the different methods for a European
call option in the different models.

##  7\. Conclusion

We have developed in this article a novel deep learning method for option
pricing in diffusion models. Starting from the option pricing PDE, using a
backward differentiation time-stepping scheme and results from the calculus of
variations, we can rewrite the PDE as an energy minimization problem for a
suitable energy functional. The time-stepping provides a good initial guess
for the next step in the optimization procedure, while the energy formulation
yields an equation that only has first derivatives, instead of second ones.
The energy minimization problem offers a natural candidate for a cost
function, and the PDE is solved by training a neural network using stochastic
gradient descent. We have considered as examples for our method the
Black–Scholes model, the Heston model and the lifted Heston model, and derived
the appropriate representations for the PDE in all these cases.

Overall, the TDGF method developed here performs well in the numerical
experiments, with an error in the order of $10^{-3}$. Moreover, the training
times for the TDGF method are faster compared to the popular DGM method, while
the computational times for the lifted Heston model are much faster than the
standard COS method. However, when the dimension of the variance processes in
the lifted Heston model increases and the correlation increases
simultaneously, then the error also increases. This increase is probably due
to the hypoellipticity of the equations, and the development of suitable
numerical methods will be the subject of future research.

## References

  * Abi Jaber [2019] E. Abi Jaber.  Lifting the Heston model.  _Quantitative Finance_ , 19:1995–2013, 2019.
  * Abi Jaber and El Euch [2019a] E. Abi Jaber and O. El Euch.  Markovian structure of the Volterra Heston model.  _Statistics & Probability Letters_, 149:63–72, 2019a.
  * Abi Jaber and El Euch [2019b] E. Abi Jaber and O. El Euch.  Multifactor approximation of rough volatility models.  _SIAM Journal on Financial Mathematics_ , 10:309–349, 2019b.
  * Abi Jaber and Li [2024] E. Abi Jaber and S. Li.  Volatility models in practice: Rough, path-dependent or Markovian?  Preprint, arXiv:2401.03345, 2024.
  * Akrivis and Smyrlis [2004] G. Akrivis and Y.-S. Smyrlis.  Implicit–explicit BDF methods for the Kuramoto–Sivashinsky equation.  _Applied numerical mathematics_ , 51(2-3):151–169, 2004.
  * Bayer and Breneis [2023a] C. Bayer and S. Breneis.  Efficient option pricing in the rough Heston model using weak simulation schemes.  Preprint, arXiv:2310.04146, 2023a.
  * Bayer and Breneis [2023b] C. Bayer and S. Breneis.  Markovian approximations of stochastic Volterra equations with the fractional kernel.  _Quantitative Finance_ , 23(1):53–70, 2023b.
  * Bayer and Breneis [2023c] C. Bayer and S. Breneis.  Weak Markovian approximations of rough Heston.  Preprint arXiv:2309.07023, 2023c.
  * Bayer et al. [2016] C. Bayer, P. K. Friz, and J. Gatheral.  Pricing under rough volatility.  _Quantitative Finance_ , 16:887–904, 2016.
  * Bayer et al. [2019] C. Bayer, P. K. Friz, A. Gulisashvili, B. Horvath, and B. Stemper.  Short-time near-the-money skew in rough fractional volatility models.  _Quantitative Finance_ , 19:779–798, 2019.
  * Bayer et al. [2020a] C. Bayer, C. Ben Hammouda, and R. Tempone.  Hierarchical adaptive sparse grids and quasi-Monte Carlo for option pricing under the rough Bergomi model.  _Quantitative Finance_ , 20:1457–1473, 2020a.
  * Bayer et al. [2020b] C. Bayer, P. K. Friz, P. Gassiat, J. Martin, and B. Stemper.  A regularity structure for rough volatility.  _Mathematical Finance_ , 30:782–832, 2020b.
  * Bayer et al. [2023] C. Bayer, P. K. Friz, M. Fukasawa, J. Gatheral, A. Jacquier, and M. Rosenbaum.  _Rough Volatility_.  Society for Industrial and Applied Mathematics, 2023.
  * Bennedsen et al. [2017] M. Bennedsen, A. Lunde, and M. S. Pakkanen.  Hybrid scheme for Brownian semistationary processes.  _Finance and Stochastics_ , 21:931–965, 2017.
  * Bergomi [2016] L. Bergomi.  _Stochastic Volatility Modeling_.  CRC Press, Boca Raton, FL, 2016.
  * Black and Scholes [1973] F. Black and M. Scholes.  The pricing of options and corporate liabilities.  _Journal of Political Economy_ , 81(3):637–654, 1973.
  * Bonesini et al. [2023] O. Bonesini, G. Callegaro, M. Grasselli, and G. Pagès.  From elephant to goldfish (and back): memory in stochastic Volterra processes.  Preprint, arXiv:2306.02708, 2023.
  * Cuchiero and Teichmann [2020] C. Cuchiero and J. Teichmann.  Generalized Feller processes and Markovian lifts of stochastic Volterra processes: the affine case.  _Journal of Evolution Equations_ , 20(4):1301–1348, 2020.
  * Delft High Performance Computing Centre [DHPC] Delft High Performance Computing Centre (DHPC).  DelftBlue Supercomputer (Phase 1).  https://www.tudelft.nl/dhpc/ark:/44463/DelftBluePhase1, 2022.
  * Düring and Fournié [2012] B. Düring and M. Fournié.  High-order compact finite difference scheme for option pricing in stochastic volatility models.  _Journal of Computational and Applied Mathematics_ , 236:4462–4473, 2012.
  * E and Yu [2018] W. E and B. Yu.  The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems.  _Communications in Mathematics and Statistics_ , 6:1–12, 2018.
  * El Euch and Rosenbaum [2019] O. El Euch and M. Rosenbaum.  The characteristic function of rough Heston models.  _Mathematical Finance_ , 29:3–38, 2019.
  * Fang and Oosterlee [2009] F. Fang and C. W. Oosterlee.  A novel pricing method for European options based on Fourier-cosine series expansions.  _SIAM Journal on Scientific Computing_ , 31:826–848, 2009.
  * Fukasawa and Hirano [2021] M. Fukasawa and A. Hirano.  Refinement by reducing and reusing random numbers of the hybrid scheme for Brownian semistationary processes.  _Quantitative Finance_ , 21:1127–1146, 2021.
  * Gatheral [2012] J. Gatheral.  _The Volatility Surface: A Practitioner’s Guide_.  Wiley, 2012.
  * Gatheral [2022] J. Gatheral.  Efficient simulation of affine forward variance models.  _Risk_ , February, 2022.
  * Gatheral et al. [2018] J. Gatheral, T. Jaisson, and M. Rosenbaum.  Volatility is rough.  _Quantitative Finance_ , 18:933–949, 2018.
  * Georgoulis et al. [2023] E. H. Georgoulis, M. Loulakis, and A. Tsiourvas.  Discrete gradient flow approximations of high dimensional evolution partial differential equations via deep neural networks.  _Communications in Nonlinear Science and Numerical Simulation_ , 117:106893, 2023.
  * Georgoulis et al. [2024] E. H. Georgoulis, A. Papapantoleon, and C. Smaragdakis.  A deep implicit-explicit minimizing movement method for option pricing in jump-diffusion models.  Preprint, arXiv:2401.06740, 2024.
  * Han et al. [2018] J. Han, A. Jentzen, and W. E.  Solving high-dimensional partial differential equations using deep learning.  _Proceedings of the National Academy of Sciences_ , 115(34):8505–8510, 2018.
  * Heston [1993] S. L. Heston.  A closed-form solution for options with stochastic volatility with applications to bond and currency options.  _The Review of Financial Studies_ , 6(2):327–343, 1993.
  * Hilber et al. [2005] N. Hilber, A.-M. Matache, and C. Schwab.  Sparse wavelet methods for option pricing under stochastic volatility.  _Journal of Computational Finance_ , 8(4), 2005.
  * Hilber et al. [2013] N. Hilber, O. Reichmann, C. Schwab, and C. Winter.  _Computational Methods for Quantitative Finance: Finite Element Methods for Derivative Pricing_.  Springer Science & Business Media, 2013.
  * Horvath et al. [2019] B. Horvath, A. Jacquier, and C. Lacombe.  Asymptotic behaviour of randomised fractional volatility models.  _Journal of Applied Probability_ , 56:496–523, 2019.
  * Horvath et al. [2021] B. Horvath, A. Muguruza, and M. Tomas.  Deep learning volatility: a deep neural network perspective on pricing and calibration in (rough) volatility models.  _Quantitative Finance_ , 21:11–27, 2021.
  * Jacquier and Oumgari [2023] A. Jacquier and M. Oumgari.  Deep curve-dependent PDEs for affine rough volatility.  _SIAM Journal on Financial Mathematics_ , 14(2):353–382, 2023.
  * Jacquier and Pannier [2022] A. Jacquier and A. Pannier.  Large and moderate deviations for stochastic Volterra systems.  _Stochastic Processes and their Applications_ , 149:142–187, 2022.
  * Jacquier and Zuric [2023] A. Jacquier and Z. Zuric.  Random neural networks for rough volatility.  Preprint, arXiv:2305.01035, 2023.
  * Jentzen et al. [2023] A. Jentzen, B. Kuckuck, and P. von Wurstemberger.  _Mathematical Introduction to Deep Learning: Methods, Implementations, and Theory_.  (forthcoming), 2023.
  * Liao and Ming [2019] Y. Liao and P. Ming.  Deep Nitsche method: Deep Ritz method with essential boundary conditions.  Preprint, arXiv:1912.01309, 2019.
  * McCrickerd and Pakkanen [2018] R. McCrickerd and M. S. Pakkanen.  Turbocharging Monte Carlo pricing for the rough Bergomi model.  _Quantitative Finance_ , 18:1877–1886, 2018.
  * Park et al. [2023] M. S. Park, C. Kim, H. Son, and H. J. Hwang.  The deep minimizing movement scheme.  _Journal of Computational Physics_ , 494:112518, 2023.
  * Raissi et al. [2019] M. Raissi, P. Perdikaris, and G. E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational Physics_ , 378:686–707, 2019.
  * Reisinger and Wittum [2007] C. Reisinger and G. Wittum.  Efficient hierarchical approximation of high-dimensional option pricing problems.  _SIAM Journal on Scientific Computing_ , 29:440–458, 2007.
  * Rømer [2022] S. E. Rømer.  Empirical analysis of rough and classical stochastic volatility models to the SPX and VIX markets.  _Quantitative Finance_ , 22:1805–1838, 2022.
  * Sirignano and Spiliopoulos [2018] J. Sirignano and K. Spiliopoulos.  DGM: A deep learning algorithm for solving partial differential equations.  _Journal of Computational Physics_ , 375:1339–1364, 2018.
  * Van der Meer et al. [2022] R. Van der Meer, C. W. Oosterlee, and A. Borovykh.  Optimally weighted loss functions for solving PDEs with neural networks.  _Journal of Computational and Applied Mathematics_ , 405:113887, 2022.
  * Zhu et al. [2021] Q. Zhu, G. Loeper, W. Chen, and N. Langrené.  Markovian approximation of the rough Bergomi model for Monte Carlo option pricing.  _Mathematics_ , 9(5):528, 2021.


# Improving multi-speaker TTS prosody variance with a residual encoder and
normalizing flows

###### Abstract

Text-to-speech systems recently achieved almost indistinguishable quality from
human speech. However, the prosody of those systems is generally flatter than
natural speech, producing samples with low expressiveness. Disentanglement of
speaker id and prosody is crucial in text-to-speech systems to improve on
naturalness and produce more variable syntheses. This paper proposes a new
neural text-to-speech model that approaches the disentanglement problem by
conditioning a Tacotron2-like architecture on flow-normalized speaker
embeddings, and by substituting the reference encoder with a new learned
latent distribution responsible for modeling the intra-sentence variability
due to the prosody. By removing the reference encoder dependency, the speaker-
leakage problem typically happening in this kind of systems disappears,
producing more distinctive syntheses at inference time. The new model achieves
significantly higher prosody variance than the baseline in a set of
quantitative prosody features, as well as higher speaker distinctiveness,
without decreasing the speaker intelligibility. Finally, we observe that the
normalized speaker embeddings enable much richer speaker interpolations,
substantially improving the distinctiveness of the new interpolated speakers.

Index Terms: text-to-speech, disentanglement, prosody, Tacotron

##  1 Introduction

In the last five years speech technologies have improved considerably. The
text-to-speech (TTS hereafter) field has been largely benefited by the rise of
deep learning [1], which allowed these systems to achieve near-human
performance at synthesizing speech that is almost indistinguishable from
human’s. One of the first achievements was WaveNet [2], a neural vocoder based
on dilated causal convolutions that was able to surpass its predecessors in
naturalness in 2016. With the rise of the attention mechanism and its variants
[3, 4, 5], Tacotron [6] and Tacotron 2 [7, 8] proposed a sequence to sequence
architecture as an end-to-end TTS solution. These models were able to map
input text (or phonemes) to a spectrogram that, given the right neural vocoder
(WaveNet for instance), would be converted into a sequence of waveform
samples.

TTS is fundamentally a generative modeling problem because a given sentence
can be mapped to multiple utterances with different prosody and speaker
characteristics [9]. Tacotron, in its initial form, is a supervised model that
performs a hard mapping between the input text and the output spectrograms.
Therefore, the utterances generated using Tacotron-like architectures,
although natural, they follow the average prosody of the training set, not
allowing to generate utterances of a sentence with multiple speaking styles.

The authors of [10] attempt to remediate the lack of expressiveness of the
model by introducing a reference encoder consisting of a latent distribution
that is conditioned on a reference mel-spectrogram (usually the target
spectrogram, at training time). This distribution is learned by the model
through a bottleneck that is intended to capture prosody aspects of the target
spectrogram, preventing phonetic or speaker information to flow through. In
practice, specially when using this approach in multi-speaker settings, the
model tends to leak speaker information to the output. This represents a
problem at inference time, when a synthetic neutral reference is provided
(generally synthesized using a production system similar to Amazon Polly),
given that the synthesized speaker identity tends to resemble the reference
instead of the target voice. This problem is known as speaker leakage and it
was already reported in [10], where the authors emphasize the importance of
properly tuning the size of the reference bottleneck to amend it.

In [11], the authors propose using a multi-task version of Tacotron to enhance
the prosody of the syntheses. The approach described in the paper consists of
jointly learning the target mel-spectrogram as well as the probability
distribution of the phrase break patterns for each word. The work of [12]
proposes a novel training schema for Tacotron where deep style features are
extracted using the SER framework ([13, 14]). These features are later used
for minimizing the style differences between the real and synthesized samples.

In this work, we start from a multi-speaker Tacotron architecture with a
reference encoder (similar to [10]) and propose two modifications: (1) we
replace the pre-trained speaker embedding with a normalized speaker embedding
using normalizing flows [15], that allows sampling from the learned Gaussian
distribution and (2) we substitute the reference encoder with a residual
encoder, which learns a latent distribution conditioned on the input phonetic
information, allowing the model to capture and encode the residual attributes
not present in the linguistic input and the speaker embedding.

Inspired by the work of [16], we show how our proposed model achieves
significantly higher prosody variation by measuring the difference in variance
between the baseline and the proposed model across a set of features derived
from the fundamental frequency, the energy, the signal to noise ratio and the
speaking rate. These features capture most of the prosody variables (pitch,
speed, loudness and timbre) [16] and allow performing objective comparisons
between systems.

Our contribution can be summarized in the following points: (1) we propose a
new architecture that significantly increases the prosody variability of the
syntheses with no intelligibility and distinctiveness degradation, (2) the
proposed architecture removes the dependency on a production voice system as
it does not use a reference spectrogram, (3) when interpolating between
speaker embeddings the synthetic speakers from the proposed model are
significantly more distinct than the ones of the baseline.

##  2 Methods

A Tacotron2-based sequence to sequence model with location-based attention
[17] is used as a baseline throughout this study. This baseline architecture
is represented in figure 1-top. A couple of encoder branches are included over
the initial Tacotron formulation to allow the model to synthesize multiple
speakers: the speaker branch and the reference branch. The first one takes as
input a pre-trained speaker embedding vector representing the speaker
characteristics. This vector is obtained from the output of a speaker
verification model trained to minimize a triplet loss. This speaker
verification model is pre-trained using pairs of utterances, similar to [18].
For every speaker, the vectors corresponding to all their utterances are pre-
computed and then averaged to form the speaker vectors. As a result, we have a
fixed size vector for every speaker in the data set. Additionally, we include
a reference branch that is conditioned on the target spectrogram (at training
time) and is intended to learn a latent distribution summarizing the prosody
information, similar to the proposal of [10].

Figure 1: Top: the baseline architecture. Bottom: the proposed architecture.
The blue blocks represent pretrained parts of the network and the blocks with
round edges are tensor operations.

The proposed architecture builds upon the baseline model. We remove the
dependency of the reference branch on the target spectrogram, because it tends
to cause speaker/phonetic leakage (read the section 4 of the samples included
in [10]). It also slows down the inference process as it requires a production
TTS system to provide the references. Instead, we introduce a learnable
variational latent space conditioned on the phonetic information and a set of
learnable free parameters that we name residual branch. The purpose of this
branch is to encode the prosody information not present in the input
linguistic features or in the speaker embedding vector in a new latent space
(i.e. the residual prosody variance). Moreover, we normalize the speaker
embedding vectors using normalizing flows based on the work of [15], so that
the normalized vectors follow a Gaussian distribution.That change enables
sampling from the speaker embedding latent space instead of just using the
average embedding vector. This increases the coverage of the speaker embedding
space, allowing smoother interpolations between speaker embedding vectors. The
architecture proposed is depicted in figure 1-bottom

###  2.1 Speaker embedding normalization using normalizing flows

Normalizing flows [19] are powerful methods for producing tractable
distributions from which one can sample. Through a sequence of invertible
transformations they transform a complex input distribution to a tractable
probability distribution. The output distribution is usually chosen to be an
isotropic unit Gaussian to allow for smooth interpolations and efficient
sampling.

Glow [15], a flow-based generative model, has shown significant improvements
in computer vision generative modeling by adding invertible 1x1 convolutions
to the sequence of transformations applied to the input. To train this model
on one-dimensional speaker embedding vectors, we have switched the
2-dimensional convolutions with 1-dimensional ones.

We train a normalizing flow over the pre-trained speaker embeddings described
previously. The trained model attains a Gaussian distributed latent space of
speaker embeddings from which one can easily sample to create embeddings
representing new, unseen speakers. The proposed architecture samples from the
Gaussian distribution defined by the normalized embeddings, both at training
time and at inference time.

###  2.2 Residual branch

The new residual branch is designed to learn the prosody-induced variance that
cannot be explained by the linguistic features and the speaker embedding
vector alone. It takes as input the linguistic features and a set of learnable
free parameters111Notice that the dense layers that produce the residual
latent distribution parameters do not have bias terms, so that all the bias is
learned in the free parameters vector., conditioning on the target sentence
and on a global prior. The motivation behind this design is to help the
network learn different ways of uttering a given input sentence (linguistic
features conditioning), and the potential global biases existing in our
training dataset (free parameters).

The output of the phonetic encoder is piped into a bi-directional recurrent
neural network [20, 21] in order to get a representation independent of time.
The last output of the recurrent neural network of both, forward and reverse
passes (represented in equations 1 and 2 as $\mathbf{o_{f}}$ and
$\mathbf{o_{r}}$, respectively), are concatenated to form a fixed-size
phonetic representation $[\mathbf{o_{f}},\mathbf{o_{r}}]$. That vector is
concatenated with a vector of free parameters $\mathbf{v_{f}}$ and the result
is passed through two stacked dense layers to form the parameters of the
residual latent distribution (as shown in equation 3, where $f$ represents the
ReLU activation function). The $h_{residual}$ vector is split in two vectors
$h_{residual}^{\mu}$ and $h_{residual}^{\sigma}$ from where a latent vector
$z$ is sampled using the re-parametrization trick [22]. We add a Kullback-
Leibler loss to assure that the distribution of the latent representation
approximates a Gaussian distribution.

$\mathbf{o_{f}}=\text{RNN}_{f}(\mathbf{h_{ph}})$ (1)

$\mathbf{o_{r}}=\text{RNN}_{r}(\mathbf{h_{ph}})$ (2)

$\mathbf{h_{residual}}=\mathbf{W_{2}}(f(\mathbf{W_{1}}\cdot([\mathbf{o_{f}},\mathbf{o_{r}},\mathbf{v_{f}}])+\mathbf{b1}))+\mathbf{b_{2}}$
(3)

##  3 Experiments

We have used a combination of two internal datasets to train the models, one
of them containing 2860 non-professional speakers, with 200 utterances per
speaker on average, and the other containing 10 professional speakers, with
13,000 studio-recorded utterances per speaker – totaling to 2870 speakers and
more than 700,000 utterances. We process audio utterances at 16kHz and extract
80 dimensional mel-spectrograms. We define a frame as a $50ms$ sequence, with
an overlap of $12.5ms$. We used the Universal Neural vocoder to synthesize the
wav samples [23]. The linguistic features have been extracted using an
internal linguistic front-end which takes the text as input and extracts the
phonemes, that are used as input for the model. The speaker embedding, the
free parameters and the residual distribution parameters $\mu$ and $\sigma$
vectors were defined to have a length of 192 elements.

We have trained both models on one million batches, with a batch size of 24.
Then we synthesized a set of 50 unseen sentences for each of the 2870
speakers. The syntheses have been evaluated in 3 ways: intelligibility,
distinctiveness and prosody variability.

Additionally, we sampled 50 speakers randomly from the pool of 2870 speakers
and generated interpolations between all the possible pairs, in order to study
how the proposed model and the baseline behave in this setting.

###  3.1 Intelligibility

We used the AWS transcribe system to transcribe each of the syntheses, and
then measure the Word Error Rate (WER) between the target sentence and the
transcription [24]. These metrics are aggregated as the average WER per
speaker. We achieved a median Word Error Rate of $8.5\%$ in the baseline, and
$8.3\%$ in the proposed architecture. Although significative, this difference
is very small, thus we can conclude that both models are roughly equivalent in
terms of intelligibility.

Figure 2: Speaker intelligibility of the baseline and residual models,
measured as Word Error Rate of the target and the transcribed synthesis,
showing a slightly smaller error for the proposed architecture (Wilcoxon’s
$p=0.0305$).

###  3.2 Distinctiveness

To measure how distinct the synthesized speakers are, we use our internal
speaker verification system. Each utterance is compared with 4 samples
randomly drawn from the full pool of samples. We measure the False Acceptance
Rate (FAR), which quantifies the percentage of speaker pairs incorrectly
identified as the same speakers by the speaker verification model. Figure 3
shows how the FAR metric varies for both, the proposed model and the baseline
as we vary the classification threshold. Table 1 shows the FAR score for four
arbitrarily picked thresholds. Lower values of FAR mean higher
distinctiveness.

Figure 3: False Acceptance Rate distribution at different speaker verification model thresholds. Table 1: False Acceptance Rate metric for different thresholds Threshold | FAR-Baseline model | FAR-Proposed model
---|---|---
85% | 31.49% | 30.27%
90% | 29.50% | 28.21%
95% | 26.23% | 24.77%
99% | 19.05% | 16.93%

###  3.3 Prosody variability

To measure how much prosody and quality variance the new model introduces with
respect to the baseline, we have defined a set of objective metrics inspired
on the work of [16]. Those metrics are listed in table 2.

Figure 4: Baseline model vs proposed model when sampling from the residual latent distribution (blue) or from the speaker embedding distribution (orange), keeping the other constant. The results are represented as proportion of speaker-sentence pairs for which the proposed model shows significantly higher variance than the baseline, for each of the previously defined prosody features. Table 2: Metrics used to quantify the prosody variation and quality across different models. Metric |  Description
---|---
f0 mean |  Mean of the f0 (calculated using SPTK [25] considering only the vocal sounds (non-zero components), with frame skip of 12.5ms.
f0 range |  Difference between the 5th and 95th percentiles of f0.
speaking rate |  Calculated as the number of phonemes in the sentence divided by total duration of the synthesis.
f0 slope |  Calculated as the slope of a linear fit using least squares in the f0 plot.
snr |  Signal to noise ratio calculated using SOX [26] as the difference in RMS dBs between the loudest and the quietest windows, using windows size of 50ms.
power mean |  $20\log 10(\hat{x})$, where $\hat{x}$ is the average absolute amplitude with frame skip of 12.5ms. Only considering the vocal sounds.
power range |  Difference between the 5th and 95th percentiles of the power along time.
power slope |  Calculated as the slope of a linear fit using least squares in the power plot.

For this experiment, we have picked 50 speakers and 50 sentences. Then, 30
samples were synthesized for every speaker-sentence pair varying the random
seed, so different latent vectors are drawn from the latent distributions in
each repetition. We repeat this procedure four times: (1) with the baseline
model, (2) with the proposed model multiplying the $\sigma$ parameter of the
speaker embedding distribution by zero (so that the latent vector becomes the
mean of the speaker embedding latent distribution), (3) with the proposed
model multiplying the $\sigma$ parameter of the residual distribution by zero,
and (4) allowing sampling from both branches in the proposed model. We then
compare the syntheses of (1) vs (2), (1) vs (3) and (1) vs (4) .

For the comparison, as we are interested to measure differences of variance
for the variables defined in the table, we first use bootstrap to approximate
the variance distribution for every variable in every speaker-sentence pair,
and then we run a one-way Wilcoxon Signed-Rank sum test between the groups. We
use a significance level of $\alpha=5\%$ to which we apply the Bonferroni
correction ($\alpha=0.05/2500=0.002\%$). The results of the tests are
summarized in the figure 4.

We have informally observed by listening to the syntheses that when we sample
from the residual distribution (keeping the speaker embedding distribution
constant), the variations in the syntheses are related with prosody aspects
like the speaking rate, syllable duration or intonation, keeping, the speaker
identity is untouched. When sampling from the speaker embedding distribution,
we notice small variations on the speaker identity.

###  3.4 Speakers interpolation

We have tested the proposed model and the baseline using interpolated speakers
to study how the speaker embedding normalization affects the intelligibility
and distinctiveness of the new speakers. For that, 50 speakers have been
chosen. We have interpolated all the pairs of speakers using polar
interpolation, generating 1225 new voice profiles. We have synthesized 50
sentences using those new voices and then evaluated distinctiveness and
intelligibility (see results in figure 5).

Figure 5: Top: False Acceptance Rate distribution at different speaker
verification thresholds, for all the interpolated speakers. Bottom: speaker
intelligibility of the baseline and residual models, measured as Word Error
Rate of the target and the transcribed synthesis, showing a slightly smaller
error for the baseline architecture (Wilcoxon’s $p=5.2263\cdot 10^{-6}$)

The proposed model achieves 13.11% lower FAR than the baseline. We hypothesize
that the rationale behind this is that the normalized speaker embedding space
has a denser latent distribution than the one in the baseline. Given that the
proposed model is trained on samples drawn from the normalized speaker
embedding distribution (as opposed of using average vectors as in the baseline
model), the latent space has higher density and hence the generalization to
new speakers is better.

From the intelligibility perspective, the baseline shows an average WER of
$6.42\%$ while the proposed model achieved a $7.35\%$. We attribute that
difference to the fact that the interpolated speakers in the baseline are less
distinctive and resemble much more to one of the two actual speakers, hence
its intelligibility is naturally higher at the cost of a worse
distinctiveness. Although significative, the difference in practice is
negligible.

##  4 Conclusions

We introduced a new TTS architecture that allows increasing the prosody
variance of a multi-speaker TTS system by learning the residual prosody into a
new latent distribution. We showed that by sampling from the residual and
normalized speaker latent distribution the model produces syntheses with
significantly different prosodies as measured by a set of quantitative
metrics. The inclusion of the residual distribution also enables removing the
reference spectrogram dependency. This does not only solve potential speaker
leakage issues but also allows for a faster and potentially cheaper inference,
given that the model does not depend on a voice production system to work.

Finally, the normalization of the speaker embedding latent space allows for
better speaker interpolation when compared with the baseline model, thus
producing more diverse and unseen synthetic speakers.

## References

  * [1] B. Sisman, J. Yamagishi, S. King, and H. Li, “An overview of voice conversion and its challenges: From statistical modeling to deep learning,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 29, pp. 132–157, 2021.
  * [2] A. van den Oord, S. Dieleman, H. Zen, K. Simonyan, O. Vinyals, A. Graves, N. Kalchbrenner, A. Senior, and K. Kavukcuoglu, “Wavenet: A generative model for raw audio,” in _Arxiv_ , 2016. [Online]. Available: https://arxiv.org/abs/1609.03499
  * [3] D. Bahdanau, K. Cho, and Y. Bengio, “Neural machine translation by jointly learning to align and translate,” in _3rd International Conference on Learning Representations_ , ser. ICLR’15, July 2015.
  * [4] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, u. Kaiser, and I. Polosukhin, “Attention is all you need,” in _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , ser. NIPS’17, December, 2017, p. 6000–6010.
  * [5] S. Chaudhari, G. Polatkan, R. Ramanath, and V. Mithal, “An attentive survey of attention models,” _CoRR_ , vol. abs/1904.02874, 2019. [Online]. Available: http://arxiv.org/abs/1904.02874
  * [6] Y. Wang, R. Skerry-Ryan, D. Stanton, Y. Wu, R. J. Weiss, N. Jaitly, Z. Yang, Y. Xiao, Z. Chen, S. Bengio, Q. Le, Y. Agiomyrgiannakis, R. Clark, and R. A. Saurous, “Tacotron: Towards end-to-end speech synthesis,” 2017. [Online]. Available: https://arxiv.org/abs/1703.10135
  * [7] J. Shen, R. Pang, R. J. Weiss, M. Schuster, N. Jaitly, Z. Yang, Z. Chen, Y. Zhang, Y. Wang, R. Skerry-Ryan, R. A. Saurous, Y. Agiomyrgiannakis, and Y. Wu, “Natural tts synthesis by conditioning wavenet on mel spectrogram predictions,” 2018. [Online]. Available: https://arxiv.org/abs/1712.05884
  * [8] R. Liu, B. Sisman, J. Li, F. Bao, G. Gao, and H. Li, “Teacher-student training for robust tacotron-based TTS,” vol. abs/1911.02839, 2019. [Online]. Available: http://arxiv.org/abs/1911.02839
  * [9] P. Taylor, _Text-to-Speech Synthesis_ , 1st ed. USA: Cambridge University Press, 2009.
  * [10] R. Skerry-Ryan, E. Battenberg, Y. Xiao, Y. Wang, D. Stanton, J. Shor, R. Weiss, R. Clark, and R. A. Saurous, “Towards end-to-end prosody transfer for expressive speech synthesis with tacotron,” in _Proceedings of the 35th International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, J. Dy and A. Krause, Eds., vol. 80. Stockholmsmässan, Stockholm Sweden: PMLR, 10–15 Jul 2018, pp. 4693–4702. [Online]. Available: http://proceedings.mlr.press/v80/skerry-ryan18a.html
  * [11] R. Liu, B. Sisman, F. Bao, G. Gao, and H. Li, “Modeling prosodic phrasing with multi-task learning in tacotron-based tts,” _IEEE Signal Processing Letters_ , vol. 27, p. 1470–1474, 2020. [Online]. Available: http://dx.doi.org/10.1109/LSP.2020.3016564
  * [12] R. Liu, B. Sisman, G. Gao, and H. Li, “Expressive tts training with frame and style reconstruction loss,” 2020.
  * [13] S. Zhang, S. Zhang, T. Huang, and W. Gao, “Speech emotion recognition using deep convolutional neural network and discriminant temporal pyramid matching,” _IEEE Transactions on Multimedia_ , vol. 20, no. 6, pp. 1576–1590, 2018.
  * [14] R. Lotfian and C. Busso, “Curriculum learning for speech emotion recognition from crowdsourced labels,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 27, no. 4, pp. 815–826, 2019.
  * [15] D. P. Kingma and P. Dhariwal, “Glow: Generative flow with invertible 1x1 convolutions,” in _Advances in Neural Information Processing Systems_ , S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, Eds., vol. 31. Curran Associates, Inc., 2018, pp. 10 215–10 224.
  * [16] T. Raitio, R. Rasipuram, and D. Castellani, “Controllable neural text-to-speech synthesis using intuitive prosodic features,” 2020. [Online]. Available: https://arxiv.org/pdf/2009.06775.pdf
  * [17] M.-T. Luong, H. Pham, and C. D. Manning, “Effective approaches to attention-based neural machine translation,” 2015.
  * [18] Z. Ren, Z. Chen, and S. Xu, “Triplet based embedding distance and similarity learning for text-independent speaker verification,” in _2019 Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA ASC)_ , 2019, pp. 558–562.
  * [19] D. J. Rezende and S. Mohamed, “Variational inference with normalizing flows,” 2016\.
  * [20] M. Schuster and K. K. Paliwal, “Bidirectional recurrent neural networks,” _IEEE Transactions on Signal Processing_ , vol. 45, no. 11, pp. 2673–2681, 1997.
  * [21] A. Graves, S. Fernández, and J. Schmidhuber, “Bidirectional lstm networks for improved phoneme classification and recognition,” in _Artificial Neural Networks: Formal Models and Their Applications – ICANN 2005_. Springer Berlin Heidelberg, 2005, pp. 799–804.
  * [22] D. P. Kingma and M. Welling, _An Introduction to Variational Autoencoders_ , 2019.
  * [23] J. Lorenzo-Trueba, T. Drugman, J. Latorre, T. Merritt, B. Putrycz, R. Barra-Chicote, A. Moinet, and V. Aggarwal, “Towards achieving robust universal neural vocoding,” 2019.
  * [24] U. Kamath, J. Liu, and J. Whitaker, _Deep Learning for NLP and Speech Recognition_. Springer, 2019. [Online]. Available: https://doi.org/10.1007/978-3-030-14596-5
  * [25] S. Imai, T. Kobayashi, and K. Tokuda, “Speech signal processing toolkit (SPTK) - version 3.11,” https://sp-tk.sourceforge.net, 2017.
  * [26] “Sound eXchange (SOX),” https://sox.sourceforge.net, 2015.