# Extreme first passage times for populations of identical rare events

James MacLaurin
Department of Mathematical Sciences
New Jersey Institute of Technology
Newark USA
<EMAIL_ADDRESS>
&Jay Newby
Department of Mathematical and Statistical Sciences
University of Alberta
Edmonton
<EMAIL_ADDRESS>

###### Abstract

A collection of identical and independent rare event first passage times is
considered. The problem of finding the fastest out of $N$ such events to occur
is called an extreme first passage time. The rare event times are singular and
limit to infinity as a positive parameter scaling the noise magnitude is
reduced to zero. In contrast, previous work has shown that the mean of the
fastest event time goes to zero in the limit of an infinite number of walkers.
The combined limit is studied. In particular, the mean time and the most
likely path taken by the fastest random walker are investigated. Using
techniques from large deviation theory, it is shown that there is a
distinguished limit where the mean time for the fastest walker can take any
positive value, depending on a single proportionality constant. Furthermore,
it is shown that the mean time and most likely path can be approximated using
the solution to a variational problem related to the single-walker rare event.

##  1 Introduction

Stochasticity is an ingredient used in a mathematical model best used when it
imbues a dynamical system with behavior that a simpler deterministic model
lacks. This is best seen when the connection between stochastic and
deterministic is explicit through a law of large number limit, in which the
noise, whatever its sources, vanishes in a so called weak-noise limit. For
simplicity, we will suppose that this limit is characterized by a positive
nondimensional parameter, $\epsilon>0$, where in the limit $\epsilon\to
0^{+}$, the stochastic process converges to a deterministic dynamical system.
For any behavior that exists in the stochastic model but not the
deterministic, this limit is singular. The random time over which an event
that characterizes this behavior occurs—commonly called a first passage time
or hitting time—diverges to infinity as $\epsilon\to 0^{+}$. Events of this
type are often called rare events. Their study is informed through exponential
asymptotic analysis, namely large deviation theory [6, 4].

There are many situations, particularly those involving exponential
proliferation, where it is not the mean first passage time of a single event
that is of interest, rather it is the first out of $N$ identical rare events
to occur that is relevant [32]. For example, the dynamics of the initial
outbreak in an epidemic are highly sensitive to when and where the first
infection occurs [36]. For another example, suppose that a tumor is comprised
of cells with a single common ancestor. In this case, the first-to-arrive
random walker is the first cell to accumulate enough mutations to cause
uncontrolled proliferation [9]. Another example is given by predator-prey
dynamics where the prey is sufficiently skilled at evading predators that
predation is a rare event [13]. In the scenario where a number of predators
are introduced to potential prey at the same time (e.g., a school of fish),
any given predator would successfully capture a prey a certain way, perhaps
skillfully. We might imagine a very different sequence of events for the first
predator to eat since this is an event that is doubly rare, in some sense.

Many examples can be found in the calcium signalling literature. The binding
of calcium to buffers in exocytosis can be characterized as a first-passage-
time problem for a large number of stochastic walkers [33, 23]. Similarly
calcium waves in the cell are often kicked off by the stochastic opening of a
single channel [30].

There are many examples of bacteria entering a dormant state to avoid the
harmful effects of an antibiotic, which kill during cell division. In some
cases, switching into and out of the dormant state is governed by a rare event
in a certain gene regulatory circuit [26, 25]. Once the antibiotics are no
longer present, switching out of this state allows the dormant subpopulation
to begin proliferating once again, driven by the first bacterium to make the
switch.

The problem can be stated succinctly as follows. Given $N$ independent and
identical random walkers (assumed to be continuous time Markov processes)
together with a first passage time event characterized by the random variable
$\tau>0$, what is the mean of $\tau^{N}_{\rm
first}=\min_{n}\\{\tau_{n},n=1,\ldots,N\\}$? Early work on the statistics of
extreme diffusion problems such as these was done by Varadhan [34], and it is
known that $\mathbb{E}[\tau^{N}_{\rm first}]$ scales like $C/\log(N)$ as
$N\to\infty$, where the constant $C$ depends on the statistics of the single
first passage time. This paper focusses on what happens in the dual limit
$N\to\infty$ and $\epsilon\to 0^{+}$: i.e. when the mean of the individual
first passage times blows up (i.e., $\mathbb{E}[\tau]\to\infty$).

Holcman et al recently proposed extreme first passage times as a useful
characterization of biological processes, particularly at the level of
cellular physiology, for example, synaptic transmission [23], and the
activation of genes by transcription factors [31]. Lawley et al [14, 15, 17]
have recently studied the extreme first-passage-time statistics of a large
number of random walkers in a range of scenarios. Most of these studies
concern situations where $N$ is the dominant asymptotic parameter. That is, it
is assumed that the magnitude of the noise experienced by an individual
walkers is $O(1)$, and it is determined how the optimal trajectory and first-
hitting-time is changes as $N\to\infty$, often using the theory of extreme
value statistics [1]. A problem with this approach is that it implies that the
first-hitting-particle almost instantaneously hits the target. It goes there
so quickly that the effect of advection is negligible. Indeed some scholars
have argued that large jumps of Brownian motions such as this are unphysical
[12].

In fact, Lawley and Madrid [19] have extended these works to also consider the
effect of small noise. Relevant to our paper, they argue that the relative
benefit of increased numbers of walkers is relatively weak (scaling as
$1/\log(N)$), and so it unclear if many biological systems are in this regime.
We also note the work of Weiss [37], who determine first-hitting-time
estimates when the walkers are initially uniformly distributed throughout the
domain.

Let us assume that the single walker rare event starts at a deterministic
attractor and terminates at some set that is well separated from the
deterministic attractor but contained within its basin of attraction. One
simple example is a 1D OU process, $dX=-\mu X+\sqrt{2\epsilon}dW$, where the
attractor is at $x=0$ and the target is at some $x=L>0$. This example is
sometimes referred to as escape from a potential well.

One way to describe a rare event is through a so called maximum likelihood
trajectory (MLT) [6, 18, 20, 22, 27]. It is often true that a single state
space trajectory describes a rare event. That is, if one takes a stochastic
system at the terminus of a rare event and follows the process back in time,
examining the sequence of events that preceded the final step of the rare
event, it is often found that a single trajectory describes this sequence of
events. To make this more concrete, we can imagine a distribution of possible
trajectories and define the MLT as the mode, or peak, of this distribution. It
is possible to formulate asymptotic approximations of MLTs as solutions to a
certain variational problem [6, 4].

The single walker rare event can be described in two parts: a lag phase and a
transition phase [20]. The rare event begins with an exponentially large (in
$\epsilon^{-1}$) lag phase where for the great majority of time the process
fluctuates within a $\sqrt{\epsilon}$ domain around the attractor. During the
lag phase, there will be brief and rare larger excursions further away from
the attractor (but not large enough to attain the target $L$) that are quickly
pulled back to the attractor. After the lag phase, there is a transitional
phase characterized by a MLT that starts on the boundary of the
$\sqrt{\epsilon}$ neighborhood and terminates at the target. The total time of
the transitional phase scales like $\log(1/\sqrt{\epsilon})$ as $\epsilon\to
0$. Although the mean time of both phases is infinite in the limit
$\epsilon\to 0$, the first phase is comparatively much longer.

In this paper we answer the following question. If we understand the MLT for a
single rare event, does this trajectory also describe the minimum of $N$
identical rare events? Perhaps the extreme first passage time is an unlikely
deviation from the typical rare event, an event that is doubly rare, in some
sense. It is known that (in the limit $N\to\infty$ with $\epsilon$ fixed) the
MLT for an extreme event is a geodesic, relative to the geometry induced by
the diffusion coefficients [34]. For example, simple additive white noise
(constant diffusivity) results in a MLT that is a direct linear path from the
starting point to the target. This effectively means that there are so many
walkers that the first to arrive is perturbed by Brownian fluctuations that
send it straight to the target in a short period of time, and therefore the
effect of a "drift velocity" term on the MLT is negligible once there are
enough walkers, even if this velocity pushes the process away from the target.

We are left with the general possibility that for $N$ finite the first-to-
arrive walker follows the single walker MLT in the limit $\epsilon\to 0$. On
the other hand, for $\epsilon$ fixed, the first-to-arrive walker follows a
straight line path which will, in general, be very different than the single
walker MLT. Is it possible to describe the extreme first passage time MLT in
the dual limit $\epsilon\to 0$ and $N\to\infty$? Does the trajectory smoothly
transform from the single walker MLT to the straight line path for
sufficiently large $N$?

Suppose we have a rare event such that
$\tau_{n}\sim\text{Exponential}(\lambda_{\epsilon}),$ asymptotically as
$\epsilon\to 0^{+}$. Previous work on this problem identified two regimes
[19]. For fixed $N>0$, $\mathbb{E}[\tau^{N}_{\rm
first}]\propto\frac{1}{\lambda_{\epsilon}N}$ for $\epsilon\to 0^{+}$. For
fixed $\epsilon>0$, We know that $\mathbb{E}[\tau^{N}_{\rm first}]$ scales
like $1/\log(N)$ as $N\to\infty$. We are interested in the crossover between
the $1/N$ regime and the $1/\log(N)$ regime as $N$ grows.

In this paper, we will show that if we combine these two perspectives, a
single picture emerges wherein there are three regimes for
$\mathbb{E}[\tau^{N}_{\rm first}]$ in the double limit $\epsilon\to 0$ and
$N\to\infty$. For smaller values of $N$ we have the first regime where the
first-to-arrive walker follows the single walker MLT. In this case,
$\mathbb{E}[\tau^{N}_{\rm first}]$ is decreased like $1/N$ by eliminating time
spent in the exponentially long lag phase where the walker fluctuates near the
attractor. As $N$ is increased sufficiently to eliminate the lag phase, we
enter the second regime. Reducing $\mathbb{E}[\tau^{N}_{\rm first}]$ now
requires shortening the duration of the transition phase of the rare event.
The extreme MLT smoothly deforms away from the single walker MLT in order to
speed up the transition time. Eventually, for $N$ sufficiently large, it is no
longer possible to improve $\mathbb{E}[\tau^{N}_{\rm first}]$ by deforming the
extreme MLT because it is a straight line from the attractor to the target. In
this regime, $\mathbb{E}[\tau^{N}_{\rm first}]$ is reduced like $1/\log(N)$ by
a transition path that moves faster and faster along the direct line extreme
MLT.

To our knowledge, the second regime described above was previously unknown.
Moreover, we find that in the second regime, it is always possible take the
double limit $\epsilon\to 0$ and $N\to\infty$ in such a way that
$\mathbb{E}[\tau^{N}_{\rm first}]\to\tau_{\infty}$, where $\tau_{\infty}$ can
take any positive value.

The remainder of the paper is organized as follows. First, we include some
necessary background from large deviation theory in Section 2. We then state
our main result in Section 3. Finally, we illustrate our main result in
Section 4 with a simple example using an Ornstein-Uhlenbeck process where many
exact results are obtainable. We conclude by discussing phenomena that are
known to exist in general nonlinear systems along with numerical strategies
for computing MLTs.

##  2 Background on Large Deviations Theory

Consider the diffusion of $N$ independent walkers in $\mathbb{R}^{d}$, each
satisfying the stochastic differential equation

$dX^{j}_{t}=f(X^{j}_{t})dt+\sqrt{2\epsilon}dW^{j}_{t},$ (1)

with $1\leq j\leq N$. Here $W^{j}_{t}:=(W^{1,j}_{t},\ldots,W^{d,j}_{t})$ is
$d$-dimensional Brownian Motion. It is assumed that $0$ is a globally-
attracting fixed point of the dynamics (when $\epsilon=0$). The initial
condition is assumed to be zero for every walker.

The main goal of this paper is to determine the probabilistic distribution of
the first-hitting-time of some target $L\in\mathbb{R}^{d}$. In particular, we
wish to understand how the expected first hitting time depends on $\epsilon$
(which is asymptotically small) and $N$ (which is asymptotically large).
Concurrently, we wish to understand the most likely trajectory following by
the first walker to hit the target.

Define a small ball about $L$ to be

$B_{\delta}(L)=\big{\\{}z\in\mathbb{R}^{d}:\|z-L\|\leq\delta\big{\\}},$ (2)

for some $\delta=\delta(\epsilon,N)$ that converges to zero very slowly as
$\epsilon\to 0$ and $N\to\infty$. As long as $\delta$ converges to 0
relatively slow, the leading order asymptotic estimate of the first hitting-
time will only depend on $\epsilon$ and $N$. Let $\tau^{N}_{\rm first}$ be the
first hitting time of $L$ over all $N$ walkers, and let $\tau_{j}$ be the
first hitting time of a single walker. Then, we have

$\displaystyle\tau_{j}$ $\displaystyle=\inf\big{\\{}t\geq 0:X^{j}_{t}\in
B_{\delta}(L)\big{\\}}$ (3) $\displaystyle\begin{split}\tau^{N}_{\rm
first}&=\inf\left\\{\tau_{1},\ldots,t_{N}\right\\}\end{split}$ (4)

