# Graph and Network Theory for the analysis of Criminal Networks

Lucia Cavallaro
University of Derby
Derby  UK
<EMAIL_ADDRESS>
&Ovidiu Bagdasar
University of Derby
Derby  UK
<EMAIL_ADDRESS>
&Pasquale De Meo
University of Messina
Messina  Italy
<EMAIL_ADDRESS>
&Giacomo Fiumara
University of Messina
Messina  Italy
<EMAIL_ADDRESS>
&Antonio Liotta
Free University of Bozen-Bolzano
Bolzano  Italy
<EMAIL_ADDRESS>

###### Abstract

Social Network Analysis is the use of Network and Graph Theory to study social
phenomena, which was found to be highly relevant in areas like Criminology.
This chapter provides an overview of key methods and tools that may be used
for the analysis of criminal networks, which are presented in a real-world
case study. Starting from available juridical acts, we have extracted data on
the interactions among suspects within two Sicilian Mafia clans, obtaining two
weighted undirected graphs. Then, we have investigated the roles of these
weights on the criminal network’s properties, focusing on two key features:
weight distribution and shortest path length. We also present an experiment
that aims to construct an artificial network that mirrors criminal behaviours.
To this end, we have conducted a comparative degree distribution analysis
between the real criminal networks, using some of the most popular artificial
network models: Watts-Strogatz, Erdős-Rényi, and Barabási-Albert, with some
topology variations. This chapter will be a valuable tool for researchers who
wish to employ social network analysis within their own area of interest.

## Introduction

Graph Theory is a well established field in mathematics. However, only
recently many of its theoretical results started to be used within Social
Network Analysis (SNA), an area with significant implications for real world
scenarios. For example, one can simulate the behaviour of social networks
using strategies like link predictions LU20111150 ; Hasan2011 , temporal
networks, or spreading of influences Tassiulas2013 ; Cavallaro2020ArXiv .
Other practical applications include to deal with large Artificial Neural
Networks Mocanu2018network ; Cavallaro2019network ; Cavallaro2020network  or
targeted advertisements to people based on their friends’ interests Bellur2007
or, on the other side, containing the spread of fake news zhou2018fake .

Network Science tools may also be used in the investigation of criminal
networks. Sometimes the complex social interactions within a clan-based
society may help the feature selection process for building machine learning
models Oluwabunmi2019IEEE . Other times, it is Network Science itself that
helps conducting better performing investigation from law enforcement
agencies. To this end, criminal networks can be encoded as graphs, and various
types of analysis and simulations can be carried out for modelling criminal
behaviours.

This chapter is to intended as a short tutorial on how Network Science
strategies may be used to conduct an in-depth analysis on real criminal
networks. Here, the Sicilian Mafia scenario has been considered. Sect. Case
study: The Sicilian Mafia relates to our previous analysis on this topic. In
particular, Sect. 0.0.8 includes a comparative Degree Distribution analysis
between real criminal networks and artificial ones. Indeed, when it is
possible to find out a synthetic network reflecting the behaviour of real-
world criminal networks, law enforcement agencies (LEAs) and network
scientists can recreate those networks and simulate how interconnections among
criminal will evolve.

This chapter is structured as follows. Sect. Complex Networks presents the key
theoretical tools (required for understanding the experiments conducted in
Sect. Case study: The Sicilian Mafia), and it is divided in two parts:  (i)
tools, where the basic definitions on network science are provided; and  (ii)
popular artificial networks description (as the topologies used in Sect.
0.0.8).  Next, in Sect. Social Network Analysis in Criminal Networks a brief
review on the use of  (i) Social Network Analysis, and its implication in
(ii) Criminal Networks is defined.  Sect. Case study: The Sicilian Mafia is a
case study summarizing our work on two real criminal networks related to
Sicilian Mafia Ficara2020 ; cavallaro2020disrupting . This section includes
four parts:  (i) datasets description, based on the data extracted from
juridical acts;  (ii) weights distribution analysis, which represents an
important preliminary study to understand how the interactions among suspected
are structured in terms of interaction frequency;  (iii) shortest path
analysis, which allows to identify trusted affiliates inside the clan who can
spread confidential and illegal messages;  (iv) comparative degree
distribution analysis between real and synthetic networks, to artificially
recreate the criminal networks used here, with the purpose of conducting
further investigations through them.  Finally, the conclusions follow in Sect.
0.0.8.

## Complex Networks

In this section we introduce the main Network Science concepts underpinning
Social Network Analysis, which are later exemplified in the criminal network
case study (Sect. Case study: The Sicilian Mafia). In particular, in Sect.
0.0.1 the main definitions required for understanding the mechanics of SNA are
provided. All theoretical concepts are derived from barabasi2016network ,
which we refer to the reader for further technical details.

###  0.0.1 Tools

In this section, we start with some basic definitions of Graph Theory. Graph

###### Definition 1

A graph denoted by $G=(N,E)$, consists of a set of nodes $N$ and a set of
edges $E\subseteq N\times N$ (also called links $L$). It is a convenient way
of representing relationships between pairs of objects.

As an example, Facebook® may be viewed as a graph, where the nodes represent
users and edges represent the friendship relationship among them. Is it also
possible to define a subgraph as follows:

###### Definition 2

A subgraph $H$ of the graph $G$ is a graph whose nodes and edges are subsets
of the nodes and edges of $G$.

Furthermore, graphs may be either weighted, or unweighted:

###### Definition 3

A weighted graph $G=(N,E,W)$ is a triplet consisting of a finite set of nodes
$N$, a set of edges $E$, and a set of weights $W:E\rightarrow\mathbb{R}$
defined on each edge. If all edges weights are equal to one, then the graph is
called unweighted.

Degree

###### Definition 4

The degree of a node $n_{i}$, denoted $deg(i)$ or $k_{i}$, is the number of
incident edges to $n_{i}$. The sum of the degrees of all nodes is equal to the
double of the number of edges $E$:

$\sum_{n\>\in\>N}k_{n}=2E.$ (1)

###### Definition 5

In weighted networks, the weighted degree (also known as strength
antoniou2008statistical ; barrat2005 ) is the sum of the edges weights $w$
incident on $n_{i}$:

$k_{i}=\sum_{(i,j)\in E}w_{ij},$ (2)

where the summation spans over all edges $(i,j)$ in the network, linked to
node $n_{i}$.

###### Definition 6

For undirected networks, the average degree is defined as

$\langle k\rangle=\frac{1}{N}\sum_{i=1}^{N}k_{i}=\frac{2L}{N},$ (3)

where $N$ is the total number of nodes, $k_{i}$ is the degree of a generic
node $i$, and $L$ represents the total number of links, or edges $E$, within
the network.

#### Small-World

The Small World phenomenon smallworld1967 ; smallworld1969  is based on the
concept of the six degrees of separations, according to which two random
people in the world may be connected each other via a few acquaintances (i.e.,
it is estimated that there are six people in the middle between the source and
the destination). In Network Science, it translates into a “short” distance
between two randomly chosen nodes within a network, that is

$\langle d\rangle\approx\frac{\ln N}{\ln\,\langle k\rangle},$ (4)

where $N$ is the total number of nodes in the graph, $\langle k\rangle$ is the
network average degree, and $\langle d\rangle$ the average distance within the
network. The denominator implies that the denser the network, the smaller the
distance between the nodes is. In conclusion, the average path length or the
diameter depends logarithmically on the system size.

##### Degree Distribution

The degree distribution $p_{k}$ provides the probability that a randomly
selected node in the network has degree $k$. Since $p_{k}$ is a probability,
it must be normalized; i.e.,

$\sum_{i=1}^{\infty}p_{k}=1.$ (5)

For a network made of $N$ nodes, the degree distribution is the normalized
histogram given by:

$p_{k}=\frac{N_{k}}{N},$ (6)

where $N_{k}$ is the number of nodes having degree $k$.

The degree distribution has assumed a central role in network theory following
the discovery of scale-free networks (See “Scale-Free Property” paragraph);
moreover, $p_{k}$ determines many network phenomena, from network robustness
to the spread of viruses.

##### Weight Distribution

The degree distribution can be extended to weighted networks considering the
weighted degree (strength) distribution $P(s)$, defined as the probability
that a node may have weighted degree (strength) equal to $s$. Based on
barrat2005 , this is

$P(s)\sim s^{-\gamma},$ (7)

where $\gamma$ is a constant typical of the network.

Clustering Coefficient Clustering is used to quantify the relationship among
nodes’ neighbours. Indeed, the degree only considers the number of direct
links between nodes. The clustering coefficient $C_{i}$ measures the edge
density in the immediate neighbourhood of a node. $C_{i}\in[0,1]$ represents
the clustering coefficient of a generic node $n_{i}$:

$\begin{cases}$if $C_{i}=0,&$there are no edges among the node's
neighbours$\\\ $if $C_{i}=1,&$ each node's neighbour is connected with the
others$\end{cases}$

The local clustering coefficient is computed as follow:

$C_{i}=\frac{2L_{i}}{k_{i}(k_{i}-1)},$ (8)

where $k_{i}$ is the degree of the generic node $n_{i}$, and $L_{i}$
represents the number of links (i.e., edges) between the $k_{i}$ neighbours of
$n_{i}$.

##### Average Clustering Coefficient

The average $\langle C\rangle$ of $C_{i}\in i=1,\dots,N$ in the whole network
is given by

$\langle C\rangle=\frac{1}{N}\sum_{i=1}^{N}C_{i}.$ (9)

Adjacency Matrix A common way to represent relationships among nodes is the
adjacency matrix $A$.

###### Definition 7

The Adjacency Matrix $A[i,j]$ holds node degree (weighted or unweighted) to
the edge $(n_{i},n_{j})$ if it exists, where $n_{i}$ is the node with index
$i$ and $n_{j}$ is the node with index $j$. If there is no such edge, then
$A[i,j]=None$.

For undirected graphs $A$ is symmetric (i.e., $A[i,j]=A[j,i]$ $\forall
n_{i},n_{j}\in N$).

Path

###### Definition 8

A path is a sequence of alternating nodes and edges that flow from a starting
node to an ending one such that each edge is incident to its predecessor and
successor node. A path is called simple if each node in the path is distinct.

More formally, a path can be defined as a sequence of nodes

$P=(n_{1},n_{2},\dots,n_{m})\in N\times N\times\cdots\times N,$

such that $n_{i}$ is adjacent to $n_{i+1}$ for $1\leq i\leq m-1$. Such a path
$P$ is called a path of length $m-1$ from $n_{1}$ to $n_{m}$.

Measures based on paths strategies are the shortest path length analysis.

Distance

###### Definition 9

The distance from a node $n_{i}$ to a node $n_{j}$ in $G$, denoted
$d(n_{i},n_{j})$ is the length of a shortest path from $n_{i}$ to $n_{j}$ (if
such a path exists).

$d_{ij}=min\left(\Gamma(i,j)\right),$

where $\Gamma(i,j)$ is the set of paths connecting $i$ and $j$.

Connectedness

###### Definition 10

A graph $G$ is connected if, for any two nodes, there is a path between them.

###### Definition 11

If $G$ is not connected, its maximal connected subgraphs are called the
connected components of $G$.

###### Definition 12

If a network consists of two components, a properly placed link can connect
them, making the network connected. Such a link is called bridge.

Scale-Free Property The majority of real networks, such as the World Wide Web,
are called scale-free networks and follow the definition:

###### Definition 13

A scale-free network is a network whose degree distribution follows a power
law.

The power-law distribution has the following form

$p_{k}\sim k^{-\gamma},$ (10)

where the exponent $\gamma$ is its degree exponent.

Some artificial network models such as the Barabási-Albert (BA) Model
successfully exhibit this feature.

###  0.0.2 Artificial Networks

The need for scientists to create Artificial, or Synthetic Networks has been
born from the aim to reproduce real network properties in a controlled
environment. For this reason, several typologies of Artificial Networks have
been formulated.

Three models in particular have found special popularity within the scientific
community: The Erdős-Rényi (ER, also known as Random Network) Model, the
Watts-Strogats Model (WS; i.e., a Random Network variation), and the Barabási-
Albert (BA) Model. This last one tries to capture two important properties of
real network: the growth and the preferential attachment. Further details on
those models are provided in the following paragraphs. Random Network Model A
random network consists of $N$ nodes where each node pair
$(n_{i},n_{j}),\forall i,j\in N$ is connected with probability $p$. To
construct a random network one needs to

  1. 1.

Start with $N$ isolated nodes,

  2. 2.

Select a node pair $(n_{i},n_{j})$ and generate a random number
$rand\in[0,1]$:

$\begin{cases}$if $rand>p,&$connect the selected node pair with a link$\\\
$otherwise$,&$leave them disconnected$\end{cases}$

  3. 3.