Classically, the asymptotics of the first-hitting-time of a single walker has
been studied using Freidlin-Wentzell Theory (also referred to as the Theory of
Large Deviations [4]). We now briefly recap the relevant aspects of this
theory (mostly following the treatment in [4]).

First, for some time $T>0$, we must define the rate function $I_{T}$ that
gives the asymptotic probability of different trajectories. $I_{T}$ is a map
from $\mathcal{C}([0,T],\mathbb{R}^{d})$ to $\mathbb{R}^{+}$. It is stipulated
that $I_{T}(z)<\infty$ only if $z\in H_{1}$ (the space of all continuous
functions which possess a bounded derivative for Lebesgue-almost-every point)
to $\mathbb{R}^{+}$. In this case, the rate function is defined to be

$\displaystyle I_{T}(y)$
$\displaystyle=\int_{0}^{T}L(y(t),\dot{y}(t))dt\;\;\;\text{ where }$ (5)
$\displaystyle L(y,\dot{y})$
$\displaystyle=\frac{1}{4}\left(\dot{y}-f(y)\right)^{T}\left(\dot{y}-f(y)\right).$
(6)

$L(y,\dot{y})$ is called the Lagrangian. Formally, the Large Deviations
Principle expresses the asymptotics of certain sets in the following manner.
Suppose that
$\mathcal{A},\mathcal{O}\subset\mathcal{C}\big{(}[0,T],\mathbb{R}^{d}\big{)}$
are (respectively) open and closed subsets of the space of all continuous
functions, whose value at time $0$ is $0$. The topology on
$\mathcal{C}\big{(}[0,T],\mathbb{R}^{d}\big{)}$ is the topology of uniform
convergence - this means that every open set must be able to be written as a
union (potentially infinite) of small balls of the form

$\displaystyle\bigg{\\{}z\in\mathcal{C}\big{(}[0,T],\mathbb{R}^{d}\big{)}\;:\;\sup_{t\in[0,T]}\sup_{1\leq
j\leq d}\big{|}z^{j}(t)-y^{j}(t)\big{|}<\xi\bigg{\\}},$ (7)

for some $\xi>0$ and $y\in\mathcal{C}\big{(}[0,T],\mathbb{R}^{d}\big{)}$. The
Large Deviations results originally due to Freidlin and Wentzell state that
[4, 6]

$\displaystyle\lim_{\xi\to
0^{+}}\sup_{\epsilon\leq\xi}\epsilon\log\mathbb{P}\big{(}x^{j}_{[0,T]}\in\mathcal{A}\big{)}\leq-\inf_{z\in\mathcal{A}}I_{T}(z),$
(8) $\displaystyle\lim_{\xi\to
0^{+}}\inf_{\epsilon\leq\xi}\epsilon\log\mathbb{P}\big{(}x^{j}_{[0,T]}\in\mathcal{O}\big{)}\geq-\inf_{z\in\mathcal{A}}I_{T}(z).$
(9)

###  2.1 Large Deviation Theory of maximum likelihood trajectories

Next, we outline how the Large Deviations results can be used to estimate MLTs
(for a single walker, in the small noise limit). The theory in this sub-
section is already well-established (see [8] for a thorough review). Define
the corresponding Hamiltonian
$\mathcal{H}:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$,

$\displaystyle\mathcal{H}(x,p)=\sup_{z\in\mathbb{R}^{d}}\left\\{z^{T}p-L(x,z))\right\\}=p^{T}f(x)+p^{T}p.$
(10)

Because the Lagrangian is convex in its second argument, Legendre duality
implies that

$\displaystyle
L(y,z)=\sup_{p\in\mathbb{R}^{d}}\left\\{z^{T}p-\mathcal{H}(x,p)\right\\},$
(11)

which is consistent with (6). Suppose that we wish to estimate (i) the
probability of $x^{j}_{t}=z$ for some $z\in\mathbb{R}^{d}$, and (ii) the most
likely trajectory followed by the walker in attaining this value. In other
words, we wish to minimize the cost function up to time $t$, for
$z\in\mathbb{R}^{d}$

$\displaystyle V(z,t)=\inf_{y\in H_{1}:y(t)=z}I_{t}(y),$ (12)

and find the trajectory $\varphi\in H_{1}$ that attains this minimum. This
means that $V(z,t)$ is such that, for $\delta\ll 1$,

$\displaystyle\mathbb{P}\big{(}\|X_{t}-z\|\leq\delta\big{)}\simeq\exp\big{(}-\epsilon^{-1}V(z,t)\big{)}$
(13)

It is well known (see Chapter 5 of Dembo and Zeitouni [4]) that the cost
function satisfies the following Hamilton-Jacobi equation,

$\frac{\partial}{\partial t}V(x,t)=\mathcal{H}(x,\nabla_{x}V).$ (14)

Finally, taking $t\to\infty$ we obtain the so-called quasi-potential
$V:\mathbb{R}^{d}\to\mathbb{R}^{+}$,

$V(z)=\inf_{t\geq 0}\inf_{\stackrel{{\scriptstyle y\in
H_{1}}}{{y(t)=z}}}I_{t}(z)=\lim_{t\to\infty}V(z,t).$ (15)

The quasipotential is a stationary solution to the Hamilton-Jacobi equation
(14). For $\Lambda>0$, we define $\Psi_{x}(\Lambda)$ to be the shortest time
that it takes for the cost-function to hit $\Lambda$,

$\displaystyle\Psi_{x}(\Lambda)=\inf\big{\\{}t\geq
0\;:\;V(x,t)=\Lambda\big{\\}}.$ (16)

By Hamilton’s Principle, the path $\varphi_{t|T}$ that minimizes the rate
functional (5) satisfies

$\begin{split}\dot{\varphi}_{t|T}&=\nabla_{p}\mathcal{H}(\varphi_{t|T},p_{t}),\\\
\dot{p}_{t|T}&=-\nabla_{x}\mathcal{H}(\varphi_{t|T},p_{t|T}),\end{split}$ (17)

along with the boundary conditions,

$\varphi_{0|T}=0,\;\;\varphi_{T|T}=L.$ (18)

An important property of the ODEs (17) is that the they admit a constant of
motion. In particular, we have that
$\mathcal{H}(\varphi_{t|T},p_{t|T})\equiv\text{const}$ over any solution of
(17). See [8] for a review of numerical methods one may employ to compute
solutions of (17).

In addition to satisfying the Hamilton-Jacobi equation (14), the cost function
(12) can be integrated along the optimal paths satisfying (17). We define the
cost function along an MLT with

$S_{t|T}=V(\varphi_{t|T},t).$ (19)

Differentiating the cost function with respect to time along an optimal path
yields

$\begin{split}\dot{S}=\frac{d}{dt}V(\varphi_{t|T},t)&=(\nabla_{x}V(\varphi_{t|T},t))^{T}\frac{\partial\varphi_{t|T}}{\partial
t}+\frac{\partial V}{\partial t}(\varphi_{t|T},t)\\\
&=p_{t|T}^{T}\nabla_{p}\mathcal{H}(\varphi_{t|T},p_{t|T})-\mathcal{H}(\varphi_{t|T},p_{t|T})\\\
&=p_{t|T}^{T}p_{t|T}.\end{split}$ (20)

Hence, we can supplement the system (17) with the auxiliary equation,

$\dot{S}_{t|T}=p_{t|T}^{T}p_{t|T},\quad S_{0|T}=0.$ (21)

Notice that it immediately follows from the above equation that
$\frac{d}{dt}V\geq 0$ along optimal paths. The quasi-potential (15) is given
by the "zero energy" solutions to (17). In particular, the
$\mathcal{H}(\varphi_{t|T},p_{t|T})=0$ solutions to (17) are consistent with
the steady state solution to the Hamilton-Jacobi equation (14).

There are several more important properties of the Hamiltonian dynamical
system (17) that are relevant in our current setting. For $x$ in a
neighborhood of a stable attractor of the deterministic dynamics
($\dot{\varphi}=f(\varphi))$, the quasi-potential is determined by optimal
paths that converge backward in time to the stable attractor as
$t-t_{0}\to\infty$. To see this we first notice that
$\nabla_{p}\mathcal{H}(x,0)=f(x)$, which means that the deterministic flows
exist on a manifold of the Hamiltonian dynamical system (17). Since,
$\nabla_{x}\mathcal{H}(x,0)\equiv 0$ for all $x$, it follows that
deterministic fixed points, $x_{*}$, satisfying $f(x_{*})=0$, are also fixed
points of the Hamiltonian dynamical system. (Note that it is a consequence of
Liouville’s Theorem that there are no stable or unstable attractors in the
Hamiltonian dynamical system, only saddles and centers.) Finally, since
$\mathcal{H}(x_{*},0)=0$, the Hamiltonian function must be zero along all
trajectories converging to or diverging from $(x,p)=(x_{*},0)$.

Together, these properties lead to the following understanding of the zero-
energy optimal paths. If the starting point is the attractor $x_{*}$ and the
end point is $L$ (within the basin of attraction) then $\dot{S}_{t|T}\geq 0$
and the optimal paths describe MLTs. If, on the other hand, the starting point
$L$ and the end point is $x_{*}$ then $\dot{S}_{t|T}\equiv 0$, $p_{t|T}\equiv
0$, and $\varphi_{t|T}$ is a solution to the deterministic dynamics.

##  3 Main result

In this section, we consider the first passage time problem (4). We will show
that there are three different regimes for $\tau^{N}_{\rm first}$, depending
on the relative magnitude of $\epsilon$ and $N$ as they approach their limit.
It is the asymptotic value of

$\Lambda=\epsilon\log(N)$ (22)

that determines which regime is entered as $\epsilon\to 0^{+}$ and
$N\to\infty$. It is assumed throughout this paper that the limit
$\lim_{\epsilon\to 0^{+}}\epsilon\log(N)$ exists (and could be infinite).

  * •

Regime 1: Labor shortage. $\Lambda\leq V(L)$ so that $N\leq
N_{*}=\exp(V(L)\epsilon^{-1})$. In this regime the trajectory followed by the
first-arriving walker is with high probability approximately the same as in
the single-walker case (i.e. when $N=1$). The expected time
$\mathbb{E}[\tau^{N}_{\rm first}]$ that it takes the first walker to arrive is
the expected time for the single-walker case, divided by $N$. Furthermore
$\mathbb{E}[\tau^{N}_{\rm first}]$ diverges to $\infty$ as $\epsilon\to 0^{+}$
and $N\to\infty$.

  * •

Regime 2. Balanced labor. $V(L)<\Lambda<\infty$ so that
$N>N_{*}=\exp(V(L)\epsilon^{-1})$. In this regime, the trajectory followed by
the first-arriving walker is with high probability different from the
trajectory in the single-walker case. Furthermore the expected time that it
takes the first walker to arrive does not diverge as $N\to\infty$ and
$\epsilon\to 0^{+}$. Indeed, the expected time $\mathbb{E}[\tau^{N}_{\rm
first}]$ that it takes the first walker to arrive is $\mathbb{E}[\tau^{N}_{\rm
first}]=\Psi_{L}(\Lambda)$, where $\Psi_{L}(\Lambda)$ is defined in (16).

  * •

Regime 3. Labor surplus. $\Lambda\to\infty$. In this regime, there are so many
walkers that the first-to-arrive walker almost instantaneously travels
straight to the target. The expected first hitting-time goes to $0$. The first
walker moves to the target so quickly that the advection term $f$ has a
negligible effect on its trajectory.

###  3.1 Proof in the case of a diverging first-hitting-time

Lets first consider the case that $\Lambda<V(L)$, and we recall that the
number of walkers $N$ is scaled according to (22). The main result of this
section is that for any $\zeta>0$,

$\displaystyle\lim_{\delta\to 0^{+}}\underset{\epsilon\to
0}{\overline{\lim}}\epsilon\log\mathbb{P}\big{(}|\epsilon\log(\tau^{N}_{\rm
first})-\Lambda+V(L)|>\zeta\big{)}<0.$ (23)

Once (23) has been demonstrated, it follows immediately (thanks to the Borel
Cantelli Lemma) that with probability one,

$\epsilon\log(\tau^{N}_{\rm first})\to V(L)-\Lambda$ (24)

as one (first) asymptotes $\epsilon\to 0^{+}$ and then (second) asymptotes
$\delta\to 0$.

To show (23), it suffices to demonstrate the following identities

$\displaystyle\lim_{\delta\to 0^{+}}\underset{\epsilon\to
0}{\overline{\lim}}\epsilon\log\mathbb{P}\big{(}\epsilon\log(\tau^{N}_{\rm
first})-\Lambda+V(L)>\zeta\big{)}$ $\displaystyle<0,$ (25)
$\displaystyle\lim_{\delta\to 0^{+}}\underset{\epsilon\to
0}{\overline{\lim}}\epsilon\log\mathbb{P}\big{(}\epsilon\log(\tau^{N}_{\rm
first})-\Lambda+V(L)<-\zeta\big{)}$ $\displaystyle<0.$ (26)

First, we recap some known results concerning the asymptotics of the first
hitting time of a single walker (see Freidlin and Wentzell [6] for the
original proof, and also [4, Section 5.7]). For any $0\leq b\leq V(L)$ and any
$a>0$ (and typically $a\ll 1$),

$\displaystyle\lim_{\delta\to 0^{+}}\lim_{\zeta\to 0^{+}}\lim_{\epsilon\to
0}\left|\epsilon\log\mathbb{P}\big{(}\epsilon\log(\tau_{1})\in[b-\zeta,b+\zeta]\big{)}-b+V(L)\right|$
$\displaystyle=0$ (27) $\displaystyle\lim_{\delta\to 0}\lim_{\epsilon\to
0}\epsilon\log\mathbb{P}\big{(}\epsilon\log(\tau_{1})<b-a\big{)}$
$\displaystyle<b-V(L)$ (28) $\displaystyle\lim_{\delta\to 0}\lim_{\epsilon\to
0}\epsilon\log\mathbb{P}\big{(}\epsilon\log(\tau_{1})>b+a\big{)}$
$\displaystyle<b-V(L).$ (29)

Next we consider multiple walkers, exploiting the fact that they are
identically distributed. To prove (25), using a union of events bound,

$\displaystyle\mathbb{P}\big{(}\epsilon\log(\tau^{N}_{\rm
first})<V(L)-\Lambda-\zeta\big{)}$
$\displaystyle\leq\sum_{j=1}^{N}\mathbb{P}\big{(}\epsilon\log(\tau_{j})<V(L)-\Lambda-\zeta\big{)}$
(30)
$\displaystyle=N\mathbb{P}\big{(}\epsilon\log(\tau_{1})<V(L)-\Lambda-\zeta\big{)},$
(31)

since the walkers are identically distributed. It thus follows from (28) and
(31) that

$\underset{\epsilon\to
0}{\overline{\lim}}\epsilon\log\mathbb{P}\big{(}\epsilon\log(\tau^{N}_{\rm
first})<V(L)-\Lambda-\zeta\big{)}<0,$ (32)

and therefore (25) holds. Next, to prove (26), using the fact that the walkers
are independent and identically distributed, for any $\zeta>0$,

$\displaystyle\mathbb{P}\big{(}\epsilon\log(\tau^{N}_{\rm
first})>V(L)-\Lambda+\zeta\big{)}$
$\displaystyle=\prod_{j=1}^{N}\big{\\{}1-\mathbb{P}\big{(}\epsilon\log(\tau_{j})\leq
V(L)-\Lambda+\zeta\big{)}\big{\\}}$ (33)
$\displaystyle\leq\big{\\{}1-\exp\big{(}-\epsilon^{-1}(-\Lambda-\xi)\big{)}\big{\\}}^{N},$
(34)

for some $\xi>0$, as long as $\delta$ is sufficiently small, thanks to the
bound in (29). Taking logarithms of both sides, before multiplying both sides
by $\epsilon$, we thus find that (making use of the inequality $\log(1-u)\leq
u$ for any $u>0$ in the second line)

$\displaystyle\underset{\epsilon\to
0}{\overline{\lim}}\epsilon\log\mathbb{P}\big{(}\epsilon\log(\tau^{N}_{\rm
first})>V(L)-\Lambda+\zeta\big{)}$ $\displaystyle\leq\underset{\epsilon\to
0}{\overline{\lim}}\epsilon
N\log\big{\\{}1-\exp\big{(}-\epsilon^{-1}(\Lambda-\xi)\big{)}\big{\\}}$ (35)
$\displaystyle\leq-\underset{\epsilon\to 0}{\overline{\lim}}\epsilon
N\exp\big{(}-\epsilon^{-1}(\Lambda-\xi)\big{)}$ (36)
$\displaystyle=-\underset{\epsilon\to
0}{\overline{\lim}}\epsilon\exp\big{(}\epsilon^{-1}\xi)\big{)}$ (37)
$\displaystyle=-\infty,$ (38)

making use of the scaling of $N$ in (22).

###  3.2 Proof in the case of a non-diverging first-hitting-time

We now consider the case that $\Lambda>V(L)$. In this case, the average first-
hitting-time does not diverge as $\epsilon\to 0^{+}$. Let
$T_{\Lambda}=\Psi_{L}(\Lambda)$ (as defined in (16)), which means that

$\displaystyle\Lambda=V(L,T_{\Lambda}).$ (39)

Indeed, it must be the case that if $0<S<T$, then necessarily $V(L,S)\geq
V(L,T)$.

We are going to show that for any $\zeta>0$,

$\displaystyle\lim_{\delta\to 0}\underset{\epsilon\to
0}{\overline{\lim}}\epsilon\log\mathbb{P}\big{(}|\tau^{N}_{\rm
first}-T_{\Lambda}|>\zeta\big{)}<0.$ (40)

Applying the Borel-Cantelli Lemma to (40) one obtains that

$\tau^{N}_{\rm first}\to T_{\Lambda}$ (41)

as (firstly) (i) $\epsilon\to 0^{+}$, and then (secondly) (ii) $\delta\to
0^{+}$.

Large Deviations estimates imply that for any $a>0$ (and typically $a\ll 1$),

$\displaystyle\lim_{\delta\to 0}\lim_{\zeta\to 0}\lim_{\epsilon\to
0}\bigg{|}\epsilon\log\mathbb{P}\big{(}\tau_{1}\in[T_{\Lambda}-\zeta,T_{\Lambda}+\zeta]\big{)}+\Lambda\bigg{|}$
$\displaystyle=0$ (42) $\displaystyle\lim_{\delta\to 0}\lim_{\epsilon\to
0}\epsilon\log\mathbb{P}\big{(}\tau_{1}<T_{\Lambda}-a\big{)}$
$\displaystyle<-\Lambda$ (43) $\displaystyle\lim_{\delta\to
0}\lim_{\epsilon\to
0}\epsilon\log\mathbb{P}\big{(}\tau_{1}>T_{\Lambda}+a\big{)}$
$\displaystyle<-\Lambda$ (44)

We are going to demonstrate that for any $\zeta>0$ (typically $\zeta\ll 1$),

$\displaystyle\lim_{\delta\to 0}\underset{\epsilon\to
0}{\overline{\lim}}\epsilon\log\mathbb{P}\big{(}\tau^{N}_{\rm
first}-T_{\Lambda}>\zeta\big{)}<0,$ (45) $\displaystyle\lim_{\delta\to
0}\underset{\epsilon\to
0}{\overline{\lim}}\epsilon\log\mathbb{P}\big{(}\tau^{N}_{\rm
first}-T_{\Lambda}<-\zeta\big{)}<0.$ (46)

To prove (46), by a union-of-events bound,

$\displaystyle\mathbb{P}\big{(}\tau^{N}_{\rm first}-T_{\Lambda}<-\zeta\big{)}$
$\displaystyle\leq\sum_{j=1}^{N}\mathbb{P}\big{(}\tilde{\tau}^{j}-T_{\Lambda}<-\zeta\big{)}$
(47) $\displaystyle\leq N\mathbb{P}\big{(}\tau^{1}-T_{\Lambda}<-\zeta\big{)}$
(48) $\displaystyle\leq
N\exp\big{(}-\epsilon^{-1}\Lambda-\epsilon^{-1}\delta_{\zeta}\big{)}$ (49)

for some $\delta_{\zeta}>0$, thanks to (43). Since
$N\simeq\exp(\epsilon^{-1}\Lambda)$, we may conclude (46). For (45), since the
walkers are independent,

$\displaystyle\mathbb{P}\big{(}\tau^{N}_{\rm first}-T_{\Lambda}>\zeta\big{)}=$
$\displaystyle\prod_{j=1}^{N}\mathbb{P}\big{(}\tau^{j}-T_{\Lambda}>\zeta\big{)}$
(50) $\displaystyle\leq$
$\displaystyle\big{(}\exp(-\epsilon^{-1}\Lambda-\epsilon^{-1}\delta_{\zeta})\big{)}^{N},$
(51)

for some $\delta_{\zeta}>0$, thanks to (44). Since $N\to\infty$ as
$\epsilon\to 0$, it is immediate that

$\displaystyle\lim_{\epsilon\to
0^{+}}\epsilon\log\mathbb{P}\big{(}\tau^{N}_{\rm
first}-T_{\Lambda}>\zeta\big{)}=-\infty.$ (52)

To finish this section, we must also demonstrate that the most likely path
follows that implied by the Hamilton’s equations (17). To this end, let
$\phi_{T,z}\in\mathcal{C}^{2}([0,T],\mathbb{R}^{d})$ be the optimal path up to
time $T$. That is, $\phi_{T,z}$ is such that $\phi_{T,z}(0)=0$,
$\phi_{T,z}(T)=z$, and it satisfies the Euler-Lagrange equations (17) for all
times in $[0,T]$. Next, we note that the solution of the Euler-Lagrange
equations is indeed optimal (in the single-walker case). For any $\upsilon>0$,
it holds that

$\displaystyle\lim_{\epsilon\to
0}\epsilon\log\mathbb{P}\bigg{(}x^{j}_{T}=z,\sup_{t\in[0,T]}\|X^{j}_{t}-\phi_{T,z}(t)\|\geq\upsilon\bigg{)}<\lim_{\epsilon\to
0}\epsilon\log\mathbb{P}\bigg{(}X^{j}_{T}=z\bigg{)}.$ (53)

##  4 Example: Ornstein-Uhlenbeck process

Consider a 2D Ornstein-Uhlenbeck (OU) process given by

$dX_{t}=AX_{t}dt+\sqrt{2\epsilon}dW,\quad X_{0}=0,$ (54)

where

$A=\begin{bmatrix}-\mu&-\omega\\\ \omega&-\mu\end{bmatrix},$

and $\mu>0$. With this choice of matrix $A$, the dynamics are rotational; that
is, the origin is a stable spiral for the deterministic dynamics. We consider
the single walker rare event

$\tau=\inf\\{t>0\mid\|X_{t}-X_{0}\|=\|L\|\\},$ (55)

for $L\neq 0$. Because the probability law is invariant under a rotation of
the axes, the large deviations rate function and MLT are also invariant under
a rotation of the axes. In other words, any MLT for any point $\tilde{L}$
(where $\|\tilde{L}\|=\|L\|$, and $\tilde{L}$ can be obtained by rotating $L$
by some angle $\alpha$) up to some time $t$, can be obtained by taking the MLT
for $L$ up to time $t$, and rotating it uniformly by $\alpha$.

The linearity of the OU dynamics facilitates analytic formulae. Indeed, given
our initial condition, the process can be described as a multivariate normal
random variable for any countable set of times. The eigenvalues of $A$ are
$-\mu\pm i\omega$.

The solution to (54) can be written as

$X_{t}=\sqrt{2\epsilon}\int_{0}^{t}e^{-\mu(t-s)}R_{t-s}dW_{s}.$ (56)

where

$\displaystyle R_{t}$ $\displaystyle=\begin{bmatrix}\cos(\omega
t)&-\sin(\omega t)\\\ \sin(\omega t)&\cos(\omega t)\end{bmatrix}.$ (57)

The Hamiltonian defined by (10) is given by

$\mathcal{H}(x,p)=p^{T}Ax+\|p\|^{2}$ (58)

For our problem (17) and (21) become

$\begin{split}\dot{\varphi}_{t|T}&=A\varphi_{t|T}+2p_{t|T},\\\
\dot{p}_{t|T}&=-A^{T}p_{t|T},\\\ \dot{S}_{t|T}&=\|p_{t|T}\|^{2},\end{split}$
(59)

with

$\varphi_{0|T}=0,\,\varphi_{T|T}=L,\,S_{0|T}=0.$ (60)

The initial condition $p_{0|T}$ determines the constant of motion $H$, which
in turn, determines the final time $T$. In general, suboptimal MLTs can be
computed by increasing $H$ from zero. The quantity $p_{t|T}$ is sometimes
referred to as the conjugate momentum, and the initial value $p_{0|T}$ can be
thought of as some initial "momentum" to "push" the transition trajectory in
finite time. If there is no initial "momentum" ($p_{0|T}=0$) then the
transition trajectory will reach the exit point $L$ in infinite time. Since
all paths considered here will start with $\varphi_{0|T}=0$, it follows from
(58) that the function $\mathcal{H}$ is constant on any circle centered at
$p=0$; i.e., $\mathcal{H}(0,p)=H$ for $p\in\mathbb{R}^{2}$ such that
$\|p\|^{2}=H$.

Rather than solve (59) numerically, we can take advantage of the linear OU
dynamics and solve it exactly. One can verify that the ($T<\infty$) suboptimal
MLTs are given by

$\displaystyle\varphi_{t|T}$ $\displaystyle=\frac{e^{\mu t}-e^{-\mu t}}{e^{\mu
T}-e^{-\mu T}}R_{t-T}L,$ (61) $\displaystyle p_{t|T}$ $\displaystyle=\frac{\mu
e^{\mu t}}{e^{\mu T}-e^{-\mu T}}R_{t-T}L,$ (62) $\displaystyle S_{t|T}$
$\displaystyle=\frac{\mu\|L\|^{2}}{2}\left(\frac{e^{-2\mu T}}{\left(1-e^{-2\mu
T}\right)^{2}}\right)(e^{2\mu t}-1).$ (63)