Repeat the previous step for all pairs of distinct nodes $(n_{i},n_{j})\in
N\times N$.

The network obtained after this procedure is called a random graph or a random
network. There are two definitions of a random network: the definition
provided in the Erdős-Rényi Model, and the one of the Gilbert Model.

####  Erdős-Rényi Model

Random networks are also called Erdős-Rényi Networks from the names of the
mathematicians Paul Erdős (1913-1996) and Alfréd Rényi (1921-1970), who
studied the properties of these networks. Their model follow the structure

$G(N,L),$ (11)

where $N$ labeled nodes are connected with $L$ randomly placed links (i.e.,
edges). Paul Erdős and Alfréd Rényi used this definition in their paper
random1959 .

####  Gilbert Model

It is a variation of the Erdős-Rényi Model. It has been defined by Edgar
Nelson Gilbert (1923-2013) and follows the structure

$G(N,p),$ (12)

where each pair of $N$ labeled nodes is connected with probability $p$.

There are two main limits in Random Network Model that had to be overcome over
the years by the academic community:

  1. 1.

The local clustering coefficient in ER model is given by barabasi2016network

$C_{i}=\frac{\langle k\rangle}{N}.$

This behaviour of $C_{i}$ is contradicted by the local clustering coefficient
of real networks.

  2. 2.

The Poisson distribution that describes the degree distribution of ER networks
does not allow large differences between the worst- and best-connected nodes
in the network. This implies that hubs, frequently observed in real networks,
cannot be found in ER networks. BA model, relying on preferential attachment
and growth, successfully reproduces this fundamental feature.

In the following paragraphs those models are described.

Watts-Strogats Model Two main considerations motivated Duncan J. Watts (1971)
and Steven Strogatz (1959) to propose this model:  (i) in real networks the
average distance between two nodes depends logarithmically on $N$ (See “Small-
World” average distance $\langle d\rangle$);  (ii) the average clustering
coefficient $\langle C\rangle$ of real networks is much higher than expected
for a random network of similar $N$ and $L$ (i.e., $E$).

To construct a random network according to Watts-Strogats Model watts1998 :

  1. 1.

Start from a ring of $N$ nodes, whereas each node is connected to its
immediate previous and next neighbours; hence, each node has $\langle
C\rangle=3/4$, initially.

  2. 2.

With probability $p\in[0,1]$ each link is rewired to a randomly chosen node

$\begin{cases}$if $p\simeq 0,&$regular lattice$\\\ $if $0<p<1,&$Small-World
property$\\\ $if $p=1,&$Random Network Model (all links rewired).$\end{cases}$

The Watts-Strogatz model interpolates between a regular lattice, which has
high clustering (but lacks the Small-World phenomenon), and a random network,
which has low clustering (but displays the Small-World property). Moreover,
high nodes degrees are absent from Watts-Strogatz model.

Barabási-Albert Model This model was theorized by Albert-László Barabási
(1967) and Réka Albert (1972) barabasi1999 . It simulates a Scale-Free Network
rather than a Random one, by introducing two new concepts to the model:
network growth and the preferential attachment. The first concept assumes that
real networks continuously increase over time, so new nodes must be
considered. The second point argues that random connections defined by a fixed
probability do not reflect the behaviour of real networks; in fact, in real
scenarios, nodes tend to link to the more connected nodes.

###### Definition 14

The Preferential Attachment is the probability $\Pi(k)$ that a link of the new
node $n_{j}$ connects to node $n_{i}$ depends on the degree $k_{i}$ through
the formula

$\Pi(k_{i})=\frac{k_{i}}{\sum_{j}k_{j}}.$ (13)

To construct an artificial network with the Barabási-Albert model, the steps
are:

  1. 1.

Start with a set of $N_{0}$ nodes, the links between which are chosen
arbitrarily, as long as each node has at least one link.

  2. 2.

Growth – At each timestep a new node $n_{j}$ with $l$ links (with $l\leq
l_{0}$) that connects the new node to nodes already in the network is added.

  3. 3.

The connections between the new node with the older nodes are defined by the
Preferential Attachment probability.

## Social Network Analysis in Criminal Networks

In this section we provide an overview of state-of-the-art of Social Network
Analysis applied to Criminal Networks. We also consider the most relevant
studies concerning specifically the Sicilian Mafia criminal topologies.

###  0.0.3 Criminal Networks Analysis

Through SNA, LEAs are able to analyze criminal networks and investigate the
relations among criminals. For this reason, nowadays there is a growing
interest in the application of Graph and Network Science onto criminal
networks. For instance, SNA has been used in Chen2004  to build crime
prevention systems. However, due to the lack of data availability on those
kind of networks, there are difficulties in finding relevant quantitative
studies. Such examples are those conducted by Szymanski Szymanski2018  and
Berlusconi BerlusconiPLOS2016 , on the problem of community detection and link
prediction.

###  0.0.4 Sicilian Mafia networks

Sicilian Mafia has a particular structure that differs from common criminal
networks (such as the terrorist nets), whereby it is a common practice for
criminals to come together to achieve a common goal and then fall apart. By
contrast, in Sicilian Mafia this behaviour does not occur. Indeed, the
affiliates are bound by blind loyalty and they still pursue further goals even
after achieving a previous one. Moreover, Families last for several
generations. They also tend to diversify their objectives: from controlling
entire economic sectors (e.g., by giving “protection” to small traders and
taking control of larger factories), to influencing countries political life
(e.g., by interfering in the results of electoral competitions). Sicilian
Mafia originated in Sicily, and has now spread worldwide franchetti1925sicilia
; McGloin2005 ; Mastrobuoni2012 . The blind loyalty of affiliates makes it
even more difficult to obtain reliable information about those criminal
networks topologies: important information about such criminal network is
likely to be missing or hidden, due to the covert and stealthy nature of
criminal actions krebs2002 ; Xu2005 ; Calderoni2010 ; Campana2012 .

## Case study: The Sicilian Mafia

This section describes firstly the real criminal datasets we used for our
tests, followed by a brief summary on the strategies conducted jointly with
the results obtained so far as an example on Network Science strength, and how
it can be used to significantly help LEAs. In particular, the experiments
relate to:  (i) weight distribution,  (ii) shortest path length, and  (iii)
degree distribution.

###  0.0.5 Dataset Description

The case study example relates to two real-world datasets we built from
juridical acts111Datsets and source code are available at
https://github.com/lcucav/criminal-nets cavallaro2020disrupting :  (i) the
Meetings dataset represents the physical meetings among criminals obtained
through LEA evensdropping;  (ii) the Phone Calls dataset refers to phone calls
between individuals obtained through LEA interceptions.

This particular investigation was a prominent operation conducted during the
first decade of the 2000s and focused on two Mafia clans known as the
“Mistretta” family, and the “Batanesi”clan online-journ2019 ; Ficara2020 ;
cavallaro2020disrupting .

Both datasets led to undirected and weighted graphs; thus, edge weights $w$
are available and represent the number of times any given pair had a meeting
in the Meetings dataset, and the number of times two individuals called each
other in Phone Calls dataset. In SNA, those coefficients are also known as the
strength of the tie binding two individuals Ficara2020 . In Fig. 1, the graphs
obtained as well as the description of what each element represents (i.e.,
nodes, colours, edges weight, nodes and edges size, etc.) is shown. The main
characteristics of the datasets are summarized in Table 0.0.5.

(a) Meetings

(b) Phone Calls

Figure 1: Dataset Description. The colours represent different clans: darker
nodes are the “Mistretta” family; in grey the “Batanesi” clan is drawn; white
and light gray circled nodes and for two others Mafia families not directly
involved in the current investigation. All circled nodes represent the bosses.
Lastly, white nodes represent other subjects not classifiable in any of the
previous categories. Edges’ width depends on the number of meetings or phone
calls, while the nodes size relates with their degree. (Reproduced from Ficara
_et al._ 2020) Table 1: Characteristics of Meetings and Phone Calls networks.
(Reproduced from Ficara _et al._ 2020)

Parameter |  Meetings |  Phone Calls
---|---|---
No. Nodes |  101 |  100
No. Edges |  256 |  124
Max. Weight |  10 |  8
Max. Frequency |  200 |  100
Avg. Degree |  5.07 |  2.48
Max. Shortest Path |  7 |  14
Common nodes | 47

###  0.0.6 Weight Distribution Analysis

Fig. 2 shows the weight distribution of the Meetings and the Phone Calls
networks. As already mentioned, the weights represent the amount of meetings
and phone calls exchanged between pairs of individuals in the networks,
respectively.

It is noteworthy that in both these networks there are just a few high-weight
edges i.e.; nodes incident on those links exhibit an high number of
interaction within the network. In Ficara2020 , we motivated this behaviour as
a necessity from affiliates to focus their efforts in trying to reduce the
risk of being intercepted by external people (i.e., LEA, and other people
outside the clan). In the Meetings network, this trend is even more
accentuated; moreover, the maximum interactions weight (i.e., $w=10$) is
greater than its counterpart in the Phone Calls network (i.e., $w=8$).

Our explanation is that mobsters prefer to communicate by face-to-face
meetings, rather than calling each other, to reduce interception risks.
Furthermore, bosses often have to participate to public events to pursue their
power inside a clan, including: funerals of other affiliates, and other solemn
religious demonstrations (masses, processions, etc.). It is a well known
practice that, during those kinds of events, bosses pass messages to their
closest subordinate affiliates.

(c) Meetings

(d) Phone Calls

Figure 2: Weight Distribution in the Meetings dataset (a), and the Phone Calls
dataset (b). (Reproduced from Ficara _et al._ 2020)

###  0.0.7 Shortest Path Length Analysis

The shortest path length distribution in Fig. 3 is closely related to dynamic
properties such as velocity of messages spreading process within the network.
Generally speaking, the criminal organizations structure aims to to optimize
the interaction frequency among members, while reducing as much as possible
the interception risks. Thus, trusted members may be discovered by following
short interactions paths; indeed, those affiliates may also be acting as a
bridge (See Sect. 0.0.1“Connectedness”) to connect distant groups in the
network.

In Ficara2020  we noticed that both the weighted and the unweighted shortest
path length analyses show a higher interaction frequency among affiliates
throughout a “balanced” number of intermediates. This means that they do not
like to spread their encrypted messages with a too low (resp., high) number of
intermediates. This is to avoid, from one side, to overexpose their bosses to
police investigations. From the other side, the longer the sequence of
intermediates, the higher the chances to be intercepted by people outside the
Family.

Even through this analysis, as it was for the weights distribution, it emerged
that the clan tries to minimize the risk of interceptions, especially to avoid
exposing those mobsters who are hierarchically in a higher rank.

(a) Unweighted graphs.

(b) Weighted graphs.

Figure 3: Distribution of shortest path lengths in Meetings and Phone Calls
networks in unweighted (a) and weighted (b) graphs. (Reproduced from Ficara
_et al._ 2020)

###  0.0.8 Degree Distribution Analysis

The Degree Distribution Analysis has been conducted in order to discover an
appropriate artificial network that would virtually mirror the real-world
criminal graphs topology under scrutiny. We previously analysed the weight
distribution, but due to lack of libraries available for weighted graphs
analysis
222https://networkx.github.io/documentation/networkx-1.9/reference/generators.html,
we opted for a preliminary analysis on nodes degree distribution. To this end,
in Fig. 4 we compared our real criminal networks against five artificial
models:  (i) the Random Network by Gilbert (G-ER),  (ii) Watts-Strogatz (WS),
(iii) its variant accordingly with Newmann NEWMAN1999341  (N-WS), and  (iv)
two different configurations of the Barabási-Albert (BA) model in terms of
links added at each step; BA2 with $m=2$, and BA3 with $m=3$.  Table 0.0.8 and
Table 0.0.8 summarize the number of edges and average degree obtained with the
configurations above described in Phone Calls and Meetings graphs,
respectively.

Note that all the results herein shown represent the average results obtained
after 100 runs per each synthetic network.

We initially compared the real networks with the Erdős-Rényi (ER) topology,
but this model did not allow to customize the number of links. Then, we opted
for the Gilbert one, whereby both number of nodes $n$ and links $m$ are
defined a priori.

In the WS model, we set $n$, $k=\frac{2m}{n}$ (that represents the number of
nearest neighbors links per node), and the rewiring probability $p=0.5$, with
$p\in[0,1]$. As previously asserted in Sect. 0.0.2, if $p=1$, we turn into a
Random Network. The main difference between WS and N-WS models is that in WS,
number $p$ is the probability of rewiring each edge; whereas in N-WS, $p=0.25$
is the probability of adding a new edge for each edge. Indeed, if $p=1$, then
number of edges is doubled.

(a) Meetings

(b) Phone Calls

Figure 4: Degree Distribution in the Meetings dataset (a), and the Phone Calls
datasets (b). Circles give the actual datasets values. G-ER is the Random
Network proposed by Gilbert. WS is the Watts-Strogatz network. N-WS is the
Newmann variation of WS. BA2 and BA3 are the Barabási-Albert models with $m=2$
and $m=3$, respectively. Table 2: Characteristics of Artificial Models in the
Phone Calls network.

Model |  No. Edges |  Avg. Degree
---|---|---
G-ER |  124 |  2.48
WS |  100 |  2.00
N-WS |  123 |  2.46
BA2 |  196 |  3.92
BA3 |  291 |  5.82

Table 3: Characteristics of Artificial Models in the Meetings network.

Model |  No. Edges |  Avg. Degree
---|---|---
G-ER |  256 |  5.07
WS |  202 |  4.00
N-WS |  250 |  4.95
BA2 |  198 |  3.92
BA3 |  294 |  5.82

The actual graphs degree distributions act very differently from one another.
In particular, the fluctuations in Meetings are justified by the fact that
face-to-face encounters have been observed not only between couple of
suspects, but also among groups of more than two people at the same time. On
the other hand, phone calls have only been considered between individual
suspects. The analysis suggests that through degree distribution it was not
possible to identify an appropriate artificial network that best fits the
network characteristics of the two real-world datasets considered. This is
mainly due to the size of the networks. In fact, artificial networks seem to
work better with lager sizes; thus, are quite unstable in the first step of
their creation. For example, the emergence of hubs in BA models cannot be
highlighted because of the small size of the overall network obtained.

## Conclusions

This chapter aims to showcase the applicability of Graph Theory in
Criminology, during a time when the use of SNA by LEAs is growing
substantially. The case study herein reported as an example, explores
different approaches on criminal networks analysis by means of network science
tools.

In our study we have first created a graph from data extracted from juridical
acts; then, we started a twofold preliminary investigation: a weight
distribution analysis, and in parallel, a shortest path length analysis. These
have been conduced to identify the extent by which weighted graphs are useful
in those small networks.

Thus, we conducted a comparative degree distribution analysis between our
real-world networks and some models generated by popular artificial networks.
The aim was to identify the appropriate synthetic network which could simulate
criminal networks artificially, but in an effective manner. The strength of
this idea is that we may also be able to understand the patterns followed by
criminals to create their internal interconnections among affiliates. Our
study has found that the network size is a limitation. Indeed, there are
significant fluctuations and through degree distribution comparative analysis
it was not possible to find an appropriate artificial network that accurately
mirrors the two real criminal networks used in our tests.

To overcome this issue, in future studies we will investigate adjacency matrix
structures for both real and synthetic networks to get insights into network
topologies.

## References

  * (1) Linyuan, L., Zhou, T.: Link prediction in complex networks: A survey. Physica A: Statistical Mechanics and its Applications. 390(6) 1150–1170 (2011) doi: 10.1016/j.physa.2010.11.027
  * (2) Hasan, M.A., Zaki, M.J.: A Survey of Link Prediction in Social Networks. In: Aggarwal C. (eds) Social Network Data Analytics. Springer, Boston, MA (2011)
  * (3) Tassiulas, L., Katsaros, D., Basaras, P.: Detecting Influential Spreaders in Complex, Dynamic Networks. In: Computer, 46(4) 24–29 (2013) doi: 10.1109/MC.2013.75
  * (4) Cavallaro, L., Costantini, S., De Meo, P., Liotta, A., Stilo, G.: Network connectivity under a probabilistic node failure model. In: ArXiv e-print (Jun 2020). arXiv: 2006.13551 [cs.SI]
  * (5) Mocanu, D.C., Mocanu, E., Stone, P., Nguyen, P.H., Gibescu, M., Liotta, A.: Scalable training of artificial neural networks with adaptive sparse connectivity inspired by network science. Nature communications 9(1) 1–12 (2018)
  * (6) Cavallaro, L., Bagdasar, O., De Meo, P., Fiumara, G., Liotta, A.: Network Science Strategies for Accelerating the Training of Artificial Neural Networks. In: Numerical Computations: Theory and Algorithms NUMTA 2019. 169 (2019)
  * (7) Cavallaro, L., Bagdasar, O., De Meo, P., Fiumara, G., Liotta, A.: Artificial Neural Networks Training Acceleration Through Network Science Strategies. In: Sergeyev, Y.D., Kvasov, D.E. (eds) Numerical Computations: Theory and Algorithms. 330–336. Springer International Publishing, Cham (2020)
  * (8) Bellur, U., Kulkarni, R.: Improved Matchmaking Algorithm for Semantic Web Services Based on Bipartite Graph Matching. IEEE International Conference on Web Services (ICWS 2007), Salt Lake City, UT, 86—93 (2007)
  * (9) Zhou, X., Zafarani, R.: Fake News: A Survey of Research, Detection Methods, and Opportunities. In: ArXiv e-print (Dec 2018). arXiv: 1812.00315 [cs.CL]
  * (10) Oluwabunmi, O., Cosma, G., Liotta, A.: Clan-Based Cultural Algorithm for Feature Selection. In: 2019 International Conference on Data Mining Workshops (ICDMW). 465–472. IEEE (2019)
  * (11) Ficara, A., Cavallaro, L., De Meo, P., Fiumara, G., Catanese, S., Bagdasar, O., Liotta A.: Social Network Analysis of Sicilian Mafia Interconnections. In: Cherifi H., Gaito S., Mendes J., Moro E., Rocha L. (eds) Complex Networks and Their Applications VIII. COMPLEX NETWORKS 2019. Studies in Computational Intelligence. Springer, Cham 882 440–450. (2020) doi: 10.1007/978-3-030-36683-4_36
  * (12) Cavallaro, L., Ficara, A., De Meo, P., Fiumara, G., Catanese, S., Bagdasar, O., Song, W., Liotta, A.: Disrupting Resilient Criminal Networks through Data Analysis: The case of Sicilian Mafia. PLoS ONE (In Press) (2020) https://doi.org/10.1371/journal.pone.0236476
  * (13) Barabási, A.L., Pósfai, M.: Network science. Cambridge: Cambridge University Press; (2016). Available from: http://barabasi.com/networksciencebook/.
  * (14) Antoniou, Ioannis E and Tsompa, ET : Statistical analysis of weighted networks. Discrete dynamics in Nature and Society 2008 (2008) doi: 10.1155/2008/375452
  * (15) Barthélemy, M., Barrat, A., Pastor-Satorras, R., Vespignani, A.: Characterization and modeling of weighted networks. Physica A: Statistical Mechanics and its Applications. 346(1-2) 34-43 (2005) doi: 10.1016/j.physa.2004.08.047
  * (16) Travers, J., Milgram, S.: The Small World Problem. Psychology Today 1(1) 61–67 (May 1967)
  * (17) Milgram, S.: An Experimental Study of the Small World Problem. Sociometry. American Sociological Association 32(4) 425–443 (Dec. 1969) doi: 10.2307/2786545
  * (18) Erdő s, P., Rényi, A.: On Random Graphs I. Publicationes Mathematicae 6 290–297 (1959)
  * (19) Watts, D.J., Strogatz, S.H.: Collective dynamics of small-world networks. Nature 393 440-–442 (1998)
  * (20) Barabási, A. L., Albert, R.: Emergence of scaling in random networks. Science 286. 509–512 (1999)
  * (21) Chen H., Chung W., Xu J., Wang G., Qin Y., Chau M.: Crime Data Mining: A General Framework and Some Examples. IEEE Computer, 37 50–56. IEEE (2004). doi: 10.1109/MC.2004.1297301
  * (22) Bahulkar, A., Szymanski, B. K., Baycik, N. O., Sharkey, T. C.: Community Detection with Edge Augmentation in Criminal Networks. In 2018 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining (ASONAM). 1168–1175 (Aug. 2018). IEEE.
  * (23) Berlusconi, G., Calderoni, F., Parolini, N., Verani, M., Piccardi, C.: Link Prediction in Criminal Networks: A Tool for Criminal Intelligence Analysis. Public Library of Science. PLOS ONE. 11(4) 1–21 (Apr. 2016) doi: 10.1371/journal.pone.0154244
  * (24) Franchetti, L., Sonnino, S.: La Sicilia nel 1876. 1. Barbèra G. (1877)
  * (25) McGloin, J.M.: Policy and intervention considerations of a network analysis of street gangs. Criminology & Public Policy. 4(3), 607–635 (2005) doi: 10.1111/j.1745-9133.2005.00306.x
  * (26) Mastrobuoni, G., Patacchini, E.: Organized crime networks: An application of network analysis techniques to the American Mafia. Review of Network Economics. 11(3) (2012) doi: 10.1515/1446-9022.1324
  * (27) Krebs, V.: Mapping networks of terrorist cells. Connections, 24(3) 43–52. INSNA (2002)
  * (28) Xu, J., Chen, H.: Criminal Network Analysis and Visualization. Communications of the ACM, 48(6) 100–107. ACM (2005) doi: 10.1145/1064830.1064834
  * (29) Calderoni, F., Morselli, C.: Inside Criminal Networks. European Journal on Criminal Policy and Research. 16(1) 69—70 (2010) doi: 10.1007/s10610-010-9118-7
  * (30) Campana, P., Varese. F.: Listening to the wire: criteria and techniques for the quantitative analysis of phone intercepts. Trends in Organized Crime. 15(1) 13–30 (2012). 10.1007/s12117-011-9131-3
  * (31) Castaldo. F. (eds): Messina, arrestati il capo ed i sodali della “Famiglia mafiosa di Mistretta”. In: Grandangolo, il giornale di Agrigento (Jan. 18th 2019) Available at: https://www.grandangoloagrigento.it/mafia/messina-arrestati-il-capo-ed-i-sodali-della-famiglia-mafiosa-di-mistretta
  * (32) Newman, M.E.J., Watts, D.J.: Renormalization group analysis of the small-world network model. Physics Letters A 263(4) 341–346 (1999) doi: 10.1016/S0375-9601(99)00757-4


Learn Quasi-stationary Distributions of Finite State Markov Chain

Zhiqiang Cai111Corresponding author. School of Data Science, City University
of Hong Kong, Tat Chee Ave, Kowloon, Hong Kong SAR. Email:
<EMAIL_ADDRESS>, Ling Lin222School of Mathematics, Sun Yat-sen
University, Guangzhou 510275, China. Email<EMAIL_ADDRESS>, Xiang
Zhou333 School of Data Science and Department of Mathematics, City University
of Hong Kong, Tat Chee Ave, Kowloon, Hong Kong SAR. Email:
<EMAIL_ADDRESS>
## Abstract

We propose a reinforcement learning (RL) approach to compute the expression of
quasi-stationary distribution. Based on the fixed-point formulation of quasi-
stationary distribution, we minimize the KL-divergence of two Markovian path
distributions induced by the candidate distribution and the true target
distribution. To solve this challenging minimization problem by gradient
descent, we apply the reinforcement learning technique by introducing the
reward and value functions. We derive the corresponding policy gradient
theorem and design an actor-critic algorithm to learn the optimal solution and
the value function. The numerical examples of finite state Markov chain are
tested to demonstrate the new method.

Keywords: quasi-stationary distribution, reinforcement learning, KL-
divergence, actor-critic algorithm

##  1\. Introduction

Quasi-stationary distribution (QSD) is the long time statistical behavior of a
stochastic process that will be surely killed when this process is conditioned
to survive [13]. This concept has been widely used applications, such as in
biology and ecology [10, 20], chemical kinetics [14, 17], epidemics [2, 12,
33], medicine [11] and neuroscience [5, 21]. Many works for the rare events in
meta-stable systems also focus on this quasi-stationary distribution [16, 22].
In addition, some new Monte Carlo sampling methods, for instance, the Quasi-
stationary Monte Carlo method [28, 36] also arise by using the QSD instead of
true stationary distribution.

We are interested in the numerical computation of QSD and focus on the finite
state Markov chain in this paper. Mathematically, the quasi-stationary
distribution can be solved as the principal left eigenvector of a sub-
Markovian transition matrix. So, the traditional numerical algebra methods can
be applied to solve the quasi-stationary distribution in finite state space,
for example, the power method [38], the multi-grid method [3] and Arnoldi’s
algorithm [27]. These eigenvector methods can produce the stochastic vector
for QSD, instead of generating samples of QSD.

In search of efficient algorithms for large state space, stochastic approaches
are in favor of either sampling the QSD or computing the expression of QSD,
and these methods can be applied or extended easily to continuous state space
too. A popular approach for sampling quasi-stationary distribution is the
Fleming-Viot stochastic method [24]. The Flemming-Viot method first simulates
$N$ particles independently. When any one of the particles falls into the
absorbing state and gets killed, a new particle is uniformly selected from the
remaining $N-1$ surviving particles to replace the dead one, and the
simulation continues. When the time and $N$ tend to infinity, the particles’
empirical distribution can converge to the quasi-stationary distribution.