In fact, one can further verify (see Appendix A) that $\varphi_{t|T}$ is the
exact conditional mean $\mathbb{E}[X_{t}\mid X_{T}=L,X_{0}=0]$, which happens
to be independent of $\epsilon$ in this example.

To complement the exact formulae for the ($T<\infty$) suboptimal MLTs
$\varphi_{t|T}$, we next derive the ($T\to\infty$) optimal MLT
$\varphi_{s|\infty}$ as a special case. Recall that the optimal MLT starts at
a critical point of (59), namely $x=p=0$. Hence, it is necessary to consider
the solution converging to the critical point in the limit $t\to-\infty$. We
can write the optimal MLT in backward time $s<0$ so that
$\varphi_{0|\infty}=L$ and $\lim_{s\to-\infty}\varphi_{s|\infty}=0$. We find

$\displaystyle\varphi_{s|\infty}$ $\displaystyle=e^{\mu
s}R_{s}L,\quad-\infty<s<0,$ (64) $\displaystyle p_{s|\infty}$
$\displaystyle=\mu e^{\mu s}R_{s}L,$ (65) $\displaystyle S_{s|\infty}$
$\displaystyle=\frac{\mu}{2}e^{2\mu s}\|L\|^{2}.$ (66)

We note that

$V(L)=S_{0|\infty}=\frac{\mu}{2}\|L\|^{2},$

where $V$ is the quasi potential (15).

For our example, it is possible to exactly sample the conditional transition
paths $X^{\epsilon}_{t|T}$. It is well-known that the the OU process over a
finite set of time indices is a multivariate normal random variable [11].
Furthermore if one conditions a multivariate normal random variable on the
values of a subset of the variables, then one obtains another multivariate
normal random variable (with a different mean and covariance matrix) [16]. Let
$U_{ts|T}\in\mathbb{R}^{2\times 2}$ be the scaled conditional auto covariance
matrix, i.e.

$\displaystyle U_{ts|T}$
$\displaystyle=\epsilon^{-1}\mathbb{E}\left[X_{t}X_{s}^{T}\mid
X_{T}=L,\,X_{0}=0\right]$ (67) $\displaystyle=U_{ts}-U_{tT}U_{TT}^{-1}U_{Ts},$
(68)

where the scaled auto covariance matrix is given by, [7, Ch. 4, p. 109]

$\displaystyle U_{sr}$
$\displaystyle=2\epsilon^{-1}\mathbb{E}\left[X_{s}X_{r}^{T}\right]$ (69)
$\displaystyle=\frac{2}{\mu}\left(e^{-\mu|s-r|}-e^{-\mu|s+r|}\right)R_{s-r}.$
(70)

Suppose we discretize $[0,T]$ into $n+2$ equally space points $t_{i}$,
$i=0,\ldots,n+1$, where $t_{0}=0$ and $t_{n+1}=T$.

One would like to illustrate the analytic results of this paper by sampling
$N$ iid copies of the single walker hitting time (55). However this is not
computationally feasible since $N\gg 1$. However we can still accurately
compute random MLT trajectories by exploiting the known Gaussian expression
for the conditional Ornstein-Uhlenbeck process. To do this we discretize the
times upto time $T$ into $n$ points. We then define the $2n$-dimensional
random variable representing a sampled conditional path with

$X_{1:n|T}^{\epsilon}=\\{X_{t_{1}|T}^{\epsilon},\ldots,X_{t_{n}|T}^{\epsilon}\\}.$
(71)

Since (61) is the conditional mean, combined with (67) they define

$X_{1:n|T}^{\epsilon}:\sim\text{MultivariateNormal}(\varphi_{1:n|T},\epsilon
U_{1:n,\,1:n|T}),$ (72)

where $U_{1:n,\,1:n|T}$ can be thought of as a $2n\times 2n$ symmetric matrix
according to (67) ( $U_{1:n,\,1:n|T}$ is known as the conditional covariance
matrix). The most likely first-hitting-time can be determined using the
rigorous relationship, $\epsilon\log(N)=V(L,T)$ between $T$ and $N$ to
illustrate our results (as proved in the previous section).

We also note that for an OU process (see Appendix A for derivation), the cost
function (12) and solution to the Hamilton-Jacobi equation (14) is given by
the quadratic form

$\displaystyle
V(x,t)=\frac{1}{2}x^{T}U_{tt}^{-1}x=\frac{\mu}{2}\|x\|^{2}\left(\frac{1}{1-e^{-2\mu
t}}\right).$ (73)

Figure 1: Single walker rare event. (A) Samples of the optimal escape
trajectory (72) (red) for $\epsilon=10^{-3}$, with streamlines showing the
deterministic flows ($\dot{x}=Ax$) for reference. (The optimal trajectory
$T=\infty$ is approximated with $T=30$.) (B) The escape trajectories follow a
"tube". As the noise magnitude $\epsilon$ goes to zero, the width of the
"tube" goes to zero. (C) sampled conditional transition paths with
$\epsilon=5\times 10^{-3}$ are shown as functions of time. The paths are
separated into two segments or phases. The lag phase is shown as red and the
transition phase is shown as orange. Other parameters are $\mu=0.25$,
$\omega=3$.

To illustrate the exit trajectories, we focus on a single exit point
$L=(1,1)$, but we note that any of the points on the circle
$\|L\|=\text{const}$ (dashed line in Fig. 1A) are equally likely.

Several samples of conditional transition trajectories (72) are shown in Fig.
1. As the noise amplitude is decreased, the sampled trajectories condense into
a tube like path, the center of which is the MLT (61). Several sample
conditional paths are shown in Fig. 1C illustrating the two phases of a given
rare event path: the exponentially long lag phase (shown in red) and the
transition phase (shown in orange). Note that in general, the lag phase is
likely to be much longer than shown.

Recall that we assume the distinguished limit $\epsilon\to 0$, $N\to\infty$
with $\epsilon\log(N)=\Lambda$. In regime 1, there is a labor shortage. If
$\Lambda<V(L)$ then there are not enough walkers to overcome the exponentially
long (in $\epsilon^{-1}$) lag phase (see Fig. 1C) and
$\mathbb{E}[\tau^{N}_{\rm first}]\to\infty$ as $\epsilon\to 0$. The MLT is
essentially described by the globally optimal $T\to\infty$ MLT (64) (see Fig.
1A,B). Increasing the number of walkers in this regime has the strongest
effect on reducing the extreme first hitting time; that is,
$\mathbb{E}[\tau^{N}_{\rm first}]\propto N^{-1}$. One can define a critical
labor supply,

$N_{*}=\exp(V(L)/\epsilon)=\exp\left(\frac{\mu}{2\epsilon}\|L\|^{2}\right),$

which serves as a boundary between the regime 1 and 2.

Figure 2: Trajectories in regime 2: balanced labor. In Regime 2, the MLTs are
deformed as $N$ is increased ($T$ is decreased), with
$N=\exp(S_{T|T}/\epsilon)$, $\epsilon=10^{-2}$. (A) Multiple suboptimal MLTs
conditioned on the final time $T$. (B) Sampled transition paths for
$\epsilon=10^{-2}$. For comparison, the trajectories are aligned at the end
and shown backward in time. Other parameters are $\mu=0.25$, $\omega=3$.

In regime 2, the labor supply is balanced against the single-walker rare event
time $\tau$. If $V(L)<\Lambda<S_{T|T}<\infty$ then $N>N_{*}$ and there is a
sufficient number of walkers to accomplish the rare event in finite time as
$\epsilon\to 0^{+}$. The MLT duration $T$ is finite and does not follow the
same path as the globally optimal MLT. We have that $\mathbb{E}[\tau^{N}_{\rm
first}]\sim T=\Psi_{L}(\Lambda)$, where $\Psi_{L}(\Lambda)$ is given by (16).
The conditional mean $\varphi_{t|T}$ (61) is plotted in Fig. 2A for several
values of $N$ (We also note the corresponding value of $T$). Sample paths are
shown in Fig. 2B as functions of time. Increasing the number of walkers is
equivalent to reducing $T$, which results in an MLT that moves more directly
toward the exit point $L$. For exceedingly large $N$ (e.g., $N=10^{223}$ or
equivalently $T=0.1$), the path is nearly a straight line (see the black
curves in Fig. 2).

In regime 3, there is a labor surplus. As $\Lambda\to\infty$ we reach a regime
where adding more walkers has the smallest effect on the extreme first hitting
time; namely $\mathbb{E}[\tau^{N}_{\rm first}]\propto 1/\log(N)$. One way we
could have the limit $\Lambda\to\infty$ is taking $N\to\infty$ with $\epsilon$
fixed. We know that we have a labor surplus when the cost function begins to
scale like $S_{T|T}\sim MT^{-1}$, as $T\to 0^{+}$, for some constant $M>0$.
The small $T$ asymptotics of the cost function has been studied for a large
class of nonlinear SDEs [5, Ch. 11, p. 407], and it is known that the constant
$M$ depends on the Riemannian length between the origin and $L$. With this
scaling, we have that $\mathbb{E}[\tau^{N}_{\rm
first}]=T\sim\frac{M}{\epsilon\log(N)}$, as $T\to 0^{+}$. We therefore
conclude that the crossover timescale can be determine by the precise small
$T$ asymptotics of $S_{T|T}$. In our simple example, expanding (63) around
$T=0$ yields

$S_{T|T}\sim\frac{\|L\|^{2}}{4}\left(T^{-1}+\mu+O(T)\right),\quad T\to 0^{+}.$
(74)

Notice that the leading order term is independent of $\mu$ and $\omega$. The
above expansion suggests that the crossover timescale is $T\ll 1/\mu$ which in
turn identifies,

$N>N_{**}=\exp\left(\frac{V(L,\mu^{-1})}{\epsilon}\right)$ (75)

where the scaling of the extreme hitting time is $1/\log(N)$. The factor
$M/\epsilon=\|L\|^{2}/(4\epsilon)$ is the mean first passage time for a single
Brownian walker to reach the boundary of the circle of radius $L$, starting
from the center. It is this time that is reduced by a factor of $1/\log(N)$ as
$N\to\infty$ so that

$\mathbb{E}[\tau^{N}_{\rm
first}]\sim\frac{\left(\frac{\|L\|^{2}}{4\epsilon}\right)}{\log(N)},\quad
N\to\infty.$ (76)

##  5 Discussion

Our results start with the establishment of a connection between the number of
walkers, $N$, and the duration of the MLT, $T$. This means that adding more
walkers is equivalent to conditioning on a faster MLT. It has long been known
in the Large Deviation Theory literature that, given a fixed exit point, the
globally optimal MLT is given by $T\to\infty$, and it is these paths that
determine the asymptotics of the single walker mean exit time. Hence, the
$T\to\infty$ MLTs have received the most attention in the literature. By
establishing a link between $N$ and $T$, we have identified a new application
for the sub-optimal, finite $T$, MLTs.

The example in Section 4 involved a linear SDE. For general nonlinear
problems, the best course of action is to employ numerical methods [10, 3, 27]
and dynamical systems techniques [20, 21, 22, 27] to study the problem. The
deterministic ($\epsilon=0$) dynamics should be studied first. The stochastic
dynamics can then be understood through (17). Notable features that can occur
in general nonlinear problems include the possibility that optimal
($T\to\infty$) MLTs are not unique [21, 27] and the behavior of MLTs near
unstable critical points of the deterministic dynamics [22, 27]— particularly
saddles, which often serve as boundaries between basins of attraction.

Finally, we note that our choice of using a continuous Markov process is
purely for the sake of brevity. The theory of large deviations, including the
characterization of rare events in terms of MLTs, has been extended to many
other types of Markov process, including continuous time jump processes [5],
piece-wise deterministic Markov processes [5, 26, 2], mixtures [27, 25, 35],
and spatially extended Markov processes [4]. We emphasize that the results
presented here should be generalized by applying LDT directly to the full
process in question and not by employing an intermediate diffusion
approximation. This is particularly true for the case where fast variables are
averaged out [5, ch. 7], since this constrains the resulting MLT and leads to
a characterization of the rare event using suboptimal trajectories [29, 28,
24]. A general strategy for applying a large deviation analysis to a given
weak-noise Markov process is to start by calculating the Hamiltonian function
analogous to Eqn. (10). The most straightforward way to accomplish this is
with a WKB expansion [27]. Once the Hamiltonian is known, one can utilize
numerical methods to solve the Hamilton-Jacobi equation (14) via (17).

## References

  * [1] Alan J Bray, Satya N Majumdar, and Grégory Schehr.  Persistence and first-passage properties in nonequilibrium systems.  Advances in Physics, 62(3):225–361, 2013.
  * [2] Paul C Bressloff and Jay M Newby.  Path integrals and large deviations in stochastic hybrid systems.  Physical Review E, 89(4):042701, 2014.
  * [3] MK Cameron.  Finding the quasipotential for nongradient sdes.  Physica D: Nonlinear Phenomena, 241(18):1532–1550, 2012.
  * [4] Amir Dembo and Ofer Zeitouni.  Large deviations techniques and applications. Second Edition, volume 38.  Springer Science & Business Media, 2009.
  * [5] Mark I Freidlin and Alexander D Wentzell.  Random Perturbations of Dynamical Systems, volume 260.  Springer Science & Business Media, 2012.
  * [6] Mark Iosifovich Freidlin, Alexander D Wentzell, MI Freidlin, and AD Wentzell.  Random perturbations.  Springer, 1998.
  * [7] Crispin Gardiner.  Stochastic methods, volume 4.  Springer Berlin, 2009.
  * [8] Tobias Grafke and Eric Vanden-Eijnden.  Numerical computation of rare events via large deviation theory.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(6):063118, 2019.
  * [9] Magnus J Haughey, Aleix Bassolas, Sandro Sousa, Ann-Marie Baker, Trevor A Graham, Vincenzo Nicosia, and Weini Huang.  First passage time analysis of spatial mutation patterns reveals sub-clonal evolutionary dynamics in colorectal cancer.  PLOS Computational Biology, 19(3):e1010952, 2023.
  * [10] Matthias Heymann and Eric Vanden-Eijnden.  The geometric minimum action method: A least action principle on the space of curves.  Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 61(8):1052–1117, 2008.
  * [11] Ioannis Karatzas, Steven E Shreve, Ioannis Karatzas, and Steven E Shreve.  Brownian motion.  Brownian motion and stochastic calculus, pages 47–127, 1998.
  * [12] Joseph B Keller.  Diffusion at finite speed and random walks.  Proceedings of the National Academy of Sciences, 101(5):1120–1122, 2004.
  * [13] Venu Kurella, Justin C Tzou, Daniel Coombs, and Michael J Ward.  Asymptotic analysis of first passage time problems inspired by ecology.  Bulletin of mathematical biology, 77(1):83–125, 2015.
  * [14] Sean D Lawley.  Distribution of extreme first passage times of diffusion.  Journal of Mathematical Biology, 80(7):2301–2325, 2020.
  * [15] Sean D Lawley and Jacob B Madrid.  A probabilistic approach to extreme statistics of brownian escape times in dimensions 1, 2, and 3.  Journal of Nonlinear Science, 30(3):1207–1227, 2020.
  * [16] Georg Lindgren, Holger Rootzén, and Maria Sandsten.  Stationary stochastic processes for scientists and engineers.  CRC press, 2013.
  * [17] Samantha Linn and Sean D Lawley.  Extreme hitting probabilities for diffusion.  Journal of Physics A: Mathematical and Theoretical, 55(34), 2022\.
  * [18] Donald Ludwig.  Persistence of dynamical systems under random perturbations.  Siam Review, 17(4):605–640, 1975.
  * [19] Jacob B Madrid and Sean D Lawley.  Competition between slow and fast regimes for extreme first passage times of diffusion.  Journal of Physics A: Mathematical and Theoretical, 53(33):335002, 2020.
  * [20] Robert S Maier and Daniel L Stein.  Escape problem for irreversible systems.  Physical Review E, 48(2):931, 1993.
  * [21] Robert S Maier and Daniel L Stein.  A scaling theory of bifurcations in the symmetric weak-noise escape problem.  Journal of statistical physics, 83:291–357, 1996.
  * [22] Robert S Maier and Daniel L Stein.  Limiting exit location distributions in the stochastic exit problem.  SIAM Journal on Applied Mathematics, 57(3):752–790, 1997.
  * [23] Victor Matveev.  Comparison of deterministic and stochastic approaches for calcium dependent exocytosis.  Biophysical Journal, 114(3):283a–284a, 2018.
  * [24] J Newby and J Chapman.  Metastable behavior in markov processes with internal states: breakdown of model reduction techniques.  J. Math. Biol., 10, 2013.
  * [25] Jay Newby.  Bistable switching asymptotics for the self regulating gene.  Journal of Physics A: Mathematical and Theoretical, 48(18):185001, 2015.
  * [26] Jay M Newby.  Isolating intrinsic noise sources in a stochastic genetic switch.  Physical biology, 9(2):026002, 2012.
  * [27] Jay M Newby.  Spontaneous excitability in the morris–lecar model with ion channel noise.  SIAM Journal on Applied Dynamical Systems, 13(4):1756–1791, 2014\.
  * [28] Jay M Newby, Paul C Bressloff, and James P Keener.  Breakdown of fast-slow analysis in an excitable system with channel noise.  Physical Review Letters, 111:128121, 2013.
  * [29] Jay M Newby and James P Keener.  An asymptotic analysis of the spatially inhomogeneous velocity-jump process.  Multiscale Modeling & Simulation, 9(2):735–765, 2011.
  * [30] Lukas Ramlow, Martin Falcke, and Benjamin Lindner.  An integrate-and-fire approach to ca2+ signaling. part i: Renewal model.  Biophysical Journal, 122(4):713–736, 2023.
  * [31] Jürgen Reingruber and David Holcman.  Computational and mathematical methods for morphogenetic gradient analysis, boundary formation and axonal targeting.  Seminars in Cell & Developmental Biology, 35:189–202, 2014.  Regulated Necrosis & Modeling developmental signaling pathways & Development of connective maps in the brain.
  * [32] Z Schuss, K Basnayake, and D Holcman.  Redundancy principle and the role of extreme statistics in molecular and cellular biology.  Physics of life reviews, 28:52–79, 2019.
  * [33] Gregory D Smith, Joel E Keizer, Michael D Stern, W Jonathan Lederer, and Heping Cheng.  A simple numerical model of calcium spark formation and detection in cardiac myocytes.  Biophysical journal, 75(1):15–32, 1998.
  * [34] Srinivasa RS Varadhan.  Diffusion processes in a small time interval.  Communications on Pure and Applied Mathematics, 20(4):659–685, 1967\.
  * [35] Benjamin L Walker and Katherine A Newhall.  Numerical computation of effective thermal equilibria in stochastically switching langevin systems.  Physical Review E, 105(6):064113, 2022.
  * [36] George H. Weiss and Menachem Dishon.  On the asymptotic behavior of the stochastic and deterministic models of an epidemic.  Mathematical Biosciences, 11(3):261–265, 1971.
  * [37] George H Weiss, Kurt E Shuler, and Katja Lindenberg.  Order statistics for first passage times in diffusion processes.  Journal of Statistical Physics, 31:255–278, 1983.

##  Appendix A Appendix

We provide some extra details on the Ornstein-Uhlenbeck process of Section 4.
It is well-known that $X_{t}$ is Gaussian [11, Section 5.6], and hence one can
immediately infer that, writing $B_{\delta}(L)\subset\mathbb{R}^{2}$ to be a
closed ball of radius $\delta$ about $L$,

$\displaystyle\epsilon\log\mathbb{P}\big{(}x_{t}\in
B_{\delta}(L)\big{)}=-\frac{1}{4}L^{T}U_{tt}^{-1}L+o(\epsilon,\delta),$ (77)