In [1, 4, 15], the authors proposed to recursively update the expression of
QSD at each iteration based on the empirical distribution of a single-particle
simulation. It is shown [4] that the convergence rate can be $O(n^{-1/2})$,
where $n$ is the iteration number. This method is later improved in [7, 39] by
applying the stochastic approximation method [19] and the Polyak-Ruppert
averaging technique [29]. These improved algorithms have a choice of flexible
step size but needs a projection operator onto probability simplex, which
carries some extra computational overhead increasing with the number of
states. [36] extends the algorithm to the diffusion process.

In this paper, we focus on how to compute the expression of the quasi-
stationary distribution, which is denoted by $\alpha(x)$ on a metric space
$\mathcal{E}$. If $\mathcal{E}$ is finite, $\alpha$ is a probability vector
and if $\mathcal{E}$ is a domain in $\mathbb{R}^{d}$, then $\alpha$ is a
probability density function on $\mathcal{E}$. We assume $\alpha$ can be
numerically represented in parametric form $\alpha_{\theta}$,
$\theta\in\Theta$. This family $\\{\alpha_{\theta}\\}$ can be in tabular form
or any neural network. Then the problem of finding the QSD $\alpha$ becomes
how to compute the optimal parameter $\theta$ in $\Theta$. We call this
problem the learning problem for the QSD. In addition, we want to directly
learn the QSD, not to use the distribution family $\\{\alpha_{\theta}\\}$ to
fit the simulated samples generated by other traditional simulation methods.

Our minimization problem for QSD is similar to the variational inference
(VI)[8], which minimizes an objective functional measuring the distance
between the target and candidate distributions. However, unlike the mainstream
VI methods such as evidence lower bound (ELBO) technique[18] or particle-
based[23], flow-based methods [31], our approach is based on recent important
progresses from reinforcement learning(RL) method[35], particularly the policy
gradient method and actor-critic algorithm. We first regard the learning
process of the quasi-stationary distribution as the interaction with the
environment which is constructed by the property of QSD. Reinforcement
learning has recently shown tremendous advancements and remarkable successes
in applications (e.g.[26, 34, 30]). The RL framework provides an innovative
and powerful modelling and computation approach for many scientific computing
problems.

The essential question is how to formulate the QSD problem as an RL problem.
Firstly, for the sub-Markovian kernel $K$ of a Markov process, we can define a
Markovian kernel $K_{\alpha}$ on $\mathcal{E}$ (see Definition 2.1) and then
the QSD is defined by the equation $\alpha=\alpha K_{\alpha}$, which equals
$\alpha$ as the initial distribution and the distribution after one step.
Secondly, we consider an optimal $\alpha$ (in our parametric family of
distribution) to minimize the Kullback–Leibler divergence (i.e., relative
entropy) of two path distributions, denoted by $\mathbb{P}$ and $\mathbb{Q}$,
associated with two Markovian kernels $K_{\alpha}$ and $K_{\beta}$ where
$\beta:=\alpha K_{\alpha}$. Thirdly, inspired by the recent work [32] of using
RL for rare events sampling problems, we transform the minimization of KL
divergence between $\mathbb{P}$ and $\mathbb{Q}$ into the maximization of a
time-averaged reward function and define the corresponding value function
$V(x)$ at each state $x$. This completes our modelling of RL for the quasi-
stationary distribution problem. Lastly, we derive the policy gradient theorem
(Theorem 2) to compute the gradient w.r.t. $\theta$ of the averaged reward for
the learning dynamic for the averaged reward. This is known as the “actor”
part. The “critic” part is to learn the value function $V$ in its parametric
form $V_{\psi}$. The actor-critic algorithm uses the stochastic gradient
descent to train the parameter $\theta$ for the action $\alpha_{\theta}$ and
the parameter $\psi$ for the value function $V_{\psi}$. See Algorithm 1.

Our contribution is that we first devise a method to transform the QSD problem
into the RL problem. Similar to [32], our paper also uses the KL-divergence to
define the RL problem. However, our paper fully adapts the unique property of
QSD that is a fixed point problem $\alpha=\alpha K_{\alpha}$ to define the RL
problem.

Our learning method allows the flexible parametrization of the distributions
and uses the stochastic gradient method to train the optimal distribution. It
is easy to implement optimization with scale up to large state spaces. The
numerical examples we tested have shown our methods converge faster than the
other existing methods [15, 7].

Finally, we remark that our method works very well for QSD of the strict sub-
Markovian kernel $K$, but not applicable to compute the invariant distribution
when $K$ is Markovian. This is because we transform the problem into the
variational problem between two Markovian kernels $K_{\alpha}$ and $K_{\beta}$
(where $\beta=\alpha K_{\alpha}$). Note
$K_{\alpha}(x,y)=K(x,y)+(1-K(x,\mathcal{E}))\alpha(y)$ (Definition 2.1) and
our method is based on the fact that $\alpha=\beta$ if and only if
$K_{\alpha}=K_{\beta}$. If $K$ is Markovian kernel, then $K_{\alpha}\equiv K$
for any $\alpha$, and our method can not work. So $K(x,\mathcal{E})$ has to be
strictly less than 1 for some $x\in\mathcal{E}$.

This paper is organized as follows. Section 2 is a short review of the quasi-
stationary distribution and some basic simulation methods of QSD. In Section
3, we first formulate the reinforcement learning problem by KL-divergence and
derive the policy gradient theorem (Theorem 2). Using the above formulation,
we then develop the actor-critic algorithm to estimate the quasi-stationary
distribution. In Section 4, the efficiency of our algorithms is illustrated by
four examples compared with the simulation methods in [39].

##  2\. Problem Setup and Review

###  2.1. Quasi-stationary Distribution

We start with an abstract setting. Let $\mathcal{E}$ be a finite state
equipped with the Borel $\sigma$-field $\mathcal{B}(\mathcal{E})$, and let
$\mathcal{P}(\mathcal{E})$ be the space of probabilities over $\mathcal{E}$. A
sub-Markovian kernel on $\mathcal{E}$ is defined as a map
$K:\mathcal{E}\times\mathcal{B}(\mathcal{E})\mapsto[0,1]$ such that for all
$x\in\mathcal{E},A\mapsto K(x,A)$ is a nonzero measure with
$K(x,\mathcal{E})\leq 1$ and for all $A\in\mathcal{B}(\mathcal{E}),x\mapsto
K(x,A)$ is measurable. Particularly, if $K(x,\mathcal{E})=1$ for all
$x\in\mathcal{E}$, then $K$ is called a Markovian kernel. Throughout the
paper, assume that $K$ is strictly sub-Markovian, i.e., $K(x,\mathcal{E})<1$
for some $x$.

Let $X_{t}$ be a Markov chain with values in
$\mathcal{E}\cup\left\\{\partial\right\\}$ where $\partial\notin\mathcal{E}$
denotes an absorbing state. We define the extinction time

$\tau:=\inf\left\\{t>0:\quad X_{t}=\partial\right\\}.$

We define the quasi-stationary distribution(QSD) $\alpha$ as the long time
limit of the conditional distribution, if there exists a probability
distribution $\nu$ on $\mathcal{E}$ such that:

$\alpha(A):=\lim_{t\rightarrow\infty}P_{\nu}\left(X_{t}\in
A\mid\tau>t\right),\qquad A\in\mathcal{B}(\mathcal{E}).$ (1)

where $P_{\nu}$ refers to the probability distribution of $X_{t}$ associated
with the initial distribution $\nu$ on $\mathcal{E}$. Such a conditional
distribution well describes the behavior of the process before extinction and
it is easy to see that $\alpha$ satisfies the following fixed point problem

${P}_{\alpha}\left(X_{t}\in A\mid\tau>t\right)=\alpha(A)$ (2)

where $P_{\alpha}$ refers to the probability distribution of $X_{t}$
associated with the initial distribution $\alpha$ on $\mathcal{E}$. (2) is
equivalent to the following stationary condition that

$\alpha=\frac{\alpha K}{\alpha K\mathbf{1}},\quad\mbox{or
}\alpha(y)=\frac{\sum_{x}\alpha(x)K(x,y)}{\sum_{x}\alpha(x)K(x,\mathcal{E})}$
(3)

where $\alpha$ is a row vector and $\mathbf{1}$ denotes the column vector with
all entries being one and

$K(x,\mathcal{E})=\sum_{x^{\prime}\in\mathcal{E}}K(x,x^{\prime}).$

For any sub-Markovian kernel $K$, we can associate with $K$ a Markovian kernel
$\tilde{K}$ on $\mathcal{E}\cup\\{\partial\\}$ defined by