where the corrections terms are negligible in taking $\epsilon\to 0^{+}$
(first), and then $\delta\to 0^{+}$ (second). We now claim that the above
quadratic form is in fact the cost-function, and that the solution of (59) -
(60) is in fact just the conditional mean (defined just below in (78).

To see this, for any $s\in[0,t]$, define $m_{s|t}\in\mathbb{R}^{2}$ to be the
conditional mean, that is

$\displaystyle m_{s|t}$ $\displaystyle=\mathbb{E}\left[x_{s}\mid
x_{t}=L,\,x_{0}=0\right]$ (78) $\displaystyle=U_{st}U_{tt}^{-1}L$ (79)

In fact one easily checks that both $m_{s|t}$ and $U_{rs|t}$ are independent
of $\epsilon$.

Next, we claim that the cost-function assumes the form (for any
$L\in\mathbb{R}^{2}$ and any $t>0$)

$\displaystyle V(L,t)=\frac{1}{2}L^{T}U_{tt}^{-1}L.$ (80)

To see why this holds, for any finite subset $\pi\subset[0,t]$, define the
event

$\displaystyle\mathcal{A}_{\pi,\tilde{\delta}}=\big{\\{}y\in\mathcal{C}([0,t],\mathbb{R}^{2}),\;\sup_{s\in\pi}\|y_{s}-m_{s|t}\|\geq\tilde{\delta}\big{\\}}$
(81)

We can thus conclude that for any $\tilde{\delta}>0$,

$\underset{\delta\to 0}{\overline{\lim}}\underset{\epsilon\to
0}{\overline{\lim}}\epsilon\log\mathbb{P}\big{(}x_{t}\in
B_{\delta}(L)\big{)}-\underset{\delta\to
0}{\overline{\lim}}\underset{\epsilon\to
0}{\overline{\lim}}\epsilon\log\mathbb{P}\big{(}x_{t}\in
B_{\delta}(L),x\in\mathcal{A}_{\pi,\tilde{\delta}}\big{)}\\\
=\frac{1}{2}\sup_{y\in\mathcal{A}_{\pi,\tilde{\delta}}}\sum_{r,s\in\pi}\sum_{a,b\in\\{1,2\\}}y_{r}(a)y_{s}(b)\tilde{V}_{rs}(a,b)>0$
(82)

where
$\tilde{V}=\big{(}\tilde{V}_{rs}(a,b)\big{)}_{r,s\in\pi,a,b\in\\{1,2\\}}\in\mathbb{R}^{2|\pi|\times
2|\pi|}$ is the covariance matrix with elements
$\big{(}U_{rs|t}(a,b)\big{)}_{r,s\in\pi,a,b\in\\{1,2\\}}$. Indeed the RHS of
(82) is strictly positive because $\tilde{V}$ must be strictly positive
definite (the inverse of a positive definite symmetric matrix is codiagonal
with the original matrix, with each eigenvalues the inverse of one of the
original eigenvalues).

It follows from (82) that

$\displaystyle\lim_{\delta\to 0^{+}}\lim_{\epsilon\to
0^{+}}\epsilon\log\mathbb{P}\big{(}x_{t}\in
B_{\delta}(L)\big{)}=\lim_{\delta\to 0^{+}}\lim_{\epsilon\to
0^{+}}\epsilon\log\mathbb{P}\big{(}\sup_{s\in[0,t]}\|x_{s}-\phi_{s|t}\|\leq\delta\big{)}.$
(83)

Because the Large Deviations Rate Function $I_{T}$ is lower-semi-continuous,
this implies that

$\displaystyle\inf_{z\in H_{1}\;:z_{t}=L}I_{T}(z)=I_{T}(\phi)$ (84)

which immediately implies (80).

It is also known that [11] the covariance matrix satisfies the ODE

$\displaystyle\frac{d}{dt}U_{tt}$
$\displaystyle=AU_{tt}+U_{tt}A^{T}+2\mathbb{I}$ (85) $\displaystyle U_{00}$
$\displaystyle=0\mathbb{I}.$ (86)

Because the eigenvalues of $A$ have negative real parts, there must exist a
unique invariant distribution to the stochastic dynamics (54). This
distribution is centered and Gaussian, with covariance matrix
$U_{*}\in\mathbb{R}^{2\times 2}$ satisfying the identity

$\displaystyle AU_{*}+U_{*}A^{T}+2\mathbb{I}=0.$ (87)

That is, there is a unique positive-definite symmetric matrix $U_{*}$
satisfying the above identity, and $U_{tt}$ converges to $U_{*}$ as
$t\to\infty$. We thus find that the quasi-potential is

$\displaystyle V(L)=\frac{1}{2}L^{T}U_{*}^{-1}L.$ (88)


# A Survey of Knowledge Graph Embedding and Their Applications

Shivani Choudhary, Tarun Luthra, Ashima Mittal, Rajat Singh
Indian Institute of Technology Delhi
Hauz Khas, Delhi-110016, India
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

Knowledge Graph embedding provides a versatile technique for representing
knowledge. These techniques can be used in a variety of applications such as
completion of knowledge graph to predict missing information, recommender
systems, question answering, query expansion, etc. The information embedded in
Knowledge graph though being structured is challenging to consume in a real-
world application. Knowledge graph embedding enables the real-world
application to consume information to improve performance. Knowledge graph
embedding is an active research area. Most of the embedding methods focus on
structure-based information. Recent research has extended the boundary to
include text-based information and image-based information in entity
embedding. Efforts have been made to enhance the representation with context
information. This paper introduces growth in the field of KG embedding from
simple translation-based models to enrichment-based models. This paper
includes the utility of the Knowledge graph in real-world applications.

##  1 Introduction

Knowledge graph(KG) has received a lot of traction in the recent past and led
to much research in this area. Most of the research is focused on the
generation of Knowledge graphs and consumption of the information enshrined in
the Knowledge Graph. Some of the earlier works to create KG are YAGO(Suchanek
et al., 2007), Freebase(Bollacker et al., 2008), DBpedia(Lehmann et al., 2012)
and WikiData(Vrandečić and Krötzsch, 2014). The evolution of the Knowledge
Graph starts with the seminal paper from Berners-Lee (Berners-Lee et al.,
2001). The knowledge graph has evolved in three phases. In the first phase,
the Knowledge representation was brought to the level of Web standard. The
core focus shifted to Data management, linked data, and its application in the
second phase. In the third phase, the focus shifted on the real-world
application (Bonatti et al., 2018). The real-world application ranges from
Semantic parsing(Berant et al., 2013; Heck and Huang, 2014), recommender
system(Sun et al., 2020; Wang et al., ), question answering(Saxena et al., ),
named entity disambiguation(Lin et al., ), information extraction (Xiong and
Callan, ; Xiong and Callan, ; Liu et al., 2018; Dietz et al., 2019) etc. The
knowledge graph is a representation of structured relational information in
the form of Entities and relations between them. It is a multi-relational
graph where nodes are entities and edges are relations. Entities are real-
world objects or abstract information. The representation entities and
relations between them are represented as triple. For e.g. (New Delhi,
IsCapitalOf, India) is an example of a triple. New Delhi and Delhi are
entities and IsCapitalOf is relation. Though the representation looks
scientific but consuming these in the real-world application is not an easy
task. The consumption of information enshrined in the knowledge graph will be
very easy if it can be converted to numerical representation. Knowledge graph
embedding is a solution to incorporate the knowledge from the knowledge graph
in a real-world application.
The motivation behind Knowledge graph embedding (Bordes et al., ) is to
preserve the structural information, i.e., the relation between entities, and
represent it in some vector space. This makes it easier to manipulate the
information. Most of the work in Knowledge graph embedding (KGE) is focused on
generating a continuous vector representation for entities and relations and
apply a relationship reasoning on the embedding. The relationship reasoning is
supposed to optimize some scoring function to learn the embedding. Researchers
have used different approaches to learn the embedding Path-based
learning(Toutanova et al., 2016), Entity-based learning, textual-based
learning, etc. A lot of work was focused were translation model (Bordes et
al., ) and semantic-based model(Bordes et al., 2014). The representation of
the triple results in a lot of information lost because it fails to take
”textual information” into account. With the proposal of Graph attention
network (Velicković et al., 2017) the representation of the entities has
become more contextualized. In recent years, the proposal of multi-modal graph
has extended the spectrum to a new level. In multi-modal knowledge graph,
Knowledge graph can have multi-modal information like image and text etc. (Sun
et al., 2020; Wei et al., 2019).
Previous work on survey has focused on the KG Embedding (Wang et al., ), KG
Embedding and application(Ji et al., 2020) , KG embedding with textual data
(Lu et al., 2020), KG embedding based on the deep-learning (Wang et al.,
2020). This work shall focus on the KG embedding from translation-based model,
semantic-based model, embedding with enriched representation from textual data
and multi-modal data and their application. In section 2, we shall provide the
details of KGE; in section 3, we shall present the application area. In the
summary section, we shall try to put emerging areas of research in KGE.

##  2 Knowledge Graph embedding

Knowledge Graph embedding is an approach to transform Knowledge Graphs (nodes,
edges, and their feature vectors) into a low dimensional continuous vector
space that preserves various graph structure, information, etc. These
approaches are broadly classified into two groups: translation models and
semantic matching models.

###  2.1 Translation Models

The translation-based model uses distance-based measures to generate the
similarity score for a pair of entities and their relationships. The
translation-based model aims to find a vector representation of entities with
relation to the translation of the entities. It maps entities to a low-
dimensional vector space.

####  2.1.1 TransE (Bordes et al., 2014)

The first model proposed was TransE. It is an energy based model. If a triple
h, r, t holds then the vector representation h and r should be as close as
possible. It can be graphically represented as Figure-1. Mathematically, it
can be stated like $\textit{h}+\textit{r}\approx\textit{t}$. The energy of a
triple is d(h + r,t) for some similarity measure d. To learn the embedding,
minimization of ranking based loss function over the training set.

$\mathcal{L=\sum_{\textit{h,r,t}\in
S}\sum_{\textit{$\hat{h}$,r,$\hat{t}$}\in\hat{S}}[\gamma+\textit{d(h+r),t}-\textit{d($\hat{h}$,r,$\hat{t}$)}]}$
(1)

d(ĥ,r,t̂) represents the set of the corrupt triples. This loss function will
be optimized so that the valid triples are ranked above the corrupt triples.
This model fails in case of the one to many relation and many to many
relation. To overcome this deficit new model TransH is proposed.

Figure 1: The intuition of TransE model. Image credit Wang et al.

####  2.1.2 TransH Wang et al. (2014)

It was proposed to address the limitations of TransE. This model enables an
entity to have distributed representations based on their involvement in the
relation. The representation of $h_{\bot}$, $t_{\bot}$ are projected in a
relation specific hyperplane. The relation between then is $d_{r}$. As in the
Figure-2, the vectors h and t are projected in the relation hyper-plane. The
loss function and intuition remain similar to TransE.

Figure 2: TransH model. Image taken from Wang et al.

####  2.1.3 TransR Lin et al. (2015)

It proposes that an entity may have multiple attributes and various relations.
Each relation may focus on different attributes of entities. TransR models
entities and relation in different embedding space. It means that two
different spaces: entity space and relation space are modelled. Each entity is
mapped into relation space. The translation construct is applied on the
projected representation in the relationship space. The Figure-3 presents the
intuition behind TransR model.

Figure 3: TransR model. Image taken from Lin et al. (2015)

To further refine the representation, a new model TransD was proposed by Ji et
al. (2015). TransR captures the possibility of relations and their embedding
from the relation space. But the TransD, extends it to the entity space as
well. Here Entity-relation pair is considered as the first-class object.

####  2.1.4 RotatE Sun et al. (2019a)

In knowledge graphs, often, the relation is symmetric/anti-symmetric,
inversion and composition. e.g., ”Marriage” is a symmetric relation, ”My niece
is my sister’s daughter” is a composition, etc. The models discussed above are
not capable of predicting these relations. The model proposed here is based on
the intuition that the relation from head to tail is modeled as rotation in a
complex plane. It is motivated by Euler’s identity.

$e^{\iota\theta}=\sin\theta+\iota\cos\theta$ (2)

For a triplet (h,r,t), the relation among them can be represented as t = h
$\circ$ r. Where h, r, t is the k-dimension embedding of the head, relation
and tail, restricting the $|{r_{i}}|$ = 1. It means we are in the unit circle
and $\circ$ represents element wise product. For each dimension
$t_{i}=r_{i}h_{i}$ subjected to constraint $|{r_{i}}|$ = 1. Under these
condition the relationship is symmetric $\iff$ for all the values of i,
$e^{0/\iota\theta}=\pm 1$, The relationship is inverse $\iff$
$r_{j}=\bar{r}_{i}$, i.e. both are the complex conjugate. Two relations are
composite $\iff$ $r_{j}=r_{i}\circ r_{k}$. It means the relation $r_{j}$ can
be obtained by a combined rotation of $r_{i}$ and $r_{k}$,
$\theta_{j}=\theta_{i}+\theta_{k}$. The scoring function measures the angular
distance.

$d_{r}(h,t)=\left\|h\circ r-t\right\|$ (3)

The Figure-4 shows the comparison of RotaE and TransE.

Figure 4: Representation of RotaE in comparison to TransE. In RotatE, the
relationship between h and t is represented as angle of rotation. This image
depicts the relationship in 1-dimension space. (Sun et al., 2019a)

####  2.1.5 HakE Zhang et al. (2019)

The approaches we have discussed so far fails to capture the semantic
hierarchies. In HakE, the authors have proposed to model the hierarchy in the
entities as concentric circles in polar coordinate. The entity with smaller
radius belongs to higher up in the hierarchy. The angle between them
represents the variation in the meaning. To represent a point on the circle,
we should have $(r,\theta)$. Similarly, this model has two components, one to
map the modulus and the other one to map the angle. The modulus part considers
the depth of the tree as moduli. Let $h_{m},t_{m}$ and $r_{m}$ are the
representation in modulus space then $h_{m}\circ r_{m}=t_{m}$, where
$h_{m},t_{m}$ is a k-dim vector. The distance function is similar to RotatE
with modification to consider only modulo part.

$d_{r,m}(h_{m},t_{m})=\left\|h_{m}\circ r_{m}-t_{m}\right\|_{2}$ (4)

Similarly, the phase part can be formulated as $(h_{p}+r_{p})mod2\pi=t_{p}$
where $h_{p},r_{p},t_{p}\in[0,2\pi)^{k}$. The distance function is

$d_{r,p}(h_{p},t_{p})=\left\|\sin((h_{p}\circ r_{p}-t_{p})/2)\right\|_{1}$ (5)

By combining both the part, an entity can be mapped in the polar coordinate
space.

###  2.2 Semantic Matching Models

Semantic Matching is one of the core tasks in Natural Language Processing. As
we have seen that Translational Distance Models use distance-based scoring
functions to calculate the similarity between the different entities and
relations thus build the embedding accordingly. On the other hand, Semantic
Matching Models use similarity-based scoring function. There are several
Knowledge Graph Embedding algorithms comes under this model. Some of the
algorithms are described below.

####  2.2.1 RESCAL Nickel et al. (2011)

RESCAL follows Statistical Relational Learning Approach which is based on a
Tensor Factorization model that takes the inherent structure of relational
data into account. A TensorKolda and Bader (2009) is a multidimensional array.
More formally we can say that a first order Tensor is a vector, second order
Tensor is a matrix and, Tensor with more than two order is called as higher
order Tensor. Tensor Factorization is expressing a Tensor as a sequence of
elementary operations acting on other, often a simpler Tensors. Statistical
Relational Learning inherits from Probability Theory and Statistics to address
uncertainty and complexity of relational structures.
Nickel et al. (2011) models the Knowledge Graph triplet of the form (head,
relation, tail) into three-way tensor, $\mathcal{X}$ as shown in Figure 5.

Figure 5: Tensor model for relational data. $E_{1}...E_{n}$ denote the
entities and $R_{1}...R_{m}$ denote the relation in the domainNickel et al.
(2011)

In $\mathcal{X}$, two modes holds the concatenated entities (head and tail),
and the third mode holds relation (relation). A Tensor entity
$\mathcal{X}_{ijk}$ = 1 denotes that there exist a relation and if
$\mathcal{X}_{ijk}$ = 0 denotes that there is unknown relation. It is assumed
that data is given as n * n * m Tensor. Where n is the number of entities and
m is the number of relations. RESCAL ”explains triples via pairwise
interaction of latent features”. It performs the rank-r factorization on each
slice of $\mathcal{X}$ (relational data) and the score of a fact (head,
relation, tail ) is given by the following bi-linear function.

$f_{r}(h,t)=\textbf{h}^{T}\textbf{M}_{r}\textbf{t}=\sum_{i=0}^{d-1}\sum_{j=0}^{d-1}[\textbf{M}_{r}]_{ij}.[\textbf{h}]_{i}.[\textbf{t}]_{j}$

where h,t $\in$ $\mathbb{R}^{d}$ are vector representation of entities, and
$M_{r}$ $\in$ $\mathbb{R}^{d*d}$ is a matrix representation of $r^{th}$
relation. Thus from this equation we are able to calculating the score of the
triple using the weighted sum of all the pairwise interactions between the
latent features of the entities $h$ and $t$ as shown in the Figure 6.

Figure 6:  RESCAL Representation. Here number of latent features for entities
are 3 and number of latent features for relations are 3. Nickel et al. (2016a)

This method require $O(d^{2})$ parameters per relation and the space
complexity of $O(nd+md^{2})$ where n is the number of entities and m is the
number of relations.

####  2.2.2 TATEC García-Durán et al. (2014)

TATEC stands for Two And Three-way Embeddings Combination. The main
disadvantage of RESCAL is that it is a Three-way model which performs fairly
good for relationships which occur frequently but it performs poor for the
rare relationships and leads to major over-fitting. The issue of major over-
fitting for rare relationships can be controlled by regularizing or reducing
the expressivity of the model and former method is not feasible. The second
method of reducing the expressivity is Two-way interaction which is
implemented in TransE and SME. Two-way interaction approaches overperform the
Three-way approaches on many datasets from which we can conclude that Two-way
interactions are more efficient for the datasets and specially for those
datasets which have more rare relationships. But the problem with the two-way
interaction is that they are limited and are not able to represent all kind of
relations with entities.
TATEC is a latent factor model which is capable of incorporating the high
capacity Three-way model with well-controlled two-way interactions and take
the advantage of both of them. Since two-way and three-way models do not use
the the same kind of data pattern and do not encode the same kind of
information in the embedding. So, in TATEC during first stage they used two
different embeddings and then combined and fine-tuned them in the later stage.
The scoring function of TATEC is given by
$f_{r}(h,t)=f_{r}^{1}(h,t)+f_{r}^{2}(h,t)$ which is a linear combination of
bi-gram and tri-gram terms, where $f_{r}^{1}(h,t)$ is a two-way interaction
score and $f_{r}^{2}(h,t)$ is a three-way interaction score. These can be
calculates as follows:
1) Two-way interactions terms can be given by:

$f_{r}^{1}(h,t)=\textbf{h}^{T}\textbf{r}+\textbf{t}^{T}\textbf{r}+\textbf{h}^{T}\textbf{D}\textbf{t}$

where D is the diagonal matrix shared across all the different relations and
does not depend on input triple and r $\in$ $\mathbb{R}^{d}$ is a vector that
depends on relationships.
2) Three-way interactions terms can be given by:

$f_{r}^{2}(h,t)=\textbf{h}^{T}\textbf{M}_{r}\textbf{t}$

The final scoring function of TATEC is given by

$f_{r}(h,t)=\textbf{h}^{T}\textbf{M}_{r}\textbf{t}+\textbf{h}^{T}\textbf{r}+\textbf{t}^{T}\textbf{r}+\textbf{h}^{T}\textbf{D}\textbf{t}$

Figure 7:  Link prediction results García-Durán et al. (2014)

Authors of García-Durán et al. (2014) compared this model to the other
existing models such as RESCAL Nickel et al. (2011), TransE (Bordes et al., ),
LFM, SE and, SME for link prediction on FB15k dataset and as a result TATEC
performs better than all other available models as shown in Figure 7.

Time complexity and the space complexity of TATEC is same as RESCAL as TATEC
extends RESCAL. The time complexity of TATEC is $O(d^{2})$ parameters per
relation and the space complexity of $O(nd+md^{2})$ where n is the number of
entities and m is the number of relations.

####  2.2.3 DistMult Yang et al. (2014)

This model compares with the NTN neural model, TransE and bi-linear models
like RESCAL. The problem with NTN is that it is the most expensive model it
incorporate both linear and bilinear relation operations. Similarly TransE
parameterizes the linear operations with one dimensional vectors. DistMult is
a simplified RESCAL, which uses the basic bilinear scoring function.

$f_{r}(h,t)=\textbf{h}^{T}\textbf{M}_{r}\textbf{t}$

These bilinear formulations are combined with different forms of
regularization to make different models. In DistMult authors considered a
simpler approach where they reduced the number of parameters by imposing
restrictions on $M_{r}$ to be a diagonal matrix. This results in a simpler
model and this model enjoys the same scalable properties of TransE as well as
it achieves better performance over TransE. Thus the final scoring function is
given as

$f_{r}(h,t)=\textbf{h}^{T}diag(\textbf{r})\textbf{t}=\sum_{i=0}^{d-1}[\textbf{r}]_{i}.[\textbf{h}]_{i}.[\textbf{t}]_{i}$

where for each relation r, r $\in\mathbb{R}^{d}$ is a vector that depends on
relationships.

Figure 8:  DistMult simple Illustration

Time complexity and the space complexity of DistMult is more efficient as
compared to RESCAL or TATEC. The time complexity of DistMult is $O(d)$
parameters per relation and the space complexity of $O(nd+md)$ where n is the
number of entities and m is the number of relations. Due to its over-
simplified mature of the model, this model is not enough powerful for the use
in case of general Knowledge Graphs because it is only able to work
efficiently with symmetric relations.

####  2.2.4 HolE Nickel et al. (2016b)

HolE stands for Holographic Embedding. HolE tried to overcome the problem of
Tensor Product used in RESCAL by using circular correlation. Tensor product
uses pairwise multiplicative interactions between feature vectors which
results in increase in dimensionality of the representation i.e.,
$\mathbb{R}^{d^{2}}$ thus increase the computational demand.

$\textbf{h}\circ\textbf{t}=\textbf{h}\otimes\textbf{t}\in\mathbb{R}^{d^{2}}$

Where a,b $\in\mathbb{R}^{d}$ are entity embeddings. Tensor products are very
rich in capturing the interactions but are computational intensive. On the
other hand HolE use Circular Correlation which can be seen as compression of
the Tensor Product. The main advantage of Circular Correlation over Tensor
Product is that it won’t increase the dimensionality of the representation.

$\textbf{h}\circ\textbf{t}=\textbf{h}*\textbf{t}\in\mathbb{R}^{d}$

Where * : $\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$
denotes the circular correlation.

$[\textbf{h}*\textbf{t}]_{i}=\sum_{k=0}^{d-1}[\textbf{h}_{k}].\textbf{t}_{(k+i)mod(d)}$

The final score of the fact in HolE is given by matching the compositional
vector ($\textbf{h}*\textbf{t}$) with the relational representation,i.e.,

$f_{r}(h,t)=\textbf{r}^{T}(\textbf{h}*\textbf{t})=\sum_{i=0}^{d-1}[\textbf{r}]_{i}\sum_{k=0}^{d-1}[\textbf{h}]_{k}.[\textbf{t}]_{(k+i)mod(d)}$

Figure 9:  HolE simple Illustration (It requires only d components) Nickel et
al. (2016b)

So HolE is more efficient as compared to RESCAL or TransE. HolE take
$O(dlogd)$ parameters per relation and the space complexity of $O(nd+md)$
where n is the number of entities and m is the number of relations. Another
advantage of HolE is that Circular Correlation is not commutative (
$\textbf{h}*\textbf{t}\neq\textbf{t}*\textbf{h}$ ) thus HolE is able to model
asymmetric relations (directed graphs) with compositional representations
which is not possible in RESCAL.

####  2.2.5 ComplExTrouillon et al. (2016)

Knowledge graphs represent the relation between entities. The entities may be
termed as subjects and objects respectively of a given relation. However not
all relations may be present in a given KG. One of the application of KG is
the ability to predict missing relations or entities.
Dot product of vector embedding of KG triplets is being successfully used for
symmetric, reflexive, anti-reflexive and even transitive relations (Bouchard
et al., 2015) however it can’t be used for anti-symmetric relations. For
example the relation capitalOf(New Delhi, India) is not symmetric since we
cannot interchange subject and object entity in this relation therefore we
need to have different embedding for an entity as subject and as object which
increases the number of parameters.
Complex embedding facilitates joint learning of subject and object entities
while preserving the asymmetry of the relation. It uses Hermitian dot product
of embedding of subject entities and object entities. The Eigen Vector
decomposition is used to identify a low rank diagonal matrix W such that there
exists X = Re($EWE^{T}$) such that X has same sign pattern as Y. The low rank
diagonal matrix W is then used to predict missing relations by applying
Re($e_{s}^{T}W\bar{e}_{o}$).

####  2.2.6 ANALOGYLiu et al. (2017)

ANALOGY is based on a multiplicative model where the relation (s,r,o) is
scored by multiplying the vector representations of of subject (s), relation
(r) and object (o). The relation present in the knowledge graph are expected
to have higher score i.e. $\phi(s,r,0)$ = $e_{s}^{T}W_{r}e_{o}$ will be high
if the triplet (s,r,o) exists.
An example of analogy may be branch:tree::petal:leaf, in such an analogy, the
relation $W_{r}$ ”is part of” may be used to predict the missing entity from
the analogy. The foundation of ANALOGY model is the linear maps of matrix
representation of relations from the triplets present in the KG. Basically,
the model uses the fact that there can be multiple paths to arrive at the
entity from a given entity through a linear map of relations
$W_{r1},W_{r2},W_{r3}$ and so on, and application of such relations in any
sequence will give the same result. Hence such linear maps may be used to
predict entities missing from the knowledge graph.
For example Let’s consider the entities teacher (t), school (s), professor (p)
and college (c). We have two relations in this set up: teachesAt (t,s) &
teachesAt(p,c) and juniorOf(t,p) & juniorOf(s,c) therefore the relation
between teacher and college is teachesAt*juniorOf = juniorOf*teachesAt.
However such linear maps are feasible only if the commutative property holds
for such relations. ANALOGY uses such linear maps to predict entities.

###  2.3 Enrichment based embedding

In the recent times, new emerging research areas are focusing on
contextualized embedding. Under this, the entity under the consideration is
enriched information from the neighbourhood information. A few notable
approaches are Graph attention network (GAT) (Velicković et al., 2017) based
information enrichment. Two methods based on the KGAT (Wei et al., 2019),
MMGAT (Sun et al., 2020) has proposed models to embed contextual information
for an entity. Under MMGAT, they proposed model to combine the embedding from
multi-modal data with an attention framework, as adopted from GAT’s attention
framework. Both the frameworks, were using translation model to learn the
representation after the enrichment. The new emerging research areas are try
to learn the structural information as well as path based information, multi-
modal data. There are other research areas in the embedding are: Text-enhanced
embedding, Logic-enhanced embedding, Image-enhanced embedding (Bianchi et al.,
2020) etc.

##  3 Applications of Knowledge graph embedding

There are many applications of KG embedding learning methods. This section
explores three of them namely- link prediction, triplet classification and
recommender systems.

The first two are In-KG applications, which are conducted within the scope of
the KG. The last is an example of Out-of-KG applications that scale to broader
domains (Wang et al., 2017).

###  3.1 Link Prediction

The set of edges in a knowledge graph is a subset of
Entities$\times$Relations$\times$Entities. The link prediction task focuses on
finding an entity that can be represented as a fact (edge) together with a
given relation and entity i.e., (entity, relation, ?) or (?, relation, entity)
where ? refers to the missing entity. For e.g. (New Delhi, isCapitalOf, ?) or
(?, isCapitalOf, India). Link prediction is a way of Knowledge graph
augmentation (Paulheim, 2017). It deduces missing information from the
knowledge graph itself.

The datasets for LP are constructed by sampling from the original knowledge
graph. Then, the links removed can be used in validation set or the test set
(Bordes et al., 2013; Dettmers et al., 2018). The structure of such graphs
play a vital role for improving the results, multiple source entities making
learning effective and multiple destination entities making learning difficult
(Rossi et al., 2021).

The LP models assigns a score to the triplet corresponding to each possible
entity to fill the question mark (?). The triplets are then ranked by a
function and entity corresponding to the lowest rank is predicted. If the
predicted facts in the ranked predictions are already present in the Knowledge
graph, they may or may not be excluded while calculating the ranks called raw
and filtered rankings respectively (Bordes et al., 2013). For e.g. if the
training knowledge graph contains the fact that (Arjuna, isSonOf, Kunti), and
the test query is (?, isSonOf, Kunti). The target answer is (Yudhishtra,
isSonOf, Kunti) and the system ranks (Arjuna, isSonOf, Kunti) and then
(Yudhishtra, isSonOf, Kunti). The raw ranking of the triplet (Yudhishtra,
isSonOf, Kunti) will be two and filtered ranking will be one.

There are several tie breaking policies that are used by the ranking system.
Assigning the minimum or the  maximum rank, or a random or the average rank to
the targeted entity (Rossi et al., 2021).

The ranks obtained are used to compute metrics such as Mean Rank (average of
all the ranks), Mean Reciprocal Rank (average of the inverse of ranks), or
Hits@M (proportion of ranks $\leq$ M).

###  3.2 Triple Classification

Triple Classification is the problem of identifying whether a given triple is
correct. It aims to give a yes or no answer to questions such as is New Delhi
capital of India? which can be written in the form of a triple (New Delhi,
isCapitalOf, India) (Socher et al., 2013).

A scoring function is used to calculate score of a triple similar to the link
prediction. If the score is greater than a certain threshold, then it is
considered a fact else a wrong triple (Wang et al., 2017).

Both the classical methods, such as micro and macro averaging, and ranking
methods such as Mean Rank are used as evaluation metrics (Guo et al., 2016).

###  3.3 Recommender Systems

Recommender system (RS) assists the user in an environment where multiple
options are available by providing a certain ordering of choices that the
recommendation algorithm infers. This inference can be based on the similarity
of the choices and behaviour pattern of different users. This type of
recommendation methods falls into the domain of collaborative filtering
methods (Adomavicius and Tuzhilin, 2005).

The CF methods suffer from problems of Data sparsity and cold start. Data
sparsity arises from the fact that only a small proportion of items are rated
by the users and most options have only limited feedback from the users. Cold
start problem is the problem of having no historical data about the new users
and options. To deal with these problems different types of side information
about a user and item are utilized by the RS (Sun et al., 2019b).

KG is utilised for side information in CF. It acts as a heterogeneous graph
that represent entities as nodes and relation as edges. The KG connects
various entities via latent relationships and also provide explainability in
recommendations (Wang et al., 2018).

The KG embedding based methods for RS use two modules - Graph embedding and
Recommendation Module. The way that these modules are coupled lead to
categorization of embedding based methods in a). two stage learning methods,
b). joint learning method and c). multi task learning method(Guo et al.,
2020).

Two stage learning methods first uses graph embedding module to obtain the
embeddings using various KG algorithms and then use recommendation module to
infer. The advantages of this method lies in its simplicity and scalability
but since the two modules are loosely coupled the embeddings might not be
suitable for recommendation tasks.

Joint learning methods train both the modules in an end to end fashion. Thus,
recommendation module guides the training in graph embedding layer.

Multi task learning method train the recommendation module with the guidance
of KG related task such as KG completion. The primary intuition behind this
method is that the bipartite graph of user and item in recommendation task
share structures with the corresponding KG entities.

##  4 Summary

Knowledge graphs provide an effective way of presenting real-world
relationships. As a result Knowledge graphs have an inherent advantage w.r.t
serving the information need. KG in itself is a growing area of research. KG
embedding is a technique to represent all the components of the KG in vector
form. These vectors represent the latent properties of the components of the
graph. Various models for embedding methods are based on different
combinations of vector algebra which present an interesting area of research.
In this work, we have surveyed the embedding methods that started this active
area of research, state-of-art models and the new frontiers which are being
explored in the KG embedding.
The KG embedding methods progressed from translation-based models which are
based on vector addition. In this work, we have presented how translation-
based models improved over time to overcome shortcomings of the earlier
models. While translation-based models used vector addition, semantic models
can be clubbed together as multiplicative models. We have included the
transition from basic semantic models to the more advanced semantic models
which may be used to explain different types of real world relationships such
as symmetric, anti-symmetric,inverse or composition.
New research areas have broaden the scope from structural embedding to more
contextual embedding by encoding additional information in the learned
representation. The latest area of research in this field is enrichment based
embedding models. In this work, we have introduced those briefly.
Vector space representation has paved a way to use the information from
Knowledge graph directly into the real world application. In this work, we
have described a few real world applications of KG embedding such as link
prediction, triple classification and recommender systems.

## References

  * Adomavicius and Tuzhilin (2005) Gediminas Adomavicius and Alexander Tuzhilin. 2005.  Toward the next generation of recommender systems: A survey of the state-of-the-art and possible extensions.  _IEEE transactions on knowledge and data engineering_ , 17(6):734–749.
  * Berant et al. (2013) Jonathan Berant, Andrew Chou, Roy Frostig, and Percy Liang. 2013.  Semantic Parsing on Freebase from Question-Answer Pairs.  Technical report.
  * Berners-Lee et al. (2001) Tim Berners-Lee, James Hendler, and Ora Lassila. 2001.  The Semantic Web A new form of Web content that is meaningful to computers will unleash a revolution of new possibilities.  _Sci. Am._ , 284(5):1–5.
  * Bianchi et al. (2020) F. Bianchi, G. Rossiello, Luca Costabello, M. Palmonari, and Pasquale Minervini. 2020.  Knowledge graph embeddings and explainable ai.  _ArXiv_ , abs/2004.14843.
  * Bollacker et al. (2008) Kurt Bollacker, Colin Evans, Praveen Paritosh, Tim Sturge, and Jamie Taylor. 2008\.  Freebase: A collaboratively created graph database for structuring human knowledge.  In _Proc. ACM SIGMOD Int. Conf. Manag. Data_ , pages 1247–1249, New York, New York, USA. ACM Press.
  * Bonatti et al. (2018) P. Bonatti, S. Decker, A. Polleres, and V. Presutti. 2018.  Knowledge graphs: New directions for knowledge representation on the semantic web (dagstuhl seminar 18371).  _Dagstuhl Reports_ , 8:29–111.
  * Bordes et al. (2014) Antoine Bordes, Xavier Glorot, Jason Weston, and Yoshua Bengio. 2014.  A semantic matching energy function for learning with multi-relational data: Application to word-sense disambiguation.  _Mach. Learn._ , 94(2):233–259.
  * (8) Antoine Bordes, Nicolas Usunier, Alberto Garcia-Durán, Jason Weston, and Oksana Yakhnenko.  Translating Embeddings for Modeling Multi-relational Data.  Technical report.
  * Bordes et al. (2013) Antoine Bordes, Nicolas Usunier, Alberto Garcia-Duran, Jason Weston, and Oksana Yakhnenko. 2013.  Translating embeddings for modeling multi-relational data.  In _Neural Information Processing Systems (NIPS)_ , pages 1–9.
  * Bouchard et al. (2015) Guillaume Bouchard, Sameer Singh, and Théo Trouillon. 2015.  On approximate reasoning capabilities of low-rank vector spaces.
  * Dettmers et al. (2018) Tim Dettmers, Pasquale Minervini, Pontus Stenetorp, and Sebastian Riedel. 2018.  Convolutional 2d knowledge graph embeddings.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 32.
  * Dietz et al. (2019) Laura Dietz, Chenyan Xiong, Jeff Dalton, and Edgar Meij. 2019.  Special issue on knowledge graphs and semantics in text analysis and retrieval.  _Inf. Retr. J._ , 22(3-4):229–231.
  * García-Durán et al. (2014) Alberto García-Durán, Antoine Bordes, and Nicolas Usunier. 2014.  Effective blending of two and three-way interactions for modeling multi-relational data.  In _Machine Learning and Knowledge Discovery in Databases_ , pages 434–449, Berlin, Heidelberg. Springer Berlin Heidelberg.
  * Guo et al. (2020) Q. Guo, F. Zhuang, C. Qin, H. Zhu, X. Xie, H. Xiong, and Q. He. 2020\.  A survey on knowledge graph-based recommender systems.  _IEEE Transactions on Knowledge and Data Engineering_ , pages 1–1.
  * Guo et al. (2016) Shu Guo, Quan Wang, Lihong Wang, Bin Wang, and Li Guo. 2016.  Jointly embedding knowledge graphs and logical rules.  In _Proceedings of the 2016 conference on empirical methods in natural language processing_ , pages 192–202.
  * Heck and Huang (2014) Larry Heck and Hongzhao Huang. 2014.  Deep learning of knowledge graph embeddings for semantic parsing of Twitter dialogs.  In _2014 IEEE Glob. Conf. Signal Inf. Process. Glob. 2014_ , pages 597–601. Institute of Electrical and Electronics Engineers Inc.
  * Ji et al. (2015) Guoliang Ji, Shizhu He, Liheng Xu, Kang Liu, and Jun Zhao. 2015.  Knowledge graph embedding via dynamic mapping matrix.  In _Proceedings of the 53rd Annual Meeting of the Association for Computational Linguistics and the 7th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , pages 687–696, Beijing, China. Association for Computational Linguistics.
  * Ji et al. (2020) Shaoxiong Ji, Shirui Pan, Erik Cambria, Pekka Marttinen, and Philip S. Yu. 2020\.  A survey on knowledge graphs: Representation, acquisition and applications.  _arXiv_ , pages 1–27.
  * Kolda and Bader (2009) Tamara G Kolda and Brett W Bader. 2009.  Tensor decompositions and applications.  _SIAM review_ , 51(3):455–500.
  * Lehmann et al. (2012) Jens Lehmann, Robert Isele, Max Jakob, Anja Jentzsch, Dimitris Kontokostas, Pablo N Mendes, Sebastian Hellmann, Mohamed Morsey, Patrick Van Kleef, Sören Auer, and Christian Bizer. 2012.  DBpedia-A Large-scale, Multilingual Knowledge Base Extracted from Wikipedia.  Technical report.
  * (21) Brian Lin, Shane Ong, Cory Williams, Yun Bin, and Matteo Zhang.  Named Entity Disambiguation with Knowledge Graphs.  Technical report.
  * Lin et al. (2015) Yankai Lin, Zhiyuan Liu, Maosong Sun, Yang Liu, and Xuan Zhu. 2015.  Learning entity and relation embeddings for knowledge graph completion.  In _Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence_ , AAAI’15, page 2181–2187. AAAI Press.
  * Liu et al. (2017) Hanxiao Liu, Yuexin Wu, and Yiming Yang. 2017.  Analogical inference for multi-relational embeddings.
  * Liu et al. (2018) Zhenghao Liu, Chenyan Xiong, Maosong Sun, and Zhiyuan Liu. 2018.  Entity-duet neural ranking: Understanding the role of knowledge graph semantics in neural information retrieval.  _ACL 2018 - 56th Annu. Meet. Assoc. Comput. Linguist. Proc. Conf. (Long Pap._ , 1:2395–2405.
  * Lu et al. (2020) Fengyuan Lu, Peijin Cong, and Xinli Huang. 2020.  Utilizing Textual Information in Knowledge Graph Embedding: A Survey of Methods and Applications.  _IEEE Access_ , 8:92072–92088.
  * Nickel et al. (2016a) Maximilian Nickel, Kevin Murphy, Volker Tresp, and Evgeniy Gabrilovich. 2016a.  A review of relational machine learning for knowledge graphs.  _Proceedings of the IEEE_ , 104(1):11–33.
  * Nickel et al. (2016b) Maximilian Nickel, Lorenzo Rosasco, and Tomaso Poggio. 2016b.  Holographic embeddings of knowledge graphs.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 30.
  * Nickel et al. (2011) Maximilian Nickel, Volker Tresp, and Hans-Peter Kriegel. 2011.  A three-way model for collective learning on multi-relational data.  In _Icml_.
  * Paulheim (2017) Heiko Paulheim. 2017.  Knowledge graph refinement: A survey of approaches and evaluation methods.  _Semantic web_ , 8(3):489–508.
  * Rossi et al. (2021) Andrea Rossi, Denilson Barbosa, Donatella Firmani, Antonio Matinata, and Paolo Merialdo. 2021.  Knowledge graph embedding for link prediction.  _ACM Transactions on Knowledge Discovery from Data_ , 15(2):1–49.
  * (31) Apoorv Saxena, Aditay Tripathi, and Partha Talukdar.  Improving Multi-hop Question Answering over Knowledge Graphs using Knowledge Base Embeddings.  Technical report.
  * Socher et al. (2013) Richard Socher, Danqi Chen, Christopher D Manning, and Andrew Ng. 2013.  Reasoning with neural tensor networks for knowledge base completion.  In _Advances in neural information processing systems_ , pages 926–934. Citeseer.
  * Suchanek et al. (2007) Fabian M. Suchanek, Gjergji Kasneci, and Gerhard Weikum. 2007.  Yago: A core of semantic knowledge.  In _16th Int. World Wide Web Conf. WWW2007_ , pages 697–706.
  * Sun et al. (2020) Rui Sun, Xuezhi Cao, Yan Zhao, Junchen Wan, Kun Zhou, Fuzheng Zhang, Zhongyuan Wang, and Kai Zheng. 2020.  Multi-modal Knowl-edge Graphs for Recommender Systems.
  * Sun et al. (2019a) Zhiqing Sun, Zhi-Hong Deng, Jian-Yun Nie, and Jian Tang. 2019a.  Rotate: Knowledge graph embedding by relational rotation in complex space.  _CoRR_ , abs/1902.10197.
  * Sun et al. (2019b) Zhu Sun, Qing Guo, Jie Yang, Hui Fang, Guibing Guo, Jie Zhang, and Robin Burke. 2019b.  Research commentary on recommendations with side information: A survey and research directions.  _Electronic Commerce Research and Applications_ , 37:100879.
  * Toutanova et al. (2016) Kristina Toutanova, Xi Victoria Lin, Wen Tau Yih, Hoifung Poon, and Chris Quirk. 2016.  Compositional learning of embeddings for relation paths in knowledge bases and text.  _54th Annu. Meet. Assoc. Comput. Linguist. ACL 2016 - Long Pap._ , 3:1434–1444.
  * Trouillon et al. (2016) Théo Trouillon, Johannes Welbl, Sebastian Riedel, Éric Gaussier, and Guillaume Bouchard. 2016.  Complex embeddings for simple link prediction.
  * Velicković et al. (2017) Petar Velicković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio. 2017.  Graph attention networks.  _arXiv_ , pages 1–12.
  * Vrandečić and Krötzsch (2014) Denny Vrandečić and Markus Krötzsch. 2014.  Wikidata: A free collaborative knowledgebase.  _Commun. ACM_ , 57(10):78–85.
  * Wang et al. (2018) Hongwei Wang, Fuzheng Zhang, Jialin Wang, Miao Zhao, Wenjie Li, Xing Xie, and Minyi Guo. 2018.  Ripplenet.  _Proceedings of the 27th ACM International Conference on Information and Knowledge Management_.
  * Wang et al. (2017) Quan Wang, Zhendong Mao, Bin Wang, and Li Guo. 2017.  Knowledge graph embedding: A survey of approaches and applications.  _IEEE Trans. Knowl. Data Eng._ , 29(12):2724–2743.
  * Wang et al. (2020) Shirui Wang, Wenan Zhou, and Chao Jiang. 2020.  A survey of word embeddings based on deep learning.  _Computing_ , 102(3):717–740.
  * (44) Xiang Wang, Xiangnan He, Yixin Cao, Meng Liu, and Tat-Seng Chua.  KGAT: Knowledge Graph Attention Network for Recommendation.
  * Wang et al. (2014) Zhen Wang, J. Zhang, Jianlin Feng, and Z. Chen. 2014.  Knowledge graph embedding by translating on hyperplanes.  In _AAAI_.
  * Wei et al. (2019) Yinwei Wei, Xiangnan He, Xiang Wang, Richang Hong, Liqiang Nie, and Tat Seng Chua. 2019.  MMGCN: Multi-modal graph convolution network for personalized recommendation of micro-video.  _MM 2019 - Proc. 27th ACM Int. Conf. Multimed._ , pages 1437–1445.
  * (47) Chenyan Xiong and Jamie Callan.  Query Expansion with Freebase.
  * Yang et al. (2014) Bishan Yang, Wen tau Yih, Xiaodong He, Jianfeng Gao, and Li Deng. 2014.  Embedding entities and relations for learning and inference in knowledge bases.
  * Zhang et al. (2019) Zhanqiu Zhang, Jianyu Cai, Yongdong Zhang, and Jie Wang. 2019.  Learning hierarchy-aware knowledge graph embeddings for link prediction.  _CoRR_ , abs/1911.09419.