$\left\\{\begin{array}[]{l}\tilde{K}(x,A)=K(x,A)\\\
\tilde{K}(x,\\{\partial\\})=1-K(x,\mathcal{E})\\\
\tilde{K}(\partial,\\{\partial\\})=1.\end{array}\right.$

for all $x\in\mathcal{E},A\in\mathcal{B}(\mathcal{E})$. The kernel $\tilde{K}$
can be understood as the Markovian transition kernel of the Markov chain
$(X_{t})$ on $\mathcal{E}\cup\\{\partial\\}$ whose transitions in
$\mathcal{E}$ is specified by $K$, but is “killed” forever once it leaves
$\mathcal{E}$.

In this paper, we assume $\mathcal{E}$ is a finite state space and the process
in consideration has a unique QSD. Assume that $K$ is irreducible, then
existence and uniqueness of the quasi-stationary distribution can be obtained
by the Perron-Frobenius theorem [25].

An important Markovian kernel is the following $K_{\alpha}$ which is defined
on $\mathcal{E}$ only and has a “regenerative probability” $\alpha$.

######  Definition 2.1.

For any given $\alpha\in\mathcal{P}(\mathcal{E})$ and a sub-Markovian kernel
$K$ on $\mathcal{E}$, we define $K_{\alpha}$, a Markovian kernel on
$\mathcal{E}$, as follows

$K_{\alpha}(x,A):=K(x,A)+\left(1-K(x,\mathcal{E})\right)\alpha(A)$ (4)

for all $x\in\mathcal{E}$ and $A\in\mathcal{B}(\mathcal{E})$.

$K_{\alpha}$ is a Markovian kernel because $K_{\alpha}(x,\mathcal{E})=1$. It
is easy to sample $X_{t+1}\sim K_{\alpha}(X_{t},\cdot)$ from any state
$X_{t}\in\mathcal{E}$: Run the transition as normal by using $\tilde{K}$ to
have a next state denoted by $Y$, then $X_{t+1}=Y$ if $Y\in\mathcal{E}$,
otherwise, sample $X_{t+1}$ from $\alpha$.

We know that $\alpha$ is the quasi-stationary distribution of $K$ if and only
if it is the stationary distribution of $K_{\alpha}$, i.e.,

$\alpha=\alpha K_{\alpha}.$ (5)

It is easy to see $\alpha=\beta$ if and only if $K_{\alpha}=K_{\beta}$ for any
two distributions $\alpha$ and $\beta$. Also, for every $\alpha^{\prime}$,
$K_{\alpha^{\prime}}$ has a unique invariant probability denoted by
$\Gamma(\alpha^{\prime})$. Then
$\alpha^{\prime}\mapsto\Gamma(\alpha^{\prime})$ is continuous in
$\mathcal{P}(\mathcal{E})$ (i.e. for the topology of weak convergence) and
there exists $\alpha\in\mathcal{P}(\mathcal{E})$ such that
$\alpha=\Gamma(\alpha)$ or, equivalently, $\alpha$ is a QSD for $K$.

###  2.2. Review of simulation methods for quasi-stationary distribution

According to the above subsection, the QSD $\alpha$ satisfies the fixed point
problem

$\alpha=\Gamma(\alpha),$ (6)

where $\Gamma(\alpha)$ is the stationary distribution of $K_{\alpha}$ on
$\mathcal{E}$. In general, (6) can be solved recursively by
$\alpha_{n+1}\leftarrow\Gamma(\alpha_{n})$.

The Fleming–Viot (FV) method [24] evolves $N$ particles independently of each
other as a Markov process associated with the transition kernel $K_{\alpha}$
until one succeeds in jumping to the absorbing state $\partial$. At that time,
this killed particle is immediately reset to $\mathcal{E}$ as an initial state
uniformly chosen from one of the remaining $N-1$ particles. The QSD $\alpha$
is approximated by the empirical distribution of the $N$ particles in total
and these particles can be regarded as samples from the quasi-stationary
distribution $\alpha$ like the MCMC method.

[6] proposed a simulation method by only using one particle at each iteration
to update $\alpha$. At iteration $n$, given an
$\alpha_{n}\in\mathcal{P}(\mathcal{E})$, one can run the discrete-time Markov
chain $X^{(n+1)}$ as normal on $\partial\cup\mathcal{E}$ with initial
$X^{(n+1)}_{0}\sim\alpha_{n}$, then $\alpha_{n+1}$ is computed as the
following weighted average of empirical distributions:

$\begin{split}\alpha_{n+1}(x):=\alpha_{n}(x)+\frac{1}{n+1}\sum_{k=0}^{\tau^{(n+1)}-1}\frac{I\left(X_{k}^{(n+1)}=x\mid
X_{0}^{(n+1)}\sim\alpha_{n}\right)-\alpha_{n}(x)}{\frac{1}{n+1}\sum_{j=1}^{n+1}\tau^{(j)}}\end{split}$
(7)

where $n\geq 0$ and $I$ is the indicator function,
$\tau^{(j)}=\min\left\\{k\geq 0\mid X_{k}^{(j)}\in\partial\right\\}$ is the
first extinction time for the process $X^{(j)}$. This iterative scheme has the
convergence rate $O(\frac{1}{\sqrt{n}})$.

In [7, 39], the above method is extended to the stochastic approximations
framework

$\alpha_{n+1}(x)=\Theta_{H}\left[\alpha_{n}+\epsilon_{n}\sum_{k=0}^{\tau^{(n+1)}-1}\left(I\left(X_{k}^{(n+1)}=x|X_{0}^{(n+1)}\sim\alpha_{n}\right)-\alpha_{n}(x)\right)\right]$
(8)

where $\Theta_{H}$ denotes the $L_{2}$ projection into the probability simplex
and $\epsilon_{n}$ is the step size satisfying $\sum{\epsilon_{n}}=\infty$ and
$\sum{\epsilon^{2}_{n}}<\infty$. Specifically, if
$\epsilon_{n}=O(\frac{1}{n^{r}})$ for $0.5<r<1$, under sufficient condition,
they have
$\sqrt{n^{r}}\left(\alpha_{n}-{\alpha}\right)\xrightarrow{\text{d}}\mathcal{N}(0,V)$
for some matrix $V$[7, 39]. If the Polyak-Ruppert averaging technique is
applied to generate

$\nu_{n}:=\frac{1}{n}\sum_{k=1}^{n}\alpha_{k},$ (9)

then the convergence rate of $\nu_{n}\to\alpha$ becomes $\frac{1}{\sqrt{n}}$
[7, 39].

The simulation schemes (7) and (8) need to sample the initial states according
to $\alpha_{n}$ and to add the empirical distribution and $\alpha_{n}$ at each
$x$ pointwisely. So they are suitable for finite state space where $\alpha$ is
a probability vector saved in the tabular form. In (8) there is no need to
record all exit times $\tau^{(j)},j=1,\dots,n$, but the additional projection
operation in (8) is computationally expensive since the cost is $O(m\log m)$
where $m=|\mathcal{E}|$ [9, 37].

##  3\. Learn Quasi-stationary Distribution

We focus on the computation of the expression of the quasi-stationary
distribution. Particularly, when this distribution is parametrized in a
certain way by $\theta$, we can extend the tabular form for finite-state
Markov chain to any flexible form, even in the neural networks for probability
density function in $\mathbb{R}^{d}$. But we do not pursue this representation
and expressivity issue here, and restrict our discussion to finite state space
only to illustrate our main idea first. In finite state space, $\alpha(x)$ for
$x\in\mathcal{E}=\\{1,\dots,m\\}$, can be simply described as a softmax
function with $m-1$ parameter $\theta_{i}:\alpha(i)\propto
e^{\theta_{i}},1\leq i\leq m-1$, ($\theta_{m}=0$). This introduces no
representation error. For the generalization to continuous space $\mathcal{E}$
in jump and diffusion processes, or even for a huge finite state space, a good
representation of $\alpha_{\theta}(x)$ is important in practice.

In this section, we shall formulate our QSD problem in terms of reinforcement
learning (RL) so that the problem of seeking optimal parameters becomes a
policy optimization problem. We derive the policy gradient theorem to
construct a gradient descent method for the optimal parameter. We then show
how to design the actor-critic algorithms based on stochastic optimization.

###  3.1. Formulation of RL and Policy Gradient Theorem

Before introducing the RL method of our QSD problem, we develop a general
formulation in terms of RL by introducing the KL-divergence between two path
distributions.

Let $P_{\theta}$ and $Q_{\theta}$ be two families of Markovian kernels on
$\mathcal{E}$ in parametric forms with the same set of parameters
$\theta\in\Theta$. Assume both $P_{\theta}$ and $Q_{\theta}$ are ergodic for
any $\theta$. Let $T>0$ and denote a path up to time $T$ by
$\omega_{0}^{T}=(X_{0},X_{1},\ldots,X_{T})\in\mathcal{E}^{T+1}$. Define the
path distributions under the Markov chain kernel $P_{\theta}$ and
$Q_{\theta}$, respectively:

$\mathbb{P}_{\theta}(\omega_{0}^{T}):=\prod_{t=1}^{T}P_{\theta}(X_{t}\mid
X_{t-1}),\quad\mathbb{Q}_{\theta}(\omega_{0}^{T}):=\prod_{t=1}^{T}Q_{\theta}(X_{t}\mid
X_{t-1}).$ (10)

Define the KL divergence from $\mathbb{P}_{\theta}$ to $\mathbb{Q}_{\theta}$
on $\mathcal{E}^{T+1}$

$\mathsf{D}_{KL}(\mathbb{P}_{\theta}\mid\mathbb{Q}_{\theta}):=\sum_{\omega_{0}^{T}}{\mathbb{P}_{\theta}(\omega_{0}^{T})\ln{\frac{\mathbb{P}_{\theta}(\omega_{0}^{T})}{\mathbb{Q}_{\theta}(\omega_{0}^{T})}}}=-\mathbf{E}_{P_{\theta}}\sum_{t=1}^{T}R_{\theta}(X_{t-1},X_{t}),$
(11)

where the expectation $\mathbf{E}_{P_{\theta}}$ is for the path
$(X_{0},X_{1},\ldots,X_{T})$ generated by the transition kernel $P_{\theta}$
and

$R_{\theta}(X_{t-1},X_{t}):=-\ln{\frac{P_{\theta}(X_{t}\mid
X_{t-1})}{Q_{\theta}(X_{t}\mid X_{t-1})}}.$ (12)

is called the (one-step) reward.

Define the  average reward  $r(\theta)$ as the time averaged negative KL
divergence in the limit of $T\to\infty$:

$\displaystyle r(\theta)$
$\displaystyle:=-\lim_{T\to\infty}{\frac{1}{T}\mathsf{D}_{KL}(\mathbb{P}_{\theta}\mid\mathbb{Q}_{\theta})}=-\lim_{T\to\infty}\frac{1}{T}\mathbf{E}_{P_{\theta}}\sum_{t=1}^{T}R_{\theta}(X_{t-1},X_{t}).$
(13)

Due to ergodicity of $P_{\theta}$,
$r(\theta)=\sum_{x_{0},x_{1}}R_{\theta}(x_{0},x_{1})P_{\theta}(x_{1}|x_{0})\mu_{\theta}(x_{0})$
where $\mu_{\theta}$ is the invariant measure of $P_{\theta}$. $r(\theta)$ is
independent of initial state $X_{0}$. Obviously $r(\theta)\leq 0$ for any
$\theta$.

######  Property 1.

The following are equivalent:

  1. (1)

$r(\theta)$ reaches its maximal value $0$ at $\theta^{*}$.

  2. (2)

$\mathbb{P}_{\theta^{*}}=\mathbb{Q}_{\theta^{*}}$ in
$\mathcal{P}(\mathcal{E}^{T+1})$ for any $T>0$.

  3. (3)

$P_{\theta^{*}}=Q_{\theta^{*}}$.

  4. (4)

$R_{\theta^{*}}\equiv 0$.

###### Proof.

We only need to show $(1)\Longrightarrow(3)$. It is easy to see

$r(\theta)=-\sum_{x_{0}}~{}\mathsf{D}_{KL}(P_{\theta}(\cdot|x_{0})~{}|~{}Q_{\theta}(\cdot|x_{0}))~{}\mu_{\theta}(x_{0}).$

If $r(\theta)=0$, since $\mu_{\theta}>0$, then

$\mathsf{D}_{KL}(P_{\theta}(\cdot|x_{0})~{}|~{}Q_{\theta}(\cdot|x_{0}))=0~{}~{}\quad\forall
x_{0}.$

So we have $P_{\theta}=Q_{\theta}$.

∎

The above property establishes the relationship between the RL problem and QSD
problem.

We show our theoretic main result below as the foundation of our algorithm to
be developed later. This theorem can be regarded as one type of the policy
gradient theorem for policy gradient method in reinforcement learning [35].

Define the value function ([35] Chapter 13):

$V\left(x\right):=\lim_{T\rightarrow\infty}\sum_{t=1}^{T}\mathbf{E}_{P_{\theta}}\left[R_{\theta}(X_{t-1},X_{t})-r(\theta)\mid
X_{0}=x\right].$ (14)

Certainly, $V$ also depends on $\theta$, though we do not write $\theta$
explicitly.

######  Theorem 2 (policy gradient theorem).

We have the following two properties

  1. (1)

At any $\theta$, for any $x\in\mathcal{E}$, the following Bellman-type
equation holds for the value function $V$ and the average reward $r(\theta)$:

$V(x)=\mathbf{E}_{Y\sim P_{\theta}(\cdot\mid
x)}\left[V(Y)+R_{\theta}(x,Y)-r(\theta)\right].$ (15)

  2. (2)

The gradient of the average reward ${r}(\theta)$ is

$\displaystyle\nabla_{\theta}r(\theta)$
$\displaystyle=\mathbf{E}\left[\nabla_{\theta}\ln Q_{\theta}(Y\mid X)\right]+$
$\displaystyle\qquad\mathbf{E}\left[\bigg{(}V(Y)-V(X)+R_{\theta}(X,Y)-r(\theta)\bigg{)}\nabla_{\theta}\ln
P_{\theta}(Y\mid X)\right],$ (16)

where the expectations are for the joint distribution
$(X,Y)\sim\mu_{\theta}(x)P_{\theta}(y\mid x)$ where $\mu_{\theta}$ is the
stationary measure of $P_{\theta}$.

###### Proof.

We shall prove the Bellman equation first and then we use the Bellman equation
to derive the gradient of the average reward ${r}(\theta)$. For any
$x_{0}\in\mathcal{E}$, by writing $\omega_{0}^{T}=(x_{0},\ldots,x_{T})$ and
defining

$\Delta R_{\theta}(\omega_{0}^{T})=\sum_{t=1}^{T}(R(x_{t-1},x_{t})-r(\theta))$

we have

$\displaystyle V\left(x_{0}\right)$
$\displaystyle=\lim_{T\rightarrow\infty}\mathbf{E}_{P_{\theta}}\left[\Delta
R_{\theta}(\omega_{0}^{T})\mid X_{0}=x\right]$
$\displaystyle=\lim_{T\to\infty}\sum_{x_{2},\ldots,x_{T}}\sum_{x_{1}}\left(\left(\prod_{t=2}^{T}P_{\theta}(x_{t}\mid
x_{t-1})\right)P_{\theta}(x_{1}\mid x_{0})\Delta R(\omega_{0}^{T})\right)$
$\displaystyle=\lim_{T\to\infty}\sum_{x_{1}}\left(P_{\theta}(x_{1}\mid
x_{0})\sum_{x_{2},\ldots,x_{T}}\left(\prod_{t=2}^{T}P_{\theta}(x_{t}\mid
x_{t-1})\left[\Delta R(\omega_{1}^{T})+\Delta
R(\omega_{0}^{1})\right]\right)\right)$
$\displaystyle=\sum_{x_{1}}\left(P_{\theta}(x_{1}\mid
x_{0})\left(\lim_{T\to\infty}\left[\sum_{x_{2},\ldots,x_{T}}\prod_{t=2}^{T}P_{\theta}(x_{t}\mid
x_{t-1})\Delta R(\omega_{1}^{T})\right]+\Delta
R(\omega_{0}^{1})\right)\right)$
$\displaystyle=\sum_{x_{1}}P_{\theta}(x_{1}\mid
x_{0})\left[V(x_{1})+R_{\theta}(x_{0},x_{1})\right]-r(\theta),$ (17)

which proves (15), i.e.,

$r(\theta)=\mathbf{E}_{Y\sim P_{\theta}(\cdot\mid
x)}\left[V(Y)+R_{\theta}(x,Y)-V(x)\right],\qquad\forall x\in\mathcal{E}.$

Next, we compute the gradient of $r(\theta)$. By the trivial equality

$\sum_{x_{1}}P_{\theta}(x_{1}\mid x_{0})\nabla_{\theta}\ln
P_{\theta}(x_{1}\mid x_{0})=\nabla_{\theta}\sum_{x_{1}}P_{\theta}(x_{1}\mid
x_{0})=0,$ (18)

and the definition (12), we can write the gradient of $r(\theta)$ as follows

$\displaystyle\nabla_{\theta}r(\theta)=$
$\displaystyle\sum_{y}\nabla_{\theta}P_{\theta}(y\mid
x)\left[V(y)+R_{\theta}(x,y)-V(x)\right]$
$\displaystyle+\sum_{y}P_{\theta}\left(y\mid
x\right)\left[\nabla_{\theta}V(y)-\nabla_{\theta}V(x)+\nabla_{\theta}\ln{Q_{\theta}(y\mid
x)}\right].$

We here keep the term $V(x)$ in the first line, even though it has no
contribution here (in fact, to add any constant to $V(x)$ is also fine). Since
this equation holds for all states $x$ on the right-hand side, we take the
expectation w.r.t. $\mu_{\theta}$, the stationary distribution of
$P_{\theta}$. So, we have

$\displaystyle\nabla_{\theta}r(\theta)=$
$\displaystyle\sum_{x,y}\mu_{\theta}(x)\nabla_{\theta}P_{\theta}\left(y\mid
x\right)\left[V(y)+R_{\theta}(x,y)-V(x)\right]$
$\displaystyle+\sum_{x,y}\mu_{\theta}(x)P_{\theta}\left(y\mid
x\right)\left[\nabla_{\theta}V(y)-\nabla_{\theta}V(x)+\nabla_{\theta}\ln{Q_{\theta}(y\mid
x)}\right]$ $\displaystyle=$
$\displaystyle\sum_{x,y}\mu_{\theta}(x)\nabla_{\theta}P_{\theta}\left(y\mid
x\right)\left[V(y)+R_{\theta}(x,y)-V(x)\right]$
$\displaystyle+\sum_{y}\mu_{\theta}(y)\nabla_{\theta}V(y)-\sum_{x}\mu_{\theta}(x)\nabla_{\theta}V(x)+\sum_{x,y}\mu_{\theta}(x)P_{\theta}\left(y\mid
x\right)\nabla_{\theta}\ln{Q_{\theta}(y\mid x)}$ $\displaystyle=$
$\displaystyle\sum_{x,y}\mu_{\theta}(x)P_{\theta}\left(y\mid
x\right)\bigg{[}V(y)+R_{\theta}(x,y)-V(x)\bigg{]}~{}\nabla_{\theta}\ln{P_{\theta}\left(y\mid
x\right)}$ $\displaystyle+\sum_{x,y}\mu_{\theta}(x)P_{\theta}\left(y\mid
x\right)\nabla_{\theta}\ln{Q_{\theta}(y\mid x)}.$

In fact, we can add any constant number $b$ (independent of $x$ and $y$)
inside the squared bracket of the last line without changing the equality, due
to the following fact similar to (18):
$\sum_{x,y}\mu_{\theta}(x)\nabla_{\theta}P_{\theta}\left(y\mid
x\right)=\sum_{y}\mu_{\theta}(y)\nabla_{\theta}\sum_{x}P_{\theta}\left(x\mid
y\right)=0$. (2) is a special case of $b=r(\theta)$. ∎

######  Remark 1.

As shown in the proof, (2) holds if $r(\theta)$ at the right-hand side is
replaced by any constant number $b$. $b=r(\theta)$ is a good choice to reduce
the variance since $r(\theta)$ can be regarded as the expectation of
$R_{\theta}$.

######  Remark 2.

If $P_{\theta}=Q_{\theta}$, then the first term of (2) vanishes due to (18)
and the second term of (2) vanishes due to (15).

######  Remark 3.

The name of “policy” here refers to the role of $\theta$ as the policy for
decision makers to improve the reward $r(\theta)$.

###  3.2. Learn QSD

Now we discuss how to connect the QSD with the results in the previous
subsection. In view of equation (5), we introduce $\beta:=\alpha K_{\alpha}$
as the one-step distribution if starting from the initial $\alpha$, i.e.,

$\beta(y):=\sum_{x\in\mathcal{E}}{\alpha(x)K_{\alpha}(x,y)},~{}~{}\quad\forall
y$ (19)

By (5), $\alpha$ is a QSD if and only if $\beta=\alpha$. However, we do not
directly compare these two distributions $\alpha$ and $\beta$. Instead, we
consider their Markovian kernels induced by (4): $K_{\alpha}$ and $K_{\beta}$.
Our approach is to consider the KL divergence similar to (11) between two
kernels $K_{\alpha}$ and $K_{\beta}$ since $\alpha=\beta$ if and only if
$K_{\alpha}=K_{\beta}$. In this way, one can view $K_{\alpha}$ and $K_{\beta}$
(note $\beta=\alpha K_{\alpha}$) as two transition matrices $P_{\theta}$ and
$Q_{\theta}$ in the previous section, in which the parameter $\theta$ here is
in fact the distribution $\alpha$.

To have a further representation of the distribution $\alpha$, which is a
(probability mass) function on $\mathcal{E}$, we propose a parametrized family
for $\alpha$ in the form $\alpha_{\theta}$ where $\theta$ is a generic
parameter. In the simplest case, $\alpha_{\theta}$ takes the so-called soft-
max form $\alpha_{\theta}(i)=\frac{e^{\theta_{i}}}{\sum_{j\geq
1}e^{\theta_{j}}}$ if $\mathcal{E}=\\{1,\ldots,N\\}$ for
$\theta=(\theta_{1},\ldots,\theta_{N-1},\theta_{N}\equiv 0).$ This
parametrization represents $\alpha$ without any approximation error for finite
state space and the effective space of $\theta$ is just $\mathbb{R}^{N-1}$.
For certain problems, particularly with large state space, if one has some
prior knowledge about the structure of the function $\alpha$ on $\mathcal{E}$,
one might propose other parametric forms of $\alpha_{\theta}$ with the
dimension of $\theta$ less than the cardinality $|\mathcal{E}|$ to improve the
efficiency, although the extra representation error in this way has to be
introduced.

For any given $\alpha_{\theta}\in\mathcal{P}(\mathcal{E})$, the corresponding
Markovian kernel $K_{\alpha_{\theta}}$ is then defined in (4) and
$\beta_{\theta}=\alpha_{\theta}K_{\alpha_{\theta}}i$ is defined by (19).
$K_{\beta_{\theta}}$ is like-wisely defined by (4) again. To use the
formulation in Section 3.1, we choose $P_{\theta}=K_{\alpha_{\theta}}$ and
$Q_{\theta}=K_{\beta_{\theta}}$. Define the objective function as before:

$\displaystyle{r}(\theta)$
$\displaystyle:=-\lim_{T\to\infty}{\frac{1}{T}\mathsf{D}_{KL}(\mathbb{P}_{\theta}\mid\mathbb{Q}_{\theta})}=-\lim_{T\to\infty}\frac{1}{T}\mathbf{E}_{P_{\theta}}\sum_{t=1}^{T}{R}_{\theta}(X_{t-1},X_{t}).$

where

${R}_{\theta}(x,y)=-\ln{\frac{K_{\alpha_{\theta}}(x,y)}{K_{\beta_{\theta}}(x,y)}}.$

The value function $V(x)$ is defined like-wisely. Theorem 2 now gives the
expression of gradient

$\begin{split}\nabla_{\theta}r(\theta)&=\mathbf{E}[\big{(}R_{\theta}(X,Y)-r(\theta)+V(Y)-V(X)\big{)}\nabla_{\theta}\ln
K_{\alpha_{\theta}}(X,Y)\\\
&\qquad+\nabla_{\theta}\ln{K_{\beta_{\theta}}(X,Y)}]\end{split}$ (20)

where $(X,Y)\sim\mu_{\theta}(x)K_{\alpha_{\theta}}(x,y)$ where $\mu_{\theta}$
is the stationary measure of $K_{\alpha_{\theta}}$.

The optimal $\theta^{*}$ for the QSD $\alpha_{\theta}$ is to maximize
$r(\theta)$ and this can be solved by the gradient descent algorithm

$\theta_{t+1}=\theta_{t}+\eta_{t}^{\theta}\nabla_{\theta}r(\theta_{t}).$ (21)

where $\eta_{t}^{\theta}>0$ is the step size. In practice, the stochastic
gradient is applied

$\nabla_{\theta}r(\theta_{t})\approx{\nabla_{\theta}\ln
K_{\alpha_{\theta}}(X_{t},X_{t+1})\times\delta(X_{t},X_{t+1})+\nabla_{\theta}\ln{K_{\beta_{\theta}}(X_{t},X_{t+1})}}$

where $X_{t},X_{t+1}$ are sampled based on the Markovian kernel
$K_{\alpha_{\theta}}$ ( see Algorithm (1)) and the differential temporal (TD)
error $\delta_{t}$ is

$\delta_{t}=\delta(X_{t},X_{t+1})=R_{\theta}(X_{t},X_{t+1})-r(\theta_{t})+V(X_{t+1})-V(X_{t}).$
(22)

Next, we need to address a remaining issue to address: how to compute the
value function $V$ and $r(\theta_{t})$ in the TD error (22). Besides, we also
need to show the details of computing $\nabla_{\theta}K_{\alpha_{\theta}}$ and
$\nabla_{\theta}K_{\beta_{\theta}}$.

###  3.3. Actor-Critic Algorithm

With the stochastic gradient method (21), we can obtain the optimal policy
$\theta^{*}$. We refer to (21) as the learning dynamics for the policy and it
is generally known as actor. To calculate the value function $V$ appearing in
$\nabla r(\theta)$, we need to have a new learning dynamics, which is called
critic. Then the overall policy-gradient method is termed as the actor-critic
method.

We start with the Bellman equation (15) for the value function and consider
the mean-square-error loss

$\operatorname{MSE}[V]=\frac{1}{2}\sum_{x}\nu(x)\left(\sum_{y}K_{\alpha_{\theta}}(x,y)\left[V(y)+R_{\theta}(x,y)-r(\theta)\right]-V(x)\right)^{2}$

where $\nu$ is any distribution supported on $\mathcal{E}$.
$\operatorname{MSE}[V]=0$ if and only if $V$ satisfies the Bellman equation
(15), i.e. $V$ is the value function. To learn $V$, we introduce the function
approximation for the value function, $V_{\psi}$, with the parameter $\psi$
and consider to minimize

$\operatorname{MSE}(\psi)=\frac{1}{2}\sum_{x}\nu(x)\left(\sum_{y}K_{\alpha_{\theta}}(x,y)\left[V(y)+R_{\theta}(x,y)-r(\theta)\right]-V_{\psi}(x)\right)^{2}$

by the semi-gradient method ([35], Chapter 9):

$\displaystyle\nabla_{\psi}\operatorname{MSE}(\psi)$
$\displaystyle=-\sum_{x,y}\nu(x)K_{\alpha_{\theta}}(x,y)\left[V(y)+R_{\theta}(x,y)-r(\theta)-V_{\psi}(x)\right]\nabla_{\psi}V_{\psi}(x)$
$\displaystyle\approx-\sum_{x,y}\nu(x)K_{\alpha_{\theta}}(x,y)\left[V_{\psi}(y)+R_{\theta}(x,y)-r(\theta)-V_{\psi}(x)\right]\nabla_{\psi}V_{\psi}(x)$

Here the term $V(y)$ is frozen first and then approximated by $V_{\psi}$ since
it could be treated as a prior guess of the value function for the future
state.

Then for the gradient descent iteration
$\psi_{t+1}=\psi_{t}-\eta_{t}^{\psi}\nabla_{\psi}\operatorname{MSE}_{V}(\psi_{t})$
where $\eta_{t}^{\psi}$ is the step size, we can have the stochastic gradient
iteration

$\psi_{t+1}=\psi_{t}+\eta_{t}^{\psi}~{}\delta(X_{t},X_{t+1})~{}\nabla_{\psi}V_{\psi_{t}}(X_{t})$
(23)

where the differential temporal (TD) error $\delta$ defined above in (22):

$\delta_{t}=\delta(X_{t},X_{t+1})=R_{\theta_{t}}(X_{t},X_{t+1})-r(\theta_{t})+V_{\psi_{t}}(X_{t+1})-V_{\psi_{t}}(X_{t}).$

Here for simplicity, $(X_{t},X_{t+1})$ are the same samples as in the actor
method for $\theta_{t}$. This means that the distribution $\nu$ above is
chosen as $\mu$ used for the gradient $\nabla_{\theta}r(\theta)$.

Next, we consider the calculation of the reward $r(\theta)$. Again by the
Bellman equation (15)

$\sum_{x}\mu(x)\sum_{y}K_{\alpha_{\theta}}(x,y)(R_{\theta}(x,y)-r(\theta)+V(y)-V(x))=0$

Let $r_{t}$ be the estimate of the reward $r(\theta_{t})$ at time $t$. We can
update our estimate of the reward every time a transition occurs as

$r_{t+1}=r_{t}+\eta_{t}^{r}\times\delta_{t}$ (24)

where $\delta_{t}$ is the TD error before

$\delta_{t}=\delta(X_{t},X_{t+1})=R_{\theta_{t}}(X_{t},X_{t+1})-r_{t}+V_{\psi_{t}}(X_{t+1})-V_{\psi_{t}}(X_{t}).$

In conclusion, (21)(23)(24) together consist of the actor-critic algorithm,
which is summarized in Algorithm 1. We remark that Algorithm 1 can be easily
adapted to use the mini-batch gradient method where several copies of
$(X_{t},X_{t+1})$ are sampled and the average is used to update the
parameters. The stationary distribution $\mu_{\theta}$ of
$K_{\alpha_{\theta}}$ is sampled by running the corresponding Markov chain for
several steps with “warm start”: the initial for $\theta_{t+1}$ is set as the
final state generated from the previous iteration at $\theta_{t}$. The length
of this “burn-in” period can be set as just one step in practice for
efficiency.

Initialization

$t=0$; $\theta=\theta_{0}$; $\psi=\psi_{0}$; $r_{t}=r_{0}$;

Sample $X_{0}\sim\mu_{\theta_{0}},$ the stationary distribution of
$K_{\alpha_{\theta_{0}}}$

for  _$t=0,1,2,\dots$_ do

Sample $X_{t+1}$ from the transition kernel
$K_{\alpha_{\theta_{t}}}(X_{t},X_{t+1})$

$\delta_{t}=V_{\psi_{t}}(X_{t+1})-V_{\psi_{t}}(X_{t})+R_{\theta_{t}}(X_{t},X_{t+1})-r_{t}$

$\theta_{t+1}=\theta_{t}+\eta_{t}^{\theta}~{}(\delta_{t}\nabla_{\theta}\ln
K_{\alpha_{\theta_{t}}}(X_{t},X_{t+1})+\nabla_{\theta}\ln
K_{\beta_{\theta_{t}}}(X_{t},X_{t+1}))$

$\psi_{t+1}=\psi_{t}+\eta_{t}^{\psi}~{}\delta_{t}\nabla_{\psi}V_{\psi_{t}}(X_{t})$

$r_{t+1}=r_{t}+\eta_{t}^{r}~{}\delta_{t}$

$X_{t}\sim\mu_{\theta_{t+1}}$ the stationary distribution of
$K_{\alpha_{\theta_{t+1}}}$

$t=t+1$

Algorithm 1 (ac-$\alpha$ method) Actor Critic algorithm for quasi-stationary
distribution $\alpha_{\theta}$

######  Remark 4.

Finally, we remark the computation of
$\nabla_{\theta}\ln{K_{\alpha_{\theta}}}$ and
$\nabla_{\theta}\ln{K_{\beta_{\theta}}}$ in Algorithm 1. The details are shown
in Appendix. We comment that the main computational cost is the function
$K(x,\mathcal{E})$, which has to be pre-computed and stored. If the problem
has some special structure, the function could be approximated in parametric
form. Another special case is our example 2 where
$K(x,\mathcal{E})=0\quad\forall x\in\\{2,3,\dots,N\\}$.

##  4\. Numerical experiment

In this section, we present two examples to demonstrate Algorithm 1. We call
the algorithm (7), (8) and (9) in Section 2.2 used in [7, 39], as Vanilla
Algorithm, Projection Algorithm and Polyak Averaging Algorithm respectively.
Let 0 be the absorbing state and $\mathcal{E}=\\{1,\ldots,N\\}$ are non-
absorbing states, the Markov transition matrix on $\\{0,\ldots,N\\}$ is
denoted by

$\tilde{K}=\left[\begin{array}[]{cc}1&0\\\ &K\end{array}\right],$

where $K$ is an $N$-by-$N$ sub-Markovian matrix. For Algorithm 1, the
distribution $\alpha_{\theta}$ on $\mathcal{E}$ is always parameterized as

$\alpha_{\theta}=\frac{1}{e^{\theta_{1}}+\ldots+e^{\theta_{N-1}}+1}\left[e^{\theta_{1}},\ldots,e^{\theta_{N-1}},1\right],$

and the value function $V_{\psi}(x)$ is represented in tabular form for
simplicity:

$V_{\psi}=[\psi_{1},\ldots,\psi_{N}]$

where $\psi\in\mathbb{R}^{N}$.

###  4.1. Loopy Markov chain

We test a toy example of the three-state loopy Markov chain which was
considered in [7, 39]. The transition probability matrix for the four state
$\\{0,1,2,3\\}$ is

$\tilde{K}=\left[\begin{array}[]{cccc}1&0&0&0\\\
\epsilon&\frac{1-\epsilon}{3}&\frac{1-\epsilon}{3}&\frac{1-\epsilon}{3}\\\
\epsilon&\frac{1-\epsilon}{3}&\frac{1-\epsilon}{3}&\frac{1-\epsilon}{3}\\\
\epsilon&\frac{1-\epsilon}{3}&\frac{1-\epsilon}{3}&\frac{1-\epsilon}{3}\end{array}\right],\quad\epsilon\in(0,1).$

The state 0 is the absorbing state $\partial$ and $\mathcal{E}=\\{1,2,3\\}$.
$K$ is the sub-matrix of $\tilde{K}$ corresponding to the states
$\\{1,2,3\\}$. With the probability $\epsilon$, the process exits
$\mathcal{E}$ directly from state 1, 2 or 3. The true quasi-stationary
distribution of this example is the uniform distribution for any $\epsilon$.

In order to show the advantage of our algorithm, we consider two cases: (1)
$\epsilon=0.1$ and (2) $\epsilon=0.9$. For a larger $\epsilon$, the original
Markov chain is very easy to exit so each iteration takes less time, but the
convergence rate of Vanilla algorithm is slower.

To quantify the accuracy of the learned quasi-stationary distribution, we
compute the $L_{2}$ norm of the error between the learned quasi-stationary
distribution and the true values.

(a) (b)
(c) (d)

Figure 1. The loopy Markov chain example with $\epsilon=0.1$. The figure shows
the log-log plots of $L_{2}$-norm error of Vanilla Algorithm (a), Projection
Algorithm(b), Polyak Averaging Algorithm (c) and our actor-critic algorithm
(d). The iteration for actor-critic algorithm is defined as one step of
gradient descent (“$t$” in Algorithm 1).

In Figure 1, we compute the QSD when $\epsilon=0.1$. We set the initial value
$\theta_{0}=[-1,1],\psi_{0}=[0,0,0],r_{0}=0$, the learning rate
$\eta_{n}^{\theta}=\max\\{1/{n^{0.1}},0.2\\},\eta_{n}^{\psi}=0.0001,\eta_{n}^{r}=0.0001$
and the batch size is 4. The step size for Projection Algorithm is
$\epsilon_{n}=n^{-0.99}$. Figure 2 is for the case when $\epsilon=0.9$ We set
the initial value $\theta_{0}=[4,-2],\psi_{0}=[0,0,0],r_{0}=0$, the learning
rate $\eta_{n}^{\theta}=0.04,\eta_{n}^{\psi}=0.0001,\eta_{n}^{r}=0.0001$ and
the batch size is 32. The step size for Projection Algorithm is
$\epsilon_{n}=n^{-0.99}$.

(a) (b)
(c) (d)

Figure 2. The loopy Markov chain example with $\epsilon=0.9$. The figure shows
the log-log plots of $L_{2}$-norm error of Vanilla Algorithm (a), Projection
Algorithm(b), Polyak Averaging Algorithm (c) and our actor-critic algorithm
(d).

###  4.2. M/M/1/N queue with finite capacity and absorption

Our second example is a M/M/1 queue with finite queue capacity. The 0 state
has been set as an absorbing state. The transition probability matrix on
$\\{0,\ldots,N\\}$ takes the form

$\tilde{K}=\left[\begin{array}[]{cccccccc}1&0&0&0&0&\ldots&0&0\\\
\mu_{1}&0&\lambda_{1}&0&0&\ldots&0&0\\\
0&\mu_{2}&0&\lambda_{2}&0&\ldots&0&0\\\
0&0&\mu_{3}&0&\lambda_{3}&\ldots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&&\vdots&\vdots\\\
0&0&0&0&0&&0&\lambda_{N-1}\\\ 0&0&0&0&0&\ldots&1&0\end{array}\right]$

where $\lambda_{i}=\frac{\rho_{i}}{\rho_{i}+1}$,
$\mu_{i}=\frac{1}{\rho_{i}+1}$, $i\in\\{1,2,\cdots,N-1\\}$. $\rho_{i}>1$ means
a higher chance to jump to right than to left. A larger $\rho_{i}$ will have
less probability of exiting $\mathcal{E}$. Note that $K(x,\mathcal{E})=1$ for
$x\in\\{2,\ldots,N\\}$. So $K_{\alpha}(x,y)=K(x,y)$ for any $\alpha$ if $x\neq
1$ and
$K_{\alpha}(1,y)=K(1,y)+\mu_{1}\alpha(y)=\begin{cases}\lambda_{1}+\mu_{1}\alpha(1)&y=1,\\\
\mu_{1}\alpha(y)&2\leq y\leq N\\\ \end{cases}.$ Then
${R}_{\theta}(x,y)=-\ln{\frac{K_{\alpha_{\theta}}(x,y)}{K_{\beta_{\theta}}(x,y)}}=0$
if $x\neq 1$ and by (20), the gradient is simplified as

$\nabla_{\theta}r(\theta)=\mathbf{E}_{Y}[\big{(}R_{\theta}(1,Y)-r(\theta)+V(Y)-V(1)\big{)}\nabla_{\theta}\ln
K_{\alpha_{\theta}}(1,Y)+\nabla_{\theta}\ln{K_{\beta_{\theta}}(1,Y)}]$

where $Y$ follows the distribution $K_{\alpha}(1,\cdot)$.

Figure 3. The QSD for M/M/1/500 queue with $\rho_{i}\equiv 1.25$ (left) and
$\rho_{i}=2-\frac{3}{2N-4}(i-1)$ (right).

We consider two cases: (1) a constant $\rho_{i}=1.25$ and (2) a state-
dependent $\rho_{i}=2-\frac{3}{2N-4}(i-1)$. Note $\rho_{i}=1$ gives an equal
probability of jumping to left and to right. So in case (1), there is a
boundary layer at the most right end and in case (2), we expect to see a peak
of the QSD near $i\approx 2N/3$. Figure 3 shows the true QSD in both cases. We
set $N=500$.

In Figure 4, we consider the case when $\rho_{i}=1.25$ and compute the $L_{2}$
errors. We set the initial value $\theta_{0}^{i}=-35+\frac{35}{498}(i-1)$ for
$i\in\\{1,2,\dots,498\\}$ and $\theta^{499}_{0}=3$, $\psi_{0}=[0,0,\dots,0]$,
$r_{0}=0$, and the learning rate
$\eta_{n}^{\theta}=0.0003,\eta_{n}^{\psi}=0.0001,\eta_{n}^{r}=0.0001$ and the
batch size is 64. The step size for Projection Algorithm is
$\epsilon_{n}=n^{-0.95}$. Figure 5 plots the errors for the state-dependent
$\rho_{i}=2-\frac{3}{2N-4}(i-1)$. We set the initial value
$\theta_{0}^{i}=8+\frac{35}{250}(i-1)$ for $i\in\\{1,2,\dots,250\\}$,
$\theta_{0}^{251}=44$, $\theta_{0}^{i}=43$ for $i\in\\{252,\dots,305\\}$,
$\theta_{0}^{306}=48$, $\theta_{0}^{307}=42$ and
$\theta_{0}^{i}=43-\frac{38}{293}(i-1)$ for $i\in\\{308,309,\dots,499\\}$,
$\psi_{0}=[0,0,\dots,0],r_{0}=0$ and the learning rate is
$\eta_{n}^{\theta}=0.0002$, $\eta_{n}^{\psi}=0.0001,\eta_{n}^{R}=0.0001$ with
batch size as 128. The step size for Projection Algorithm is
$\epsilon_{n}=n^{-0.95}$. Both figures demonstrate the actor-critic algorithm
performs quite well on this example.

In Table 1, we compared the CPU time of each algorithm in the M/M/1/500 queue
when they obtain the accuracy at $2\times 10^{-1}$. We found that our
algorithm cost less time on this example.

(a) (b)
(c) (d)

Figure 4. The M/M/1/500 queue with $\rho_{i}=1.25$. The figure shows the log-
log plots of $L_{2}$-norm error of Vanilla Algorithm (a), Projection
Algorithm(b), Polyak Averaging Algorithm (c) and our actor-critic algorithm
(d).

(a) (b)
(c) (d)

Figure 5. The M/M/1/500 queue with $\rho_{i}=2-\frac{3}{2N-4}(i-1)$. The figure shows the log-log plots of $L_{2}$-norm error of Vanilla Algorithm (a), Projection Algorithm(b), Polyak Averaging Algorithm (c) and our actor-critic algorithm (d).  Table 1. The CPU time of each algorithm in the M/M/1/500 queue when they obtain the accuracy at $2\times 10^{-1}$. Algorithm | Vanilla | Projection | Polyak Averaging | ac_$\alpha$
---|---|---|---|---
Time(s) | 1038.3279 | 429.6304 | 505.2299 | 186.9280
Time(s) | 753.9503 | 259.0671 | 268.5476 | 251.5370

##  5\. Summary and Conclusion

In this paper, we propose a reinforcement learning (RL) method for the quasi-
stationary distribution (QSD) in discrete time finite-state Markov chains. By
minimizing the KL-divergence of two Markovian path distributions induced by
the candidate distribution and the true target distribution, we introduce the
formulation in terms of RL and derive the corresponding policy gradient
theorem. We devise an actor-critic algorithm to learn the QSD in its
parameterized form $\alpha_{\theta}$. This formulation of RL can get benefit
from the development of the RL method and the optimization theory. We
illustrated our actor-critic methods on two numerical examples by using simple
tabular parametrization and gradient descent optimization. It has been
observed that the performance of our method is more prominent for large scale
problem

We only demonstrate the preliminary mechanism of the idea here, and there is
much space left for improving the efficiency and extensions in future works.
The generalization from the current consideration of finite-state Markov chain
to the jump Markov process and the diffusion case is in consideration. More
importantly, for very large or high dimensional state space, modern function
approximation methods like kernel methods or neural networks should be used
for the distribution $\alpha_{\theta}$ and the value function $V_{\psi}$. The
recent tremendous advancement of optimization techniques for policy gradient
in reinforcement learning could also contribute much to the efficiency
improvement of our current formulation.

## Acknowledgement

LL acknowledges the support of NSFC 11871486. XZ acknowledges the support of
Hong Kong RGC GRF 11305318.

## Appendix

In the appendix, we discuss the computation of the gradient of
$\nabla_{\theta}\ln K_{\alpha_{\theta}}$ and $\nabla_{\theta}\ln
K_{\beta_{\theta}}$. Note $\nabla_{\theta}\alpha_{\theta}$ is straightforward
since we model $\alpha$ in its parametrization form $\theta$. By definition
(4),

$\nabla_{\theta}\ln
K_{\alpha_{\theta}}(X_{t},X_{t+1})=\frac{1-K(X_{t},\mathcal{E})}{K_{\alpha_{\theta}}(X_{t},X_{t+1})}\nabla_{\theta}\alpha_{\theta}(X_{t+1})$

and

$\nabla_{\theta}\ln
K_{\beta_{\theta}}(X_{t},X_{t+1})=\frac{1-K(X_{t},\mathcal{E})}{K_{\beta_{\theta}}(X_{t},X_{t+1})}\nabla\beta_{\theta}(X_{t+1}).$

where $K_{\beta}(x,y)=K(x,y)+(1-K(x,\mathcal{E}))\beta(y)$. The vector
$K(x,\mathcal{E})$ for any $x$ can be pre-computed and saved in tabular form.

By (19), the one-step distribution $\beta$ is computed below

$\beta(X_{t+1})=\sum_{x}\alpha(x)\bigg{[}K(x,X_{t+1})+(1-K(x,\mathcal{E}))\alpha(X_{t+1})\bigg{]}\approx\frac{1}{n}\sum_{i=1}^{n}K(Z_{i},X_{t+1})+(1-K(Z_{i},\mathcal{E}))\alpha(X_{t+1})$

Here the samples $Z_{i}\sim\alpha$ and could be approximated by the stationary
distribution $\mu$ and thus one may simply use the known sample $X_{t}$ to
replace $Z_{i}$ with $n=1$.

To find $\nabla_{\theta}\beta_{\theta}$, we use stochastic approximation again

$\displaystyle\nabla_{\theta}\beta_{\theta}(X_{t+1})$
$\displaystyle=\sum_{x}\nabla\alpha_{\theta}(x)\left[K(x,X_{t+1})+(1-K(x,\mathcal{E})){\alpha_{\theta}}(X_{t+1})\right]+\left[\sum_{x}\alpha_{\theta}(x)(1-K(x,\mathcal{E}))\right]\nabla_{\theta}\alpha_{\theta}(y),$
$\displaystyle\approx\frac{1}{n}\sum_{i=1}^{n}\nabla\ln\alpha_{\theta}(Z_{i})\left[K(Z_{i},X_{t+1})+(1-K(Z_{i},\mathcal{E})){\alpha_{\theta}}(X_{t+1})\right]+(1-K(Z_{i},\mathcal{E}))\nabla_{\theta}\alpha_{\theta}(X_{t+1}).$

## References

  * [1] D. Aldous, B. Flannery, and J. L. Palacios, Two applications of urn processes the fringe analysis of search trees and the simulation of quasi-stationary distributions of markov chains, Probability in the engineering and informational sciences, 2 (1988), pp. 293–307.
  * [2] J. R. Artalejo, A. Economou, and M. J. Lopez-Herrero, Stochastic epidemic models with random environment: quasi-stationarity, extinction and final size, Journal of mathematical biology, 67 (2013), pp. 799–831.
  * [3] M. Bebbington, Parallel implementation of an aggregation/disaggregation method for evaluating quasi-stationary behavior in continuous-time markov chains, Parallel Computing, 23 (1997), pp. 1545–1559.
  * [4] M. Benaïm and B. Cloez, A stochastic approximation approach to quasi-stationary distributions on finite spaces, Electronic Communications in Probability, 20 (2015), pp. 1–13.
  * [5] N. Berglund and D. Landon, Mixed-mode oscillations and interspike interval statistics in the stochastic fitzhugh–nagumo model, Nonlinearity, 25 (2012), p. 2303.
  * [6] J. Blanchet, P. Glynn, and S. Zheng, Empirical analysis of a stochastic approximation approach for computing quasi-stationary distributions, in EVOLVE - A Bridge between Probability, Set Oriented Numerics, and Evolutionary Computation II, O. Schütze, C. A. Coello Coello, A.-A. Tantar, E. Tantar, P. Bouvry, P. Del Moral, and P. Legrand, eds., Berlin, Heidelberg, 2013, Springer Berlin Heidelberg, pp. 19–37.
  * [7] , Analysis of a stochastic approximation algorithm for computing quasi-stationary distributions, Advances in Applied Probability, 48 (2016), pp. 792–811.
  * [8] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe, Variational inference: A review for statisticians, Journal of the American Statistical Association, 112 (2017), pp. 859–877.
  * [9] S. Boyd, S. P. Boyd, and L. Vandenberghe, Convex optimization, Cambridge university press, 2004.
  * [10] F. Buckley and P. Pollett, Analytical methods for a stochastic mainland–island metapopulation model, Ecological modelling, 221 (2010), pp. 2526–2530.
  * [11] D. C. Chan, P. K. Pollett, and M. C. Weinstein, Quantitative risk stratification in markov chains with limiting conditional distributions, Medical Decision Making, 29 (2009), pp. 532–540.
  * [12] D. Clancy and S. T. Mendy, Approximating the quasi-stationary distribution of the sis model for endemic infection, Methodology and Computing in Applied Probability, 13 (2011), pp. 603–618.
  * [13] P. Collet, S. Martínez, and J. San Martín, Quasi-stationary distributions: Markov chains, diffusions and dynamical systems, Springer Science & Business Media, 2012.
  * [14] M. M. de Oliveira and R. Dickman, Quasi-stationary distributions for models of heterogeneous catalysis, Physica A: Statistical Mechanics and its Applications, 343 (2004), pp. 525–542.
  * [15] M. M. de Oliveira and R. Dickman, How to simulate the quasistationary state, Physical Review E, 71 (2005), p. 016129.
  * [16] G. Di Gesù, T. Lelièvre, D. Le Peutrec, and B. Nectoux, Jump markov models and transition state theory: the quasi-stationary distribution approach, Faraday discussions, 195 (2017), pp. 469–495.
  * [17] M. I. Dykman, T. Horita, and J. Ross, Statistical distribution and stochastic resonance in a periodically driven chemical system, The Journal of chemical physics, 103 (1995), pp. 966–972.
  * [18] M. I. Jordan, Z. Ghahramani, T. S. Jaakkola, and L. K. Saul, An Introduction to Variational Methods for Graphical Models, Machine Learning, 37 (1999), pp. 183–233.
  * [19] H. Kushner and G. G. Yin, Stochastic approximation and recursive algorithms and applications, vol. 35, Springer Science & Business Media, 2003\.
  * [20] A. Lambert, Population dynamics and random genealogies, Stochastic Models, 24 (2008), pp. 45–163.
  * [21] D. Landon, Perturbation et excitabilité dans des modeles stochastiques de transmission de l’influx nerveux, PhD thesis, Université d’Orléans, 2012.
  * [22] T. Lelièvre and F. Nier, Low temperature asymptotics for quasistationary distributions in a bounded domain, Analysis & PDE, 8 (2015), pp. 561–628.
  * [23] Q. Liu and D. Wang, Stein variational gradient descent: A general purpose bayesian inference algorithm, in Advances in Neural Information Processing Systems, D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, eds., vol. 29, Curran Associates, Inc., 2016.
  * [24] S. Martinez and J. San Martin, Quasi-stationary distributions for a brownian motion with drift and associated limit laws, Journal of applied probability, 31 (1994), pp. 911–920.
  * [25] S. Méléard, D. Villemonais, et al., Quasi-stationary distributions and population processes, Probability Surveys, 9 (2012), pp. 340–410.
  * [26] V. Mnih, K. Kavukcuoglu, D. Silver, A. Graves, I. Antonoglou, D. Wierstra, and M. Riedmiller, Playing atari with deep reinforcement learning, arXiv preprint arXiv:1312.5602, (2013).
  * [27] P. Pollett and D. Stewart, An efficient procedure for computing quasi-stationary distributions of markov chains by sparse transition structure, Advances in Applied Probability, 26 (1994), pp. 68–79.
  * [28] M. Pollock, P. Fearnhead, A. M. Johansen, and G. O. Roberts, The scalable langevin exact algorithm: Bayesian inference for big data, arXiv preprint arXiv:1609.03436, (2016).
  * [29] B. T. Polyak and A. B. Juditsky, Acceleration of stochastic approximation by averaging, SIAM journal on control and optimization, 30 (1992), pp. 838–855.
  * [30] M. Popova, O. Isayev, and A. Tropsha, Deep reinforcement learning for de novo drug design, Science advances, 4 (2018), p. eaap7885.
  * [31] D. Rezende and S. Mohamed, Variational inference with normalizing flows, in Proceedings of the 32nd International Conference on Machine Learning, F. Bach and D. Blei, eds., vol. 37 of Proceedings of Machine Learning Research, Lille, France, 07–09 Jul 2015, PMLR, pp. 1530–1538.
  * [32] D. C. Rose, J. F. Mair, and J. P. Garrahan, A reinforcement learning approach to rare trajectory sampling, New Journal of Physics, 23 (2021), p. 013013.
  * [33] A. Sani, D. Kroese, and P. Pollett, Stochastic models for the spread of hiv in a mobile heterosexual population, Mathematical biosciences, 208 (2007), pp. 98–124.
  * [34] D. Silver, T. Hubert, J. Schrittwieser, I. Antonoglou, M. Lai, A. Guez, M. Lanctot, L. Sifre, D. Kumaran, T. Graepel, et al., A general reinforcement learning algorithm that masters chess, shogi, and go through self-play, Science, 362 (2018), pp. 1140–1144.
  * [35] R. S. Sutton and A. G. Barto, Reinforcement learning: An introduction, MIT press, 2018.
  * [36] A. Q. Wang, G. O. Roberts, and D. Steinsaltz, An approximation scheme for quasi-stationary distributions of killed diffusions, Stochastic Processes and their Applications, 130 (2020), pp. 3193–3219.
  * [37] W. Wang and M. A. Carreira-Perpinán, Projection onto the probability simplex: An efficient algorithm with a simple proof, and an application, arXiv preprint arXiv:1309.1541, (2013).
  * [38] D. S. Watkins, Fundamentals of matrix computations, vol. 64, John Wiley & Sons, 2004.
  * [39] S. Zheng, Stochastic approximation algorithms in the estimation of quasi-stationary distribution of finite and general state space Markov chains, PhD thesis, Columbia University, 2014.