# Learning Distances from Data with Normalizing Flows and Score Matching

Peter Sorrenson
Heidelberg University
&Daniel Behrend-Uriarte
Heidelberg University
Christoph Schnörr
Heidelberg University
&Ullrich Köthe
Heidelberg University
Correspondence to<EMAIL_ADDRESS>
###### Abstract

Density-based distances (DBDs) offer an elegant solution to the problem of
metric learning. By defining a Riemannian metric which increases with
decreasing probability density, shortest paths naturally follow the data
manifold and points are clustered according to the modes of the data. We show
that existing methods to estimate Fermat distances, a particular choice of
DBD, suffer from poor convergence in both low and high dimensions due to i)
inaccurate density estimates and ii) reliance on graph-based paths which are
increasingly rough in high dimensions. To address these issues, we propose
learning the densities using a normalizing flow, a generative model with
tractable density estimation, and employing a smooth relaxation method using a
score model initialized from a graph-based proposal. Additionally, we
introduce a dimension-adapted Fermat distance that exhibits more intuitive
behavior when scaled to high dimensions and offers better numerical
properties. Our work paves the way for practical use of density-based
distances, especially in high-dimensional spaces.

##  1 Introduction

Metric learning aims to determine a distance metric that accurately measures
the similarity or dissimilarity between data points. Typically, data points
are mapped into a space where simple metrics like Euclidean distance can be
applied. Although this method is computationally efficient, it is limited by
the configurations possible in Euclidean space, where arbitrary distance
relationships between sets of points cannot be represented.

A more flexible but computationally expensive approach involves defining a
Riemannian metric in the data space and solving for geodesic distances. This
allows for arbitrary distance relationships, but it requires choosing an
appropriate Riemannian metric. Fermat distances (Groisman et al., 2022), which
are a type of density-based distance (Bousquet et al., 2003), offer an elegant
solution: the metric should be conformal (a multiple of the identity matrix)
and inversely proportional to the probability density. In this way, geodesics
pass through high-density regions of the data, respecting any inherent
manifold structure. This approach is consistent with Fermat’s principle of
least time in optics, where the inverse density functions like a refractive
index.

To address the computational challenge of computing geodesics with this
metric, previous work has approximated trajectories using shortest paths in a
graph of nearest neighbors (Bijral et al., 2012). The edge weights are
approximated by raising Euclidean distances to an appropriate power, with
theoretical analysis supporting the consistency of this approximation.
However, in practice, we observe that convergence quickly plateaus, showing no
significant improvement with increasing sample size.

By employing more accurate edge weights obtained through normalizing flows, we
achieve much faster convergence. Additionally, previous approaches did not
account for the limitations of graph-based methods in higher dimensions, where
paths become increasingly rough due to the curse of dimensionality. To address
this, we use a relaxation scheme to smooth trajectories and maintain
satisfactory convergence. We find that training models via score matching is
more effective for this purpose than using normalizing flows.

In summary, our contributions are threefold: First, we introduce the use of
normalizing flows to obtain more accurate edge weights for estimating Fermat
distances, significantly improving convergence rates. Second, we address the
limitations of graph-based methods in high-dimensional spaces by implementing
a smooth relaxation scheme using score matching, which better maintains the
continuity and smoothness of the trajectories. Finally, we propose a
dimension-adapted Fermat distance that scales intuitively with increasing
dimensionality and exhibits superior numerical properties, paving the way for
practical applications of density-based distances in complex, high-dimensional
datasets. Through these innovations, we demonstrate the practicality and
effectiveness of our approach, opening up a new possibility for metric
learning using density-based distances in high-dimensional spaces.

Figure 1: Ground truth and estimated geodesics on a 2-dimensional spiral
dataset. The green line shows a ground truth geodesic, solved via relaxation
on the true density (Algorithm 1). The dashed orange line is an estimate for
the geodesic, found as the shortest path in a nearest neighbor graph with edge
weights calculated from a normalizing flow density (Algorithm 2).

##  2 Related work

Bousquet et al. (2003) introduces the idea of changing the geometry of the
data space by use of a conformal density-based metric in order to get more
informative distance functions. However they do not attempt to calculate
geodesics or distances implied by this geometry, instead linking it to a
modified SVM algorithm used for semi-supervised learning. Sajama and Orlitsky
(2005) also uses a density-based metric for semi-supervised learning. They
estimate densities using kernel density estimation and construct a graph of
nearest neighbor distances from which they compute geodesics.

Bijral et al. (2012) continues this line of work by approximating the density
as inversely proportional to a power of the Euclidean distance between close
neighbors. Hwang et al. (2016) gives convergence guarantees for the power-
weighted shortest path distances to geodesic distances under the density-
dependent metric. Several other papers give similar results, and some
practical implementations (Chu et al., 2020; Little et al., 2022; Groisman et
al., 2022; Moscovich et al., 2017; Alamgir and Von Luxburg, 2012; Mckenzie and
Damelin, 2019; Little et al., 2020). Terms used are power-weighted shortest
path distances (PWSPD) (Little et al., 2022) and Fermat distances (Groisman et
al., 2022; Trillos et al., 2023).

One crucial contribution of our work is to check that proposed methods
actually converge in practice, regardless of theoretical guarantees. Previous
work does not compare shortest paths to ground truth geodesics, even for a
uniform distribution, where geodesics are simply straight lines. We address
this major shortcoming by computing the ground truth geodesics for problems
with known density and find that existing methods show very poor convergence.

##  3 Background

###  3.1 Riemannian geometry

Riemannian geometry offers powerful tools for analyzing curved spaces, which
can be particularly useful for understanding and optimizing over non-Euclidean
spaces such as manifolds. A fundamental object in this field is the metric
tensor, where different choices of metric can lead to vastly different
geometric structures, such as positive curvature (e.g., a sphere which closes
in on itself) or negative curvature (e.g., hyperbolic space which expands
exponentially quickly), as well as more complicated mixtures of the two.

More formally, a Riemannian manifold $\mathcal{M}$ is a smooth $n$-dimensional
space that locally resembles $\mathbb{R}^{n}$. At each point $p$ on the
manifold, there is an associated tangent space $\mathcal{T}_{p}$, which is a
linear approximation of the manifold at $p$. The metric tensor $g$ is a
smoothly varying positive-definite bilinear form on the tangent space. It
defines an inner product in the tangent space, allowing us to measure lengths
and angles.

For a tangent vector $v\in\mathcal{T}_{p}$, the norm is defined as
$\|v\|=\sqrt{g(v,v)}$, providing a notion of length in the tangent space. The
length of a smooth curve $\gamma(t)$ on the manifold is given by
$\int\sqrt{g(\dot{\gamma},\dot{\gamma})}\,dt$, where $\dot{\gamma}$ is the
derivative of $\gamma(t)$ with respect to $t$. This integral measures the
total length of the curve by summing up the infinitesimal lengths along the
curve.

The distance between two points on the manifold is defined as the minimal
length of a curve connecting the points. The curves that minimize this length
are known as geodesics, which are the generalization of straight lines in
Euclidean space to curved spaces.

###  3.2 Fermat distances

Density-based distances (DBDs) provide a way to define distances between
points in a distribution by utilizing the underlying probability density
function. This approach can be particularly useful for tasks such as
clustering and pathfinding in data analysis and machine learning. If designed
well, a DBD can capture desirable properties such as:

  1. 1.

Shortest paths between points should pass through supported regions of the
data. If there is a region with no data between two points, the shortest path
should deviate around it to pass through regions with more data.

  2. 2.

Data should cluster according to its modes. If two data points belong to
different modes with little data connecting them, the distance between them
should be high.

One way to define a DBD that satisfies these properties is through Fermat
distances. Named after Fermat’s principle of least time in optics, these
distances model the shortest paths between points by considering the inverse
probability density (or a monotonic function of it) as analogous to the
refractive index.

In the language of Riemannian geometry, we define a metric tensor of the form

$g(u,v)=\frac{\langle u,v\rangle}{f(p)^{2}},$ (1)

where $p$ is the probability density, $f$ is a monotonic function, and
$\langle\cdot,\cdot\rangle$ is the Euclidean inner product. This formulation
ensures that regions with higher probability density have lower "refractive
index," guiding shortest paths through denser regions of the data. In
practice, it is convenient to choose $f(p)=p^{\beta}$, where $\beta$ is an
inverse temperature parameter that controls the influence of the density on
the metric. Since our metric enacts a simple scaling of the Euclidean inner
product, it is a conformal metric, which locally preserves angles and
simplifies key geometric objects such as the geodesic equation and scalar
curvature.

The length of a smooth curve $\gamma(t)$ is therefore given by

$L(\gamma)=\int\frac{\|\dot{\gamma}(t)\|}{p(\gamma(t))^{\beta}}\,dt$ (2)

####  3.2.1 Fermat distances which scale with dimension

Figure 2: Geodesics of the standard normal distribution, projected into 2d,
taken between two orthogonal points at distance $\sqrt{D}$ from the origin.
The paths are scaled by $\sqrt{D}$ to enable comparison. (Left) Unscaled
metric (same temperature for all dimensions) leads to sharply curved geodesics
in high dimensions. (Right) Scaled metric (temperature equal to dimension)
leads to consistent geodesics in all dimensions (note that the lines are
overlapping).

If $\beta=1$, this leads to unintuitive behavior in higher dimensions.
Consider the standard normal distribution in $D$ dimensions. Due to rotational
symmetry, the geodesic connecting $x_{1}$ and $x_{2}$ will lie in the plane
spanned by these two points and the origin. Without loss of generality,
suppose this plane aligns with the first two axis directions. The density
restricted to this plane is the same as in the two-dimensional case, up to a
constant factor. Therefore, geodesics in the $D$-dimensional standard normal
distribution are the same as in the two-dimensional standard normal
distribution, once we account for the rotation into the plane of $x_{1}$ and
$x_{2}$.

However, despite this similarity, the typical distance from the origin of
points sampled from the $D$-dimensional standard normal distribution is much
larger, with a mean of $\sqrt{D}$. As a result, geodesics will start and end
in very low-density regions relative to the two-dimensional standard normal.
This means that most of the trajectory is spent moving towards and then
outwards from the high-density region near the origin.

To address this issue, if we scale the temperature linearly with dimension
($\beta=1/D$), the results become consistent with the two-dimensional case
(see Fig. 2). Scaling the temperature with dimension also increases numerical
stability. This is because the effect of the probability density is reduced,
leading to more Euclidean-like geodesics. This effect is clearly visible in
Fig. 2, where the geodesics are much closer to straight lines in the scaled
case.

###  3.3 Geodesic equation

From the metric, we can compute the geodesic equation. This is the second-
order differential equation which all paths of minimum length (geodesics) must
satisfy. In our case, we have the relatively simple solution for a trajectory
$\gamma(t)$:

$\ddot{\gamma}-2\beta(s(\gamma)\cdot\dot{\gamma})\dot{\gamma}+\beta
s(\gamma)\|\dot{\gamma}\|^{2}=0$ (3)

where $s(x)=\frac{\partial\log p(x)}{\partial x}$ is the score of the
distribution and $\|\cdot\|$ is the ordinary Euclidean norm.

The geodesic equation has the property that the speed is constant with respect
to the metric, meaning $\sqrt{g(\dot{\gamma},\dot{\gamma})}$ is constant. In
other words, the geodesic distance traveled in a given time interval is the
same at all points of the curve. If $\gamma$ passes from a high-probability
region to a low-probability one, the distances in the low-probability region
could be orders of magnitude larger than in the high-probability region. This
makes numerically solving the equations challenging due to the vastly
different Euclidean velocities along the trajectory.

It is convenient to reparameterize the equation such that the Euclidean
velocity is constant instead, as this is much easier to handle numerically.
This is especially useful for solving via relaxation, where having equal-
length segments leads to much more stable behavior. Making this change results
in a very similar equation (removing only the factor of 2 in the second term):

$\ddot{\varphi}-\beta(s(\varphi)\cdot\dot{\varphi})\dot{\varphi}+\beta
s(\varphi)\|\dot{\varphi}\|^{2}=0$ (4)

See Section A.1 for derivations of these results.

##  4 Methods

###  4.1 Ground truth distances through relaxation

A relaxation scheme allows one to solve a differential equation constrained by
boundary conditions by iteratively updating the solution such that it better
obeys the differential equation. We divide the interval $[0,1]$ into equal
time steps of size $h=1/n$ and write $\varphi_{i}=\varphi(ih)$. Initialize the
values of the curve $\varphi$ at $t=ih$ for $i=0,\dots,n$, with the constraint
that $\varphi_{0}$ and $\varphi_{n}$ are the chosen endpoints. The initial
curve could be the straight line connecting the end points, or a more
informative guess if additional information is available.

Equation 4 leads to the following relaxation scheme. Define

$v_{i}=\frac{\varphi_{i+1}-\varphi_{i-1}}{2}$ (5)

then update with

$\varphi^{\prime}_{i}=\frac{\varphi_{i+1}+\varphi_{i-1}}{2}+\frac{\beta}{2}\left(s(\varphi_{i})\|v_{i}\|^{2}-(s(\varphi_{i})\cdot
v_{i})v_{i}\right)$ (6)

See Section A.2 for a derivation.

####  4.1.1 Implementation

In practice, some additional features are needed to ensure stable convergence
of the relaxation scheme, especially when the score function is not
particularly smooth. Firstly, we found that oscillatory features commonly
arise. To address this, we use smoothing after the update, returning a value
which is interpolated between the value at the node and the average of its
neighbors. Secondly, for rougher score functions we use a step size between 0
and 1, where 0 corresponds to no update and 1 corresponds to the ordinary
relaxation step. See Algorithm 1.

Algorithm 1 One step of relaxation. $\varphi=\\{\varphi_{i}\\}_{i=0}^{n}$ is
the current state of the trajectory, $s$ is the score function, $\beta$ the
inverse temperature, $\gamma\in[0,1)$ the smoothing parameter and
$\alpha\in(0,1]$ the step size. This step should be iterated until
convergence, which can be measured via the difference between $\varphi$ and
$\varphi^{\prime}$.

function RelaxationStep($\varphi,s,\beta,\gamma,\alpha$)

$v_{i}\leftarrow\frac{1}{2}(\varphi_{i+1}-\varphi_{i-1})$ $\triangleright$
$i=1,\dots,n-1$

$\varphi^{\prime}_{i}\leftarrow\frac{1}{2}(\varphi_{i+1}+\varphi_{i-1})+\frac{\beta}{2}\left(s(\varphi_{i})\|v_{i}\|^{2}-(s(\varphi_{i})\cdot
v_{i})v_{i}\right)$

$\varphi^{\prime}_{i}\leftarrow(1-\gamma)\varphi^{\prime}_{i}+\frac{\gamma}{2}(\varphi^{\prime}_{i+1}+\varphi^{\prime}_{i-1})$

$\varphi^{\prime}_{i}\leftarrow\alpha\varphi^{\prime}_{i}+(1-\alpha)\varphi_{i}$

return $\varphi^{\prime}$

end function

###  4.2 Weighted graph algorithm based on Euclidean distance

If we have a large number of samples from $p(x)$, we can approximate the
density between close neighbors $x_{1}$ and $x_{2}$ as constant and hence the
shortest path between them as a straight line. The density is roughly
proportional to the inverse of the Euclidean distance, raised to the power
$d$, where $d$ is the intrinsic dimension of the distribution (Bijral et al.,
2012):

$p(\gamma)\propto\frac{1}{\|x_{1}-x_{2}\|^{d}}$ (7)

Hence the distance between $x_{1}$ and $x_{2}$ according to the metric is
approximately proportional to a power of the Euclidean distance

$\operatorname{dist}(x_{1},x_{2})\propto\|x_{1}-x_{2}\|^{\beta d+1}$ (8)

Bijral et al. (2012) therefore proposes to approximate shortest paths between
more distant points by constructing a $k$-nearest neighbors graph of the
samples with edge weights given by $\|x_{1}-x_{2}\|^{\beta d+1}$ and applying
Dijkstra’s algorithm.

Hwang et al. (2016) gives consistency guarantees for this estimate, proving
that a scaled version of it converges to the ground truth distance with a rate
of $\exp(-\theta_{0}n^{1/(3d+2)})$, with $n$ the sample size, $d$ the
intrinsic dimension, and $\theta_{0}$ a positive constant. Groisman et al.
(2022) additionally show that the shortest paths converge to geodesics in a
large sample size limit.

While this approach is appealing in its simplicity, can be easily implemented
using standard libraries, and comes with theoretical guarantees, the density
estimate is too poor to provide accurate estimates of Fermat distances in
practice. In the experimental section 5.1.1 we show that paths computed with
this approach do not converge to the ground truth geodesics as sample size
increases, even in very simple 2-dimensional distributions such as the
standard normal.

We are unsure why we observe such a large deviation from the theoretical
prediction. It is possible that the unknown constant $\theta_{0}$ is very
small, leading to extremely slow convergence, despite the exponential
function. We leave investigation into this point to future work.

###  4.3 Normalizing flows

Normalizing flows (Rezende and Mohamed, 2015; Kobyzev et al., 2020;
Papamakarios et al., 2021) are a set of techniques to learn a probability
density via a parameterized diffeomorphism $f_{\theta}$ which maps from data
to latent space. The probability density is obtained through the change of
variables formula:

$p_{\theta}(x)=\pi(f_{\theta}(x))\left|\det\left(\frac{\partial
f_{\theta}}{\partial x}\right)\right|^{-1}$ (9)

where $\pi$ is a simple latent distribution. A common strategy to design
normalizing flows is through coupling blocks (Dinh et al., 2014, 2016; Durkan
et al., 2019), where $f_{\theta}$ is the composition of a series of blocks,
each with a triangular Jacobian. This makes the determinant easy to calculate,
and hence the above formula is tractable.

Normalizing flows are trained by minimizing the negative log-likelihood on a
training set:

$\mathcal{L}=\mathbb{E}_{p_{\text{data}}(x)}[-\log p_{\theta}(x)]$ (10)

###  4.4 Density-weighted graph algorithm

In order to better estimate probability densities, we use a normalizing flow
trained on the data and directly use the density predicted by the model. In
order to approximate the geodesic distance between close neighbors, we assume
that the geodesic connecting them is approximately a straight line, and
estimate its length by a discrete approximation:

$\operatorname{dist}(x_{1},x_{2})\approx\sum_{i=1}^{S}\frac{\|y_{i}-y_{i-1}\|}{p_{\theta}(y_{i-1/2})}$
(11)

Here the $y_{i}$ are points which interpolate $x_{1}$ and $x_{2}$ with
$y_{0}=x_{1}$ and $y_{S}=x_{2}$, and $y_{i-1/2}$ is the midpoint of $y_{i-1}$
and $y_{i}$.

This allows us to take a k-nearest neighbors graph of the data and replace the
edge weights by the approximate distance. This results in a sparse graph, to
which we can apply Dijkstra’s algorithm to estimate shortest paths and
distances between any two points. See Algorithm 2.

Algorithm 2 Construction of a density-weighted graph which can subsequently be
used to compute shortest paths and distances between points. $X$ is the data
matrix, $p$ the probability density, $\beta$ the inverse temperature, $k$ the
number of neighbors, and $S$ the number of segments with which to approximate
the integral. The computations can also be performed in the log space,
returning log edge weights, which is more numerically stable in higher
dimensions.

function DensityGraph($X,p,\beta,k,S$)

$\mathcal{G}\leftarrow\texttt{KNN}(X,k)$ $\triangleright$ Construct
$k$-nearest neighbor graph

for $(l,m)\in\mathcal{G}$ do$\triangleright$ Iterate over edges

$y_{i}\leftarrow(1-\frac{i}{S})X_{l}+\frac{i}{S}X_{m}$$\triangleright$ For
$i=0,\dots,S$

$y_{i-1/2}\leftarrow\frac{1}{2}(y_{i-1}+y_{i})$$\triangleright$ For
$i=1,\dots,S$

$\mathcal{G}_{lm}\leftarrow\sum_{i=1}^{S}\|y_{i}-y_{i-1}\|/p(y_{i-1/2})^{\beta}$

end for

return $\mathcal{G}$

end function

###  4.5 Relaxation with a learned score model

The score function is the gradient of the log-likelihood with respect to
spatial inputs, i.e., $s(x)=\nabla_{x}\log p(x)$. Note the difference from the
definition commonly used in statistic, where the derivative is with respect to
model parameters.

Score matching (Hyvärinen, 2005) is a family of methods to estimate this
function from data, including variants such as sliced score matching (Song et
al., 2020a), denoising score matching (Vincent, 2011) and diffusion models
(Song et al., 2020b).

##### Sliced score matching

We can learn the score of $p(x)$ by minimizing the following objective (known
as the score matching objective):

$\mathbb{E}_{p(x)}[\operatorname{tr}(\nabla_{x}s_{\theta}(x))+\tfrac{1}{2}\|s_{\theta}(x)\|^{2}]$
(12)

In high dimensions the trace term is expensive to evaluate exactly, so sliced
score matching uses a stochastic approximation (the Hutchinson trace estimator
(Hutchinson, 1989)):

$\mathbb{E}_{p(x)p(v)}[v^{T}\nabla_{x}s_{\theta}(x)v+\tfrac{1}{2}\|s_{\theta}(x)\|^{2}]$
(13)

where $v$ is sampled from an appropriate distribution, typically standard
normal.

Since only the score is needed for the relaxation algorithm described in
Section 4.1, we can use a trained score model to perform relaxation to
approximate geodesics. In practice, we initialize the path from a graph-based
method to speed up convergence and avoid local minima.

##  5 Experiments

##### Performance metric

After solving for a shortest path $\varphi$ between two points $x_{1}$ and
$x_{2}$, we compute its length $L(\varphi)$ by using the ground truth density
function, then compare this value to the ground truth distance between the two
points:

$M(\varphi)=\log\frac{L(\varphi)}{\operatorname{dist}(x_{1},x_{2})}$ (14)

Provided we have access to the ground truth density, this performance metric
can be calculated independently of how $\varphi$ is obtained, and has a value
of zero only if $\varphi$ is the ground truth geodesic. The metric will
otherwise take positive values.

###  5.1 Scaling the dataset

Figure 3: The length of shortest paths in power-weighted distance graphs do
not converge to the ground truth with increasing sample size (green). In
contrast, learned and ground-truth density weighted graphs (orange and green)
do converge, at similar rates.

In this section we perform experiments on a set of 5 two-dimensional datasets
for which we know the ground truth density (see Fig. 3). Details on how the
datasets were created can be found in Appendix B.

####  5.1.1 Power-weighted graph paths do not converge with increasing sample
size

We first estimate paths using the power-weighted graph method (Bijral et al.,
2012), applied to a k-nearest neighbor graph. We find that in all 5 datasets,
increasing the sample size leads to saturating performance (see Fig. 3, green
lines). In Fig. 6 in the appendix we show that increasing the number of
neighbors does not improve performance.

####  5.1.2 Poor convergence is due to poor density estimation

In order to investigate the root cause of the poor convergence, we designed
other graph methods which use the nearest neighbor estimator for the density,
but do not use power-weighted graph edges. We use 4 variants for the density
estimate (see Section C.1). We find that all 4 have very similar performance
to the power-weighted graph, whereas the equivalent estimators using the
ground truth density quickly converge. See Section C.1 for plots supporting
this claim.

####  5.1.3 Learning the density leads to convergent paths

Since the problem lies in poor density estimation, we hope to fix the problem
by better density estimation. To this end, we use a normalizing flow to learn
the density from samples, then use Algorithm 2 to construct a graph over which
we compute shortest paths. We find that the paths obtained converge at close
to the same rate as paths obtained by an equivalent method using the ground
truth densities (Fig. 3 orange and blue lines). Full experimental details can
be found in Appendix B.

###  5.2 Scaling the dimension

Figure 4: Mean log error with a sample size of 200,000. Graph-based paths (the
first three methods) converge very slowly in higher dimensions with
increasingly worse performance as dimension increases. When using the scaled
metric ($\beta=1/D$), convergence improves over the unscaled metric, but all
graph-based methods perform similarly in higher dimensions. Doing relaxation
with a learned score function (red line) leads to much faster convergence,
avoiding the curse of dimensionality.

Machine learning datasets often lie in high-dimensional spaces, so it is
crucial that any method applied to them scales well, avoiding the curse of
dimensionality as much as possible. We test this by taking one of our two-
dimensional datasets, namely the standard normal, and scaling it up to 25
dimensions. We perform experiments both for the scaled and unscaled version of
the metric.

####  5.2.1 Graph-based paths do not converge well in higher dimensions

Our first finding is that all graph-based methods show declining performance
as the dimension increases, for both unscaled and scaled metrics (Fig. 4).
This can be straightforwardly explained by the curse of dimensionality: with
increasing dimension but fixed sample size, it is increasingly unlikely that
points in the dataset lie close to the true trajectory. The shortest path in
the graph becomes increasingly rough in contrast to the smooth ground truth
trajectory.

####  5.2.2 Relaxation using a learned score function improves convergence
rates

In order to obtain smoother paths, we turn to the same relaxation algorithm we
use to solve for ground truth geodesics. Instead of using the ground truth
score function, we use a learned score function. An obvious choice would be to
use automatic differentiation to get the score from a trained normalizing
flow, but we find the resulting scores are too noisy to be useful, inhibiting
reliable convergence of the relaxation scheme. Instead, we learn a score model
with sliced score matching, which results in a smoother function suitable for
performing relaxation. Full experimental details can be found in Appendix B.

We find that this fixes the convergence problems, both for the unscaled and
scaled metrics, though convergence is more reliable in the scaled case (Fig.
4). Please see Fig. 8 and Fig. 9 in the appendix for more detailed evidence.

##  6 Conclusion

In conclusion, our study focuses on understanding the convergence properties
of various algorithms to ground-truth geodesics by employing toy distributions
with known densities. We investigate the standard normal distribution in
dimensions higher than 2, highlighting the challenges faced by graph-based
methods even in this simple scenario. Our findings demonstrate significant
improvements through the use of normalizing flows and score models, yet
further research is needed to address more complex distributions and real-
world data where the underlying density is unknown. Future work should aim to
unify the roles of normalizing flows and score models into a single cohesive
model, enhancing the method’s applicability and efficiency. Additionally, a
deeper theoretical analysis is required to understand why our proposed methods
converge effectively, while previous methods with consistency guarantees do
not. This work lays the groundwork for practical applications of density-based
distances, setting the stage for more robust and scalable solutions in metric
learning.

## References

  * Alamgir and Von Luxburg (2012) Morteza Alamgir and Ulrike Von Luxburg.  Shortest path distance in random k-nearest neighbor graphs.  _arXiv preprint arXiv:1206.6381_ , 2012.
  * Ardizzone et al. (2018-2022) Lynton Ardizzone, Till Bungert, Felix Draxler, Ullrich Köthe, Jakob Kruse, Robert Schmier, and Peter Sorrenson.  Framework for Easily Invertible Architectures (FrEIA), 2018-2022.  URL https://github.com/vislearn/FrEIA.
  * Bijral et al. (2012) Avleen S Bijral, Nathan Ratliff, and Nathan Srebro.  Semi-supervised learning with density based distances.  _arXiv preprint arXiv:1202.3702_ , 2012.
  * Bousquet et al. (2003) Olivier Bousquet, Olivier Chapelle, and Matthias Hein.  Measure based regularization.  _Advances in Neural Information Processing Systems_ , 16, 2003.
  * Chu et al. (2020) Timothy Chu, Gary L Miller, and Donald R Sheehy.  Exact computation of a manifold metric, via Lipschitz embeddings and shortest paths on a graph.  In _Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms_ , pages 411–425. SIAM, 2020.
  * Dinh et al. (2014) Laurent Dinh, David Krueger, and Yoshua Bengio.  Nice: Non-linear independent components estimation.  _arXiv preprint arXiv:1410.8516_ , 2014.
  * Dinh et al. (2016) Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio.  Density estimation using Real NVP.  _arXiv preprint arXiv:1605.08803_ , 2016.
  * Durkan et al. (2019) Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios.  Neural spline flows.  _Advances in neural information processing systems_ , 32, 2019.
  * Groisman et al. (2022) Pablo Groisman, Matthieu Jonckheere, and Facundo Sapienza.  Nonhomogeneous Euclidean first-passage percolation and distance learning.  _Bernoulli_ , 28(1):255–276, 2022.
  * Harris et al. (2020) Charles R Harris, K Jarrod Millman, Stéfan J Van Der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J Smith, et al.  Array programming with NumPy.  _Nature_ , 585(7825):357–362, 2020.
  * Hunter (2007) John D Hunter.  Matplotlib: A 2d graphics environment.  _Computing in science & engineering_, 9(03):90–95, 2007.
  * Hutchinson (1989) Michael F Hutchinson.  A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines.  _Communications in Statistics-Simulation and Computation_ , 18(3):1059–1076, 1989.
  * Hwang et al. (2016) Sung Jin Hwang, Steven B Damelin, and Alfred O Hero III.  Shortest path through random points.  2016\.
  * Hyvärinen (2005) Aapo Hyvärinen.  Estimation of non-normalized statistical models by score matching.  _Journal of Machine Learning Research_ , 6(4), 2005.
  * Kobyzev et al. (2020) Ivan Kobyzev, Simon JD Prince, and Marcus A Brubaker.  Normalizing flows: An introduction and review of current methods.  _IEEE transactions on pattern analysis and machine intelligence_ , 43(11):3964–3979, 2020.
  * Little et al. (2020) Anna Little, Mauro Maggioni, and James M Murphy.  Path-based spectral clustering: Guarantees, robustness to outliers, and fast algorithms.  _Journal of machine learning research_ , 21(6):1–66, 2020.
  * Little et al. (2022) Anna Little, Daniel McKenzie, and James M Murphy.  Balancing geometry and density: Path distances on high-dimensional data.  _SIAM Journal on Mathematics of Data Science_ , 4(1):72–99, 2022.
  * Mckenzie and Damelin (2019) Daniel Mckenzie and Steven Damelin.  Power weighted shortest paths for clustering Euclidean data.  _arXiv preprint arXiv:1905.13345_ , 2019.
  * Moscovich et al. (2017) Amit Moscovich, Ariel Jaffe, and Nadler Boaz.  Minimax-optimal semi-supervised regression on unknown manifolds.  In _Artificial Intelligence and Statistics_ , pages 933–942. PMLR, 2017.
  * Papamakarios et al. (2021) George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan.  Normalizing flows for probabilistic modeling and inference.  _Journal of Machine Learning Research_ , 22(57):1–64, 2021.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32, 2019.
  * Pedregosa et al. (2011) Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  Scikit-learn: Machine learning in Python.  _the Journal of machine Learning research_ , 12:2825–2830, 2011.
  * Rezende and Mohamed (2015) Danilo Rezende and Shakir Mohamed.  Variational inference with normalizing flows.  In _International conference on machine learning_ , pages 1530–1538. PMLR, 2015.
  * Sajama and Orlitsky (2005) Sajama and Alon Orlitsky.  Estimating and computing density based distance metrics.  In _Proceedings of the 22nd international conference on Machine learning_ , pages 760–767, 2005.
  * Song et al. (2020a) Yang Song, Sahaj Garg, Jiaxin Shi, and Stefano Ermon.  Sliced score matching: A scalable approach to density and score estimation.  In _Uncertainty in Artificial Intelligence_ , pages 574–584. PMLR, 2020a.
  * Song et al. (2020b) Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  _arXiv preprint arXiv:2011.13456_ , 2020b.
  * Trillos et al. (2023) Nicolás García Trillos, Anna Little, Daniel McKenzie, and James M Murphy.  Fermat distances: Metric approximation, spectral convergence, and clustering algorithms.  _arXiv preprint arXiv:2307.05750_ , 2023.
  * Vincent (2011) Pascal Vincent.  A connection between score matching and denoising autoencoders.  _Neural computation_ , 23(7):1661–1674, 2011.
  * Virtanen et al. (2020) Pauli Virtanen, Ralf Gommers, Travis E Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, et al.  SciPy 1.0: fundamental algorithms for scientific computing in Python.  _Nature methods_ , 17(3):261–272, 2020.

##  Appendix A Derivations

###  A.1 Geodesic equation

The geodesic equation which characterizes paths of zero acceleration $\gamma$
is

$\ddot{\gamma}^{i}+\Gamma^{i}_{jk}\dot{\gamma}^{j}\dot{\gamma}^{k}=0$ (15)

where $\Gamma^{i}_{jk}$ are the Christoffel symbols of the second kind

$\Gamma^{i}_{jk}=\frac{1}{2}g^{il}\left(\frac{\partial g_{lj}}{\partial
x^{k}}+\frac{\partial g_{lk}}{\partial x^{j}}-\frac{\partial g_{jk}}{\partial
x^{l}}\right)$ (16)

and $g^{il}$ is the inverse of $g$, i.e. such that
$g^{il}g_{lj}=\delta^{i}_{j}$.

Our metric has the form $g_{ij}=\lambda^{2}\delta_{ij}$, meaning that

$\displaystyle\Gamma^{i}_{jk}$
$\displaystyle=\frac{1}{2\lambda^{2}}\delta^{il}\left(\frac{\partial\lambda^{2}}{\partial
x^{k}}\delta_{lj}+\frac{\partial\lambda^{2}}{\partial
x^{j}}\delta_{lk}-\frac{\partial\lambda^{2}}{\partial
x^{l}}\delta_{jk}\right)$ (17)
$\displaystyle=\delta^{il}\left(\frac{\partial\log\lambda}{\partial
x^{k}}\delta_{lj}+\frac{\partial\log\lambda}{\partial
x^{j}}\delta_{lk}-\frac{\partial\log\lambda}{\partial
x^{l}}\delta_{jk}\right)$ (18)
$\displaystyle=\frac{\partial\log\lambda}{\partial
x^{k}}\delta^{i}_{j}+\frac{\partial\log\lambda}{\partial
x^{j}}\delta^{i}_{k}-\frac{\partial\log\lambda}{\partial
x^{l}}\delta^{il}\delta_{jk}$ (19)

In this case $\lambda=p^{-\beta}$ and hence
$\frac{\partial\log\lambda}{\partial x}=-\beta\frac{\partial\log p}{\partial
x}=-\beta s$, where $s=\frac{\partial\log p}{\partial x}$ is the score
function of the probability density. As a result

$\Gamma^{i}_{jk}=-\beta s_{k}\delta^{i}_{j}-\beta s_{j}\delta^{i}_{k}+\beta
s_{i}\delta^{il}\delta_{jk}$ (20)

Substituting into the geodesic equation leads to

$\ddot{\gamma}^{i}-\beta s_{k}\dot{\gamma}^{i}\dot{\gamma}^{k}-\beta
s_{j}\dot{\gamma}^{j}\dot{\gamma}^{i}+\beta
s_{l}\delta^{il}\delta_{jk}\dot{\gamma}^{j}\dot{\gamma}^{j}=0$ (21)

which can be written in vector notation as

$\ddot{\gamma}-2\beta(s(\gamma)\cdot\dot{\gamma})\dot{\gamma}+\beta
s(\gamma)\|\dot{\gamma}\|^{2}=0$ (22)

Let $f:[0,1]\rightarrow[0,1]$ be a strictly increasing diffeomorphism. Given a
curve $\gamma:[0,1]\rightarrow\mathbb{R}^{n}$, consider the curve
$\varphi:[0,1]\rightarrow\mathbb{R}^{n}$ defined such that
$\varphi(u)=\gamma(f^{-1}(u))$ or equivalently $\varphi(f(t))=\gamma(t)$. The
derivatives of the two curves are related by

$\dot{\varphi}\dot{f}=\dot{\gamma}$ (23)

and

$\ddot{\varphi}\dot{f}^{2}+\dot{\varphi}\ddot{f}=\ddot{\gamma}$ (24)

where derivatives of $\varphi$ are with respect to $u=f(t)$. By substituting
these expressions into the geodesic equation, we arrive at the equivalent
equation for $\varphi$:

$\ddot{\varphi}\dot{f}^{2}+\dot{\varphi}\ddot{f}-2\beta(s(\gamma)\cdot\dot{\varphi}\dot{f})\dot{\varphi}\dot{f}+\beta
s(\gamma)\|\dot{\varphi}\dot{f}\|^{2}=0$ (25)

We would like $\varphi$ to have a constant Euclidean speed (the magnitude but
not necessarily the direction of the velocity is constant), meaning

$\frac{d}{du}\left(\frac{1}{2}\|\dot{\varphi}\|^{2}\right)=\dot{\varphi}\cdot\ddot{\varphi}=0$
(26)

Multiplying the geodesic equation by $\dot{\varphi}$, this is equivalent to
requiring

$\|\dot{\varphi}\|^{2}\ddot{f}-2\beta(s(\gamma)\cdot\dot{\varphi})\|\dot{\varphi}\|^{2}\dot{f}^{2}+\beta(s(\gamma)\cdot\dot{\varphi})\|\dot{\varphi}\|^{2}\dot{f}^{2}=0$
(27)

implying that

$\ddot{f}=\beta(s(\gamma)\cdot\dot{\varphi})\dot{f}^{2}$ (28)

Resubstituting into the geodesic equation and dividing by $\dot{f}^{2}$ leads
to the result:

$\ddot{\varphi}-\beta(s(\varphi)\cdot\dot{\varphi})\dot{\varphi}+\beta
s(\varphi)\|\dot{\varphi}\|^{2}=0$ (29)

###  A.2 Relaxation scheme

The central finite difference approximations for the first and second
derivatives are

$\dot{\varphi}_{i}\approx\frac{\varphi_{i+1}-\varphi_{i-1}}{2h}$ (30)

and

$\ddot{\varphi}_{i}\approx\frac{\varphi_{i+1}-2\varphi_{i}+\varphi_{i-1}}{h^{2}}$
(31)

By substituting this into the differential equation we have

$\frac{\varphi_{i+1}-2\varphi_{i}+\varphi_{i-1}}{h^{2}}-\beta(s(\varphi_{i})\cdot\dot{\varphi}_{i})\dot{\varphi}_{i}+\beta
s(\varphi_{i})\|\dot{\varphi}_{i}\|^{2}\approx 0$ (32)

where $\dot{\varphi}_{i}$ is estimated using finite differences. We can
rearrange this to

$\varphi_{i}\approx\frac{\varphi_{i+1}+\varphi_{i-1}}{2}+\frac{h^{2}\beta}{2}\left(s(\varphi_{i})\|\dot{\varphi}_{i}\|^{2}-(s(\varphi_{i})\cdot\dot{\varphi}_{i})\dot{\varphi}_{i}\right)$
(33)

and use the equation as an update rule, updating each position of the curve
except the endpoints at each iteration.

Dividing by very small $h$ could lead to numerical instability. We can avoid
dividing and multiplying by $h$ by the following update rule. First define

$v_{i}=\frac{\varphi_{i+1}-\varphi_{i-1}}{2}$ (34)

then update with

$\varphi_{i}=\frac{\varphi_{i+1}+\varphi_{i-1}}{2}+\frac{\beta}{2}\left(s(\varphi_{i})\|v_{i}\|^{2}-(s(\varphi_{i})\cdot
v_{i})v_{i}\right)$ (35)

##  Appendix B Experimental Details

Figure 5: The 5 two-dimensional datasets.

##### Software libraries

We used the numpy [Harris et al., 2020] and scipy [Virtanen et al., 2020]
Python libraries for basic numerical and graph operations and matplotlib
[Hunter, 2007] for plotting. We used scikit-learn [Pedregosa et al., 2011] to
compute nearest neighbor graphs and pytorch [Paszke et al., 2019] for training
neural networks. We used FrEIA [Ardizzone et al., 2018-2022] for constructing
normalizing flows.

##### Two-dimensional datasets

We use the following two-dimensional datasets. All are implemented as Gaussian
mixture models. The first two are explicitly Gaussian, the final three are
GMMs with 50 components fitted using scikit-learn to data generated by adding
noise to certain geometric structures.

  * •

Standard normal

  * •

Gaussian mixture with 3 components

  * •

Circle

  * •

1 spiral

  * •

2 spirals

##### Normalizing flow training

We used the FrEIA library [Ardizzone et al., 2018-2022] to design normalizing
flows using spline coupling blocks [Durkan et al., 2019]. We used the
following hyperparameters:

  * •

Blocks: 12 for two-dimensional data, 5 for standard normal in $D$ dimensions

  * •

Spline bins: 10

  * •

Domain clamping in the splines: 5

  * •

ActNorm between blocks

  * •

Noise of $10^{-5}$ added to the training data

  * •

LR = $5\times 10^{-4}$

  * •

Weight decay = $10^{-6}$

  * •

Fully-connected coupling block subnets with 3 hidden layers of 64 dimensions
each, with ReLU activation and BatchNorm

  * •

Training iterations (as a function of training set size $n$): $2500\times
2^{n//1000}$

  * •

Batch size = 256 if $n>400$, else $\lfloor n/3\rfloor$

  * •

Adam optimizer

##### Score model training

We use fully connected networks to predict the score, with 6 hiddens layer
with the width depending on the dimension of the input. We use the following
hyperparameters:

  * •

Hidden layer width (as a function of input dimension $D$)

    * –

$3\leq D\leq 5$: 128

    * –

$6\leq D\leq 8$: 170

    * –

$9\leq D\leq 25$: 200

  * •

LR = $10^{-3}$

  * •

Weight decay = $10^{-6}$

  * •

Softplus activation function

  * •

Training iterations (as a function of training set size $n$): $2500\times
2^{n//1000}$

  * •

Batch size = 256 if $n>400$, else $\lfloor n/3\rfloor$

  * •

Adam optimizer

In all cases, the chosen model checkpoint was that which achieved the lowest
training loss.

##### Compute resources

For the two-dimensional datasets, we performed all experiments on a machine
with a NVIDIA GeForce RTX 2070 GPU and 32 GB of RAM. It took approximately 24
hours to run all experiments.

For the higher dimensional standard normal distributions, we performed all
experiments on a machine with two GPUs: a NVIDIA GeForce RTX 2070 and a NVIDIA
GeForce RTX 2080 Ti Rev. A. The machine has 32 GB of RAM. It took
approximately 24 hours to run all experiments.

In addition we did some preliminary experiments on the first machine,
totalling approximately 48 hours of compute time.

##  Appendix C Additional experimental results

Figure 6: Increasing the number of nearest neighbors does not significantly
improve the convergence of the power-weighted graph method on two-dimensional
datasets.

###  C.1 Convergence of alternative graph edge weightings

Figure 7: Alternative edge weightings based on the nearest neighbor density
estimator do not converge. The same weightings based on the ground truth
density do converge. Experiments performend on the GMM two-dimensional
dataset.

We try out other edge weights, as described in the main text. Given the
nearest-neighbor density estimates
$\hat{p}(X_{l})\propto(\min_{m}\|X_{l}-X_{m}\|)^{-d}$ for data points $X_{l}$,
we approximate the graph edges by
$\mathcal{G}_{lm}=\|X_{l}-X_{m}\|/\tilde{p}((X_{l}+X_{m})/2)$ where
$\tilde{p}$ is a density estimate for the midpoint. We use 4 variants:

  1. 1.

Inverse of mean:
$\tilde{p}((X_{l}+X_{m})/2)=(\hat{p}(X_{l})+\hat{p}(X_{m}))/2$

  2. 2.

Mean of inverse:
$1/\tilde{p}((X_{l}+X_{m})/2)=(1/\hat{p}(X_{l})+1/\hat{p}(X_{m}))/2$

  3. 3.

Max: $\tilde{p}((X_{l}+X_{m})/2)=\max(\hat{p}(X_{l}),\hat{p}(X_{m}))$

  4. 4.

Min: $\tilde{p}((X_{l}+X_{m})/2)=\min(\hat{p}(X_{l}),\hat{p}(X_{m}))$

We find that all 4 have very similar performance to the power-weighted graph,
whereas the equivalent estimators using the ground truth density quickly
converge. See Fig. 7.

Figure 8: Detailed plot of convergence behavior in higher dimensions, unscaled
metric. We see that all graph-based methods plateau in performance, whereas
the relaxation-based method continues to show good convergence. Figure 9:
Detailed plot of convergence behavior in higher dimensions, scaled metric. We
see that all graph-based methods plateau in performance, whereas the
relaxation-based method continues to show good convergence.


††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.††thanks: These authors contributed
equally to this work.

# Proximity effect induced intriguing superconductivity in van der Waals
heterostructure of magnetic topological insulator and conventional
superconductor

Peng Dong  Xiang Zhou  Xiaofei Hou  Jiadian He  Yiwen Zhang  Yifan Ding
Xiaohui Zeng  Jinghui Wang  Yueshen Wu  ShanghaiTech Laboratory for
Topological Physics & School of Physical Science and Technology, ShanghaiTech
University, Shanghai 201210, China  Kenji Watanabe  Research Center for
Functional Materials, National Institute for Materials Science, Tsukuba
305-0044, Japan  Takashi Taniguchi  International Center for Materials
Nanoarchitectonics, National Institute for Materials Science, Tsukuba
305-0044, Japan  Wei Xia<EMAIL_ADDRESS>Yanfeng Guo
ShanghaiTech Laboratory for Topological Physics & School of Physical Science
and Technology, ShanghaiTech University, Shanghai 201210, China  Yulin Chen
ShanghaiTech Laboratory for Topological Physics & School of Physical Science
and Technology, ShanghaiTech University, Shanghai 201210, China  Department of
Physics, Clarendon Laboratory, University of Oxford, Oxford OX1 3PU, UK  Wei
Li<EMAIL_ADDRESS>State Key Laboratory of Surface Physics and Department
of Physics, Fudan University, Shanghai 200433, China  Jun Li
<EMAIL_ADDRESS>ShanghaiTech Laboratory for Topological Physics &
School of Physical Science and Technology, ShanghaiTech University, Shanghai
201210, China

(August 28, 2024)

###### Abstract

Nontrivial topological superconductivity has received enormous research
attentions due to its potential for diverse applications in topological
quantum computing. The intrinsic issue concerning the correlation between a
topological insulator and a superconductor is, however, still widely open.
Here, we systemically report an emergent superconductivity in a cross-junction
composed of a magnetic topological insulator MnBi2Te4 and a conventional
superconductor NbSe2. Remarkably, the interface indicates existence of a
reduced superconductivity at surface of NbSe2 and a proximity-effect-induced
superconductivity at surface of MnBi2Te4. Furthermore, the in-plane angular-
dependent magnetoresistance measurements reveal the fingerprints of the paring
symmetry behaviors for these superconducting gaps as a unconventional nature.
Our findings extend our views and ideas of topological superconductivity in
the superconducting heterostructures with time-reversal symmetry breaking,
offering an exciting opportunity to elucidate the cooperative effects on the
surface state of a topological insulator aligning a superconductor.

Topological superconductors (TSC), which can host Majorana quasiparticles
Frolov  _et al._ (2020); Sato and Ando (2017), are generally considered as a
promising path to obey non-Abelian statistics and encode and manipulate
quantum information in a topologically protected manner Clarke and Wilhelm
(2008); Zhu  _et al._ (2021). Commonly, the superconductivity gap symmetry of
TSC behaves a nontrivial topological nature and the pair potential of Cooper
pairs exhibits an anisotropic orbital symmetry, such as a chiral _p_ -wave
superconductivity Ivanov (2001). Up to now, the most persuasive way to achieve
TSC is through the conventional _s_ -superconductor proximate to a topological
insulator (TI) Fu and Kane (2008); Sun  _et al._ (2016); Lutchyn  _et al._
(2018); Kim  _et al._ ; Mourik  _et al._ (2012); Nadj-Perge  _et al._ (2014);
Jäck  _et al._ (2019). Due to the inevitable degradation of superconducting
wavefunction into the TIs in the presence of the interfacial electronic states
Liu  _et al._ (2016); Geim and Grigorieva (2013); Rhodes  _et al._ (2019);
Novoselov  _et al._ (2016), interface-control is challenging for topological
superconducting heterojunction engineering.

Recently, the van de Waals magnet MnBi2Te4 (MBT) has attracted great attention
owing to the existence of the long-range magnetic order in the nontrivial
topological phase He (2020); Liu  _et al._ (2020); Otrokov  _et al._ (2019);
Li  _et al._ ; Deng  _et al._ (2020). MBT crystallizes an alternatively
stacking topological insulating Bi2Te3 layers and metallic MnTe layers. As a
consequence, the bulk MBT behaves as a Weyl semimetal Li  _et al._ . In
addition, MBT is an A-type antiferromagnet (AFM) in which the Mn2+ ions are
ferromagnetically ordered within a single layer (SL) while
antiferromagnetically coupled between layers when the temperature is lower
than the N$\acute{e}$el temperature. Therefore, of particular interest is that
when the topological MBT flake with time-reversal symmetry breaking is
proximate to a conventional superconductor (SC), the heterointerface will
provide an ideal path to understanding the emergence of rich topological
superconducting phenomena. Historically, it has been proposed in theory that a
heterojunction comprised of MBT and a conventional $s$-wave SC probably
induces a TSC associated with chiral Majorana edge modes. In experiments, the
superconducting thin films of Nb and NbN have been grown on MBT Chen  _et al._
(2022); Xu  _et al._ (2022), and a proximity-effect-induced superconductivity
gap of about 0.1 meV has been observed. More interestingly, a clear Coulomb
blockade oscillation has been perceived. Nevertheless, the understanding of
the surface state of MBT in proximity to a superconductor remains a complex
and challenging area of research, and many questions still lack comprehensive
answers.

In this work, we fabricated MBT/NbSe2 heterojunctions by stacking the single
crystalline MBT and NbSe2. By selecting different electrodes on the
heterojunctions, we simultaneously studied the transport properties of the
NbSe2, the MBT, the whole heterojunction, and the interface. Intriguingly, the
emergent three superconducting gaps have been perceived through the
measurements on the heterojunctions, including the intrinsic NbSe2, the
reduced superconductivity at the surface of the NbSe2, and the proximity-
effect-induced superconductivity at the surface of the MBT.

Fig. 1(a) shows the optical microscopic image and the schematic diagram of the
MBT/NbSe2 cross-junction, where the superconducting NbSe2flake ($\sim$ 30 nm
in thickness) was stacked on top of the magnetic topological MBT flake ($\sim$
40 nm), configuring a “cross-like” junction, and a 20 nm thick insulating
$h$-BN was laid on the bottom of the heterostructure as a substrate. The
electrodes on MBT and NbSe2 flakes can be selected in various channels to
reveal the electronic properties of the NbSe2, the MBT, the whole
heterojunction and its heterointerface.

Considering the interface between the MBT and NbSe2, as given in Fig. 1(b),
MBT reveals a ferromagnetic order along the out-of-plane, while electrons in
NbSe2 are condensed into Cooper pairs which behave as a singlet
superconducting state ($\uparrow\downarrow$-$\downarrow\uparrow$) below the
critical temperature ($T_{c}$). As has been demonstrated beforeBuzdin (2005);
Xi  _et al._ (2016); Hess  _et al._ (1989); Hamill  _et al._ (2021); Yokoya
_et al._ (2001), when the distance from NbSe2 increases, its superconducting
gap $\Delta_{0}$ will gradually degenerate and form a reduced
superconductivity gap $\Delta_{1}$ at the interface, as illustrated in Fig.
1(c). Notably, the gap $\Delta_{1}$ strongly depends on the electron state of
the non-SC side. Once the non-SC side is a normal metal, the superconducting
wave function from the NbSe2 can flow into the metal within a distance less
than the coherence length $\xi_{\text{N}}=\sqrt{\hbar D/k_{\text{B}}T}$, where
$D$ is the diffusion coefficient Buzdin (2005). While for a magnetic material,
the superconductivity coherence length ($\xi_{\text{M}}$) turns into
$\xi_{\text{M}}=\sqrt{\hbar D_{\text{M}}/k_{\text{B}}T_{\text{M}}}$, where
$T_{\text{M}}$ is the magnetic transition temperature. Meanwhile, the
superconducting wave function decays exponentially
$\Delta_{\text{M}}$=exp$(D_{\text{M}}\textfractionsolidus\xi_{\text{M}})$ due
to the existence of the exchange splitting from the spin-spin interaction,
resulting in a restricted superconducting gap function ($\Delta_{2}$). Since
the exchange splitting is absent in an AFM, the superconductivity coherence
length should be comparable to that of normal metal. However, due to the
A-type AFM magnetic structure of MBT, it is essential to take into account the
ferromagnetic coupling when considering the interface between MBT and NbSe2.

Figure 1: (color online). Sample geometry and the interface of MBT/NbSe2. (a)
The optical microscopic image of the MBT/NbSe2 structure. The bottom h-BN ( 20
nm) provides an atomic flat substrate, and MBT (40 nm) and NbSe2 (30 nm)
flakes are stacked as a heterojunction, in which the junction area is labeled
by a white dash line. Nine electrodes marked from 1 to 9 are connected in
different configurations of current flowing and voltage measurement for
comparison of electrical properties of the pristine NbSe2 (current 1-5,
voltage 2-3), MBT (6-9, 7-8), the whole heterojunction (1-6, 3-8), and the
interface (1-6, 5-9) channels. (b) Schematic images of MBT/NbSe2
heterojunction. MBT reveals an antiferromagnetic order along the out-of-plane,
while electrons in NbSe2 demonstrate Cooper pairs with opposite spin
directions in the superconducting state. (c) Schematic image of the
degradation of the superconducting wave function ($\psi$) along the distance
from NbSe2. Here, $\Delta_{0}$, $\Delta_{1}$, and and $\Delta_{2}$ correspond
to the superconductivity gaps of the intrinsic NbSe2, the surface
superconductivity of NbSe2, and the proximity-effect induced superconductivity
at the MBT surface. (d) The equivalent circuit of the junction.

By selecting electrodes on the MBT and NbSe2 flakes in different channels (see
Fig. 1(a)), we were able to evaluate the electronic properties of NbSe2, MBT,
the whole heterojunction, and the interface are illustrated in Fig. 2(a),
2(e), 2(i), and 2(m). The detailed analysis is introduced in the Supplementary
Materials.

For the pristine NbSe2, the temperature- and magnetic field-dependent
resistances are shown in Fig. 2(b) and 2(c), respectively, indicating the
intrinsic superconductivity of NbSe2. Specifically, a superconducting
transition temperature ($T_{\text{c}}$) of 7 K and a critical field
($B_{\text{c}}$) of 3 T at 1.8 K were observed. Moreover, the critical current
of NbSe2 is about 1.8 mA at 1.8 K, as shown in Fig. 2(d). According to the
electron-phonon coupling strength in BCS model
$2\Delta_{0}=3.52k_{\text{B}}T_{\text{c}}$, we can estimate the
superconducting gap of the pristine NbSe2 to be 1.2 meV at $T$ = 0, which is
consistent with the previous reportsHess  _et al._ (1989).

The transport properties of MBT are consistent with previous resultsOtrokov
_et al._ (2019); Zeugner  _et al._ (2019); Cui  _et al._ (2019); Lee  _et al._
(2019), where the N$\acute{e}$el temperature ($T_{\text{N}}$) is 23.7 K (see
Fig. 2(f)) and the spin-flop critical field ($B_{\text{c}}^{\text{MBT}}$) is
$\sim$3.34 T (see Fig. 2(g)). It seems that the $B_{\text{c}}^{\text{MBT}}$ is
comparable to that of the $B_{\text{c}}$ of NbSe2 ($\sim$3 T). While, one can
distinguish the contribution of MBT or NbSe2 because these transitions induce
opposite modifications on the magnetoresistance, namely, MBT causes a downturn
while NbSe2 leads to an upturn as the magnetic field is increased (see Fig.
2(c) and Fig. 2(g), respectively).

Figure 2: (color online). Transport properties for different channels of the
heterojunction. (a), (e), (i), and (m). Measurement schematic of the NbSe2
(current 1-5, voltage 2-3 as shown in Fig. 1), MBT (6-9, 7-8), the whole
heterojunction (1-6, 3-8), and the interface (1-6, 5-9) channels. (b), (f),
(j), and (n). Temperature-dependent resistance for NbSe2, MBT, MBT/NbSe2
heterojunction, and interface, respectively. Here, all $R-T$ curves are
measured under zero field. Particularly, two temperature transitions are
identified at the R-T curve of the heterojunction and interface, including the
superconducting transition Tc of NbSe2 at 7 K and antiferromagnetic transition
$T_{\text{N}}$ of MBT at 23.7 K. The corresponding magnetoresistances at 1.8 K
are given in (c), (g), (k), and (o), where the field was applied along the
out-of-plane. The magnetic transition of MBT is observed in (g) and (k) as the
spin-flop critical field ($B_{\text{c}}^{\text{MBT}}$), and the critical field
of NbSe2 exhibits in (c) as $B_{\text{N}}$. (d), (h), (l), and (p). The
$dV/dI$ vs. $I$ curves at 1.8 K for NbSe2, MBT, MBT/NbSe2 heterojunction, and
the interface, respectively. Here, the current axis for NbSe2 in (d) is timed
0.6, because the critical current of NbSe2 is 1.8 mA.

With the whole heterojunction measurement configuration (see Fig. 2(i)), the
temperature-dependent resistance shows two characteristic transitions at 23.7
and 7 K in Fig. 2(j), which correspond to the $T_{\text{N}}$ of MBT and the
$T_{\text{c}}$ of NbSe2, respectively. The magnetoresistance of the
heterojunction in Fig. 2(k) reveals a downturn at 3.34 T following a nonlinear
upturn transition at around 3 T, which can be related to the critical field of
MBT ($B_{\text{c}}^{\text{MBT}}$) and NbSe2 ($B_{\text{c}}$), respectively. On
the other hand, the $dV/dI$ vs. $I$ curve of the heterojunction at 1.8 K (see
Fig. 2(i)) demonstrates two pairs of transitions, locating at around $\pm$895
and $\pm$ 245 $\mu$A, respectively, which will be discussed in detail in the
following part.

The transport properties of the MBT/NbSe2 interface were investigated by
applying a current through two adjacent sides of the “cross-shape”
heterojunction and measuring the voltage through the other two sides (see Fig.
2(m)). The temperature dependence of resistance in Fig. 2(n) also shows two
distinct transitions at 23.7 and 7 K, corresponding to the $T_{\text{N}}$ and
$T_{\text{c}}$, respectively. However, only an upturn transition is observed
on the magnetoresistance curve as shown in Fig. 2(o), suggesting that the
superconducting transition dominates the magnetoresistance at the interface.
It is worth noting that the resistance of the interface exhibits an upturn
behavior below the superconductivity transition, which can be attributed to
the Andreev reflection (AR) effect. Such AR phenomenon also contributes to the
peak of differential resistance at the zero-bias current as shown in Fig.
2(p). On the other hand, two pairs of transitions are observed at $\pm$878 and
$\pm 248\mu$A on the $dV/dI$ \- $I$ curve, which are consistent with those of
heterojunction. However, these critical currents are significantly lower than
that of the pristine NbSe2, which is 1.8 mA at 1.8 K. These features can be
attributed to the superconducting transitions of the reduced superconductivity
at the surface of NbSe2 and the proximity-effect-induced areas on the surface
of MBT (see Fig. 1(d)).

Figure 3: (color online). Temperature-dependent critical current. Measurement
schematic of (a) the whole MBT/NbSe2 heterojunction and (d) the MBT/NbSe2
interface. The differential resistance mapping of (b) the heterojunction and
(e) the interface with the temperature ranging from 1.8 K to 7.0 K. The
selected dV/dI-I spectrums of (c) the heterojunction and (f) the interface
with several typical temperatures. $I_{\text{c}}^{0}$, $I_{\text{c}}^{1}$, and
$I_{\text{c}}^{2}$ correspond to the critical currents of the intrinsic NbSe2,
the reduced superconductivity at the surface of the NbSe2, and the proximity-
effect-induced superconductivity at the surface of the MBT.  Figure 4: (color
online). Angular-dependent critical current at 1.8 K. (b) Schematic image of
the measurement setup. Here, the magnetic field was applied within the in-
plane, and $\theta$ corresponds to the angle between the current within the
NbSe2 and field $B$. (a) and (c). The angular-dependent differential
resistance mapping under different magnetic fields for the whole
heterojunction and the interface, respectively. (d) and (e). The angular-
dependent of $I_{\text{c}}^{0}$, $I_{\text{c}}^{1}$, and $I_{\text{c}}^{2}$ at
4 T for the whole heterojunction and the interface, respectively. Here, the
$I_{\text{c}}^{0}$, $I_{\text{c}}^{1}$, and $I_{\text{c}}^{2}$ are estimated
from mapping in (a) and (b).

To understand the proximity-effect-induced superconductivity and AR through
the heterojunction, we further studied the temperature-dependent $dV/dI-I$
spectrums of the whole heterojunction and the interface. Fig. 3(b) shows the
mapping of normalized differential resistance at various temperatures of the
whole heterostructure, and the differential resistances at some typical
temperatures are given in Fig. 3(c). Three characteristic contour lines,
labeled as $I_{\text{c}}^{0}$, $I_{\text{c}}^{1}$, and $I_{\text{c}}^{2}$, can
be seen on both positive and negative current sides of the mapping in Fig.
3(c), corresponding to the peaks of the $dV/dI$-$I$ curves shown in Fig. 3(d).
It is worthy to mention that the temperature-dependent $I_{\text{c}}^{0}$,
$I_{\text{c}}^{1}$, and $I_{\text{c}}^{2}$ and their corresponding gaps
($\Delta_{\text{c}}^{0}$, $\Delta_{\text{c}}^{1}$, and
$\Delta_{\text{c}}^{2}$, respectively) are consistent with the conventional
BCS model as
$\Delta_{\text{T}}=1.74\Delta_{\text{c}}^{0}\sqrt{1-T/T_{\text{c}}}$.
Therefore, $I_{\text{c}}^{0}$, $I_{\text{c}}^{1}$, and $I_{\text{c}}^{2}$, can
be attributed to the intrinsic NbSe2, the reduced superconductivity at the
surface of the NbSe2, and the proximity-effect-induced superconductivity at
the surface of the MBT.

In comparison, we perform the same measurement on the interface (see Fig.
3(d)). The mapping results shown in Fig. 3(e) reveal the absence of
$I_{\text{c}}^{0}$, which further confirms that $I_{\text{c}}^{0}$ corresponds
to the intrinsic NbSe2 since it is excluded during the measurement setup.
Moreover, the transition of $I_{\text{c}}^{1}$ and $I_{\text{c}}^{2}$, and
particularly the AR, are considerably more intense than those of the
heterojunction. The main reason is that both NbSe2 and MBT sides are excluded
in this configuration, and the interfacial effect dominates the measurement
result. Therefore, the “cross-junction” geometry is an ideal configuration to
investigate the interfacial effects.

Considering the AR effect at an SC/non-SC interface, it is an effective method
to determine the polarization $P$ of a certain material Soulen  _et al._
(1998); Strijkers  _et al._ (2001). In our present AR effect in the MBT/NbSe2
interface, the reduction of conductance at zero-bias current indicates weak
spin-polarization from MBT. However, the MBT crystal in our case is about 40
nm ($\sim$10 SLs), which should demonstrate a long-term AFM coupling as a bulk
crystal. For comparison, we also studied the AR behavior of the
Fe5GeTe2/NbSe2, where Fe5GeTe2 reveals typical FM coupling along the $c$-axis.
Surprisingly, the $dI/dV$-$V$ spectrums at the interface of MBT/NbSe2 and
Fe5GeTe2/NbSe2 (see Fig. S5) are comparable. Therefore, we can conclude that
the surface state of MBT reveals an FM ordering and partially restricts AR.

According to the temperature- and field-dependent (see Fig. 2 and Figs. S1-S3)
results, we have confirmed the observation of proximity-effect-induced
superconductivity ($\Delta_{1}$) on the surface of MBT. In principle, the AFM
coupling and topological nontrivial state in the electron bands structure may
contribute to the superconducting coupling, resulting in a novel
superconductivity gap symmetry or even topological superconductivity, which
should be pronouncedly different from the conventional superconductivity in
NbSe2. To uncover the symmetry of the proximity-effect-induced
superconductivity in MBT, we also measured the angular-dependent
magnetoresistance spectrums of the whole MBT/NbSe2 heterojunction and the
interface.

The angular-dependent measurement geometry is shown in Fig. 4(b). Figs. 4(a)
and 4(c) show the mapping of angular-dependent $dV/dI$-$I$ spectrums for the
heterostructure and interface at 1.8 K under magnetic fields of 1, 4, and 9 T.
In particular, the $I_{\text{c}}^{0}$, $I_{\text{c}}^{1}$, and
$I_{\text{c}}^{2}$ can be identified. Moreover, the corresponding angular-
dependent normalized critical currents are plotted in Figs. 4(d) and 4(e),
where $\theta$ is the angle between the magnetic field and $a$-axis direction
of NbSe2. Surprisingly, in contrast to what is generally considered as a
conventional $s$-wave superconductor that should behave as an isotropic
superconductivity symmetry, the angular-dependent $I_{\text{c}}^{0}$ of NbSe2
reveals a cosine function with two-fold symmetry, where the maximum appears at
0∘ and 180∘ and the minimum appears at 90∘ and 270∘, respectively. Thus, an
unconventional pairing mechanism may exist in finite layers of NbSe2, where
spontaneous nematic superconductivity or strong gap-mixing triggered by a
small symmetry-breaking field may induce an anisotropic nature Huang  _et al._
(2007). However, a deeper understanding of the symmetry behavior is still
needed.

The angular dependence of $I_{\text{c}}^{1}$ is comparable for both
heterostructure and interface measurement configurations (see Figs. 4(d) and
4(e)). Nevertheless, it does not appear to follow a cosine or sine function,
but rather a combination of trigonometric function with a four-fold symmetry.
Although, the minimum appears at about 90∘ and 270∘ as that of NbSe2, the
maximum appears at 120∘ and 330∘ ($B$ = 4 T), indicating that the phase of
$I_{\text{c}}^{1}$ is shifted compared to $I_{\text{c}}^{0}$. Moreover, the
$I_{\text{c}}^{2}-\theta$ curves also reveal a two-fold symmetry for both
measurement configurations, while the pattern is slightly different from those
of $I_{\text{c}}^{0}$ and $I_{\text{c}}^{1}$. Particularly, the minimum of
$I_{\text{c}}^{2}$ appears at $\sim$45∘ and 270∘, which reveals a 45∘ phase
shift from both $I_{\text{c}}^{0}$ and $I_{\text{c}}^{1}$. However, the
maximum seems to correspond to that of $I_{\text{c}}^{0}$ and
$I_{\text{c}}^{1}$. Based on the angular-dependent profiles of
$I_{\text{c}}^{0}$, $I_{\text{c}}^{1}$, and $I_{\text{c}}^{2}$, although we
can hardly confirm the exact gap symmetry of the induced superconductivity on
the surface of MnBi2Te4 or even the intrinsic NbSe2, the superconductivity gap
symmetries of $\Delta_{1}$ and $\Delta_{2}$ are basically different from that
of NbSe2 ($\Delta_{0}$).

In summary, we studied the van der Waals “cross-shape” heterojunction between
a nontrivial TI MBT and a conventional SC NbSe2. Through different measurement
configurations, three different superconductivity gaps have been observed,
including the intrinsic NbSe2, the reduced superconductivity at the surface of
the NbSe2, and the proximity-effect-induced superconductivity at the surface
of the MBT. A clear zero-bias peak was observed in the $dV/dI-I$ curve, which
is corresponding to the AR effect at a FM/SC interface, suggesting that the
surface of MBT exhibits a ferromagnetic plane. The temperature-, field- and
angular-dependent $dV/dI-I$ spectrums with the whole heterojunction and the
interface measurement configurations have been systematically investigated.
Although we cannot identify the superconductivity pairing symmetry, an
unconventional superconducting state which is different from the conventional
SC NbSe2 should be seriously considered, as the various quantum states like
topological states and magnetic ordering may contribute to the electronic
coupling.

This research was supported in part by the Ministry of Science and Technology
(MOST) of China (No. 2022YFA1603903), the National Natural Science Foundation
of China (Grants No. 12004251, 12104302, 12104303, 12304217), the Science and
Technology Commission of Shanghai Municipality, the Shanghai Sailing Program
(Grant No. 21YF1429200), the start-up funding from ShanghaiTech University,
and Beijing National Laboratory for Condensed Matter Physics, the
Interdisciplinary Program of Wuhan National High Magnetic Field Center
(WHMFC202124), the open project from State Key Laboratory of Surface Physics
and Department of Physics, Fudan University (Grant No. KF2022-13), the
Shanghai Sailing Program (23YF1426900), and the Double First-Class Initiative
Fund of ShanghaiTech University. Growth of $h$-BN was supported by the
Elemental Strategy Initiative conducted by the MEXT, Japan (JPMXP0112101001),
JSPS KAKENHI (JP20H00354) and A3 Foresight by JSPS.

## References

  * Frolov  _et al._ (2020) S. M. Frolov, M. J. Manfra, and J. D. Sau, Nature Physics 16, 718 (2020).
  * Sato and Ando (2017) M. Sato and Y. Ando, Reports on Progress in Physics 80, 076501 (2017).
  * Clarke and Wilhelm (2008) J. Clarke and F. K. Wilhelm, Nature 453, 1031 (2008).
  * Zhu  _et al._ (2021) D. Zhu, T. Jaako, Q. He, and P. Rabl, Physical Review Applied 16, 014024 (2021).
  * Ivanov (2001) D. A. Ivanov, Physical Review Letters 86, 268 (2001).
  * Fu and Kane (2008) L. Fu and C. L. Kane, Physical Review Letters 100, 096407 (2008).
  * Sun  _et al._ (2016) H.-H. Sun, K.-W. Zhang, L.-H. Hu, C. Li, G.-Y. Wang, H.-Y. Ma, Z.-A. Xu, C.-L. Gao, D.-D. Guan, Y.-Y. Li, C. Liu, D. Qian, Y. Zhou, L. Fu, S.-C. Li, F.-C. Zhang, and J.-F. Jia, Physical Review Letters 116, 257003 (2016).
  * Lutchyn  _et al._ (2018) R. M. Lutchyn, E. P. A. M. Bakkers, L. P. Kouwenhoven, P. Krogstrup, C. M. Marcus, and Y. Oreg, Nature Reviews Materials 3, 52 (2018).
  * (9) H. Kim, A. Palacio-Morales, T. Posske, L. Rózsa, K. Palotás, L. Szunyogh, M. Thorwart, and R. Wiesendanger, Science Advances 4, eaar5251.
  * Mourik  _et al._ (2012) V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Science 336, 1003 (2012).
  * Nadj-Perge  _et al._ (2014) S. Nadj-Perge, I. K. Drozdov, J. Li, H. Chen, S. Jeon, J. Seo, A. H. MacDonald, B. A. Bernevig, and A. Yazdani, Science 346, 602 (2014).
  * Jäck  _et al._ (2019) B. Jäck, Y. Xie, J. Li, S. Jeon, B. A. Bernevig, and A. Yazdani, Science 364, 1255 (2019).
  * Liu  _et al._ (2016) Y. Liu, N. O. Weiss, X. Duan, H.-C. Cheng, Y. Huang, and X. Duan, Nature Reviews Materials 1, 16042 (2016).
  * Geim and Grigorieva (2013) A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013).
  * Rhodes  _et al._ (2019) D. Rhodes, S. H. Chae, R. Ribeiro-Palau, and J. Hone, Nature Materials 18, 541 (2019).
  * Novoselov  _et al._ (2016) K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H. Castro Neto, Science 353, aac9439 (2016).
  * He (2020) K. He, npj Quantum Materials 5, 90 (2020).
  * Liu  _et al._ (2020) C. Liu, Y. Wang, H. Li, Y. Wu, Y. Li, J. Li, K. He, Y. Xu, J. Zhang, and Y. Wang, Nature Materials 19, 522 (2020).
  * Otrokov  _et al._ (2019) M. M. Otrokov, I. I. Klimovskikh, H. Bentmann, D. Estyunin, A. Zeugner, Z. S. Aliev, S. Gaß, A. U. B. Wolter, A. V. Koroleva, A. M. Shikin, M. Blanco-Rey, M. Hoffmann, I. P. Rusinov, A. Y. Vyazovskaya, S. V. Eremeev, Y. M. Koroteev, V. M. Kuznetsov, F. Freyse, J. Sánchez-Barriga, I. R. Amiraslanov, M. B. Babanly, N. T. Mamedov, N. A. Abdullayev, V. N. Zverev, A. Alfonsov, V. Kataev, B. Büchner, E. F. Schwier, S. Kumar, A. Kimura, L. Petaccia, G. Di Santo, R. C. Vidal, S. Schatz, K. Kißner, M. Ünzelmann, C. H. Min, S. Moser, T. R. F. Peixoto, F. Reinert, A. Ernst, P. M. Echenique, A. Isaeva, and E. V. Chulkov, Nature 576, 416 (2019).
  * (20) J. Li, Y. Li, S. Du, Z. Wang, B.-L. Gu, S.-C. Zhang, K. He, W. Duan, and Y. Xu, Science Advances 5, eaaw5685.
  * Deng  _et al._ (2020) Y. Deng, Y. Yu, M. Z. Shi, Z. Guo, Z. Xu, J. Wang, X. H. Chen, and Y. Zhang, Science 367, 895 (2020).
  * Chen  _et al._ (2022) J. Chen, W. Xu, Z. Tan, Z. Pan, P. Zhu, Z.-M. Liao, and D. Yu, Nano Letters 22, 6484 (2022).
  * Xu  _et al._ (2022) W.-Z. Xu, C.-G. Chu, Z.-C. Pan, J.-J. Chen, A.-Q. Wang, Z.-B. Tan, P.-F. Zhu, X.-G. Ye, D.-P. Yu, and Z.-M. Liao, Physical Review B 105, 184515 (2022).
  * Buzdin (2005) A. I. Buzdin, Reviews of modern physics 77, 935 (2005).
  * Xi  _et al._ (2016) X. Xi, Z. Wang, W. Zhao, J.-H. Park, K. T. Law, H. Berger, L. Forró, J. Shan, and K. F. Mak, Nature Physics 12, 139 (2016).
  * Hess  _et al._ (1989) H. F. Hess, R. B. Robinson, R. C. Dynes, J. M. Valles, and J. V. Waszczak, Physical Review Letters 62, 214 (1989).
  * Hamill  _et al._ (2021) A. Hamill, B. Heischmidt, E. Sohn, D. Shaffer, K.-T. Tsai, X. Zhang, X. Xi, A. Suslov, H. Berger, L. Forró, F. J. Burnell, J. Shan, K. F. Mak, R. M. Fernandes, K. Wang, and V. S. Pribiag, Nature Physics 17, 949 (2021).
  * Yokoya  _et al._ (2001) T. Yokoya, T. Kiss, A. Chainani, S. Shin, M. Nohara, and H. Takagi, Science 294, 2518 (2001).
  * Zeugner  _et al._ (2019) A. Zeugner, F. Nietschke, A. U. B. Wolter, S. Gaß, R. C. Vidal, T. R. F. Peixoto, D. Pohl, C. Damm, A. Lubk, R. Hentrich, S. K. Moser, C. Fornari, C. H. Min, S. Schatz, K. Kißner, M. Ünzelmann, M. Kaiser, F. Scaravaggi, B. Rellinghaus, K. Nielsch, C. Hess, B. Büchner, F. Reinert, H. Bentmann, O. Oeckler, T. Doert, M. Ruck, and A. Isaeva, Chemistry of Materials 31, 2795 (2019).
  * Cui  _et al._ (2019) J. Cui, M. Shi, H. Wang, F. Yu, T. Wu, X. Luo, J. Ying, and X. Chen, Physical Review B 99, 155125 (2019).
  * Lee  _et al._ (2019) S. H. Lee, Y. Zhu, Y. Wang, L. Miao, T. Pillsbury, H. Yi, S. Kempinger, J. Hu, C. A. Heikes, P. Quarterman, W. Ratcliff, J. A. Borchers, H. Zhang, X. Ke, D. Graf, N. Alem, C.-Z. Chang, N. Samarth, and Z. Mao, Physical Review Research 1, 012011 (2019).
  * Soulen  _et al._ (1998) R. J. Soulen, J. M. Byers, M. S. Osofsky, B. Nadgorny, T. Ambrose, S. F. Cheng, P. R. Broussard, C. T. Tanaka, J. Nowak, J. S. Moodera, A. Barry, and J. M. D. Coey, Science 282, 85 (1998).
  * Strijkers  _et al._ (2001) G. J. Strijkers, Y. Ji, F. Y. Yang, C. L. Chien, and J. M. Byers, Physical Review B 63, 104510 (2001).
  * Huang  _et al._ (2007) C. L. Huang, J. Y. Lin, Y. T. Chang, C. P. Sun, H. Y. Shen, C. C. Chou, H. Berger, T. K. Lee, and H. D. Yang, Physical Review B 76, 212504 (2007).


11institutetext: Charles University in Prague, Faculty of Mathematics and
Physics
11email<EMAIL_ADDRESS>
11email<EMAIL_ADDRESS>
11email<EMAIL_ADDRESS>
http://ufal.mff.cuni.cz/ondrej-platek

# Recurrent Neural Networks for Dialogue State Tracking

Ondřej Plátek  Petr Bělohlávek  Vojtěch Hudeček  Filip Jurčíček

###### Abstract

This paper discusses models for dialogue state tracking using recurrent neural
networks (RNN). We present experiments on the standard dialogue state tracking
(DST) dataset, DSTC2 [7]. On the one hand, RNN models became the state of the
art models in DST, on the other hand, most state-of-the-art DST models are
only turn-based and require dataset-specific preprocessing (e.g.
DSTC2-specific) in order to achieve such results. We implemented two
architectures which can be used in an incremental setting and require almost
no preprocessing. We compare their performance to the benchmarks on DSTC2 and
discuss their properties. With only trivial preprocessing, the performance of
our models is close to the state-of-the-art results.111 Acknowledgment: We
thank Mirek Vodolán and Ondřej Dušek for useful comments. This research was
partly funded by the Ministry of Education, Youth and Sports of the Czech
Republic under the grant agreement LK11221, core research funding, grant GAUK
1915/2015, and also partially supported by SVV project number 260 333. We
gratefully acknowledge the support of NVIDIA Corporation with the donation of
the Tesla K40c GPU used for this research. Computational resources were
provided by the CESNET LM2015042 and the CERIT Scientific Cloud LM2015085,
provided under the programme “Projects of Large Research, Development, and
Innovations Infrastructures”.

##  1 Introduction

Dialogue state tracking (DST) is a standard and important task for evaluating
task-oriented conversational agents [18, 7, 8]. Such agents play the role of a
domain expert in a narrow domain, and users ask for information through
conversation in natural language (see the example system and user responses in
Figure 1). A dialogue state tracker summarizes the dialogue history and
maintains a probability distribution over the (possible) user’s goals (see
annotation in Figure 1). Dialogue agents as introduced in [20] decide about
the next action based on the dialogue state distribution given by the tracker.
User’s goals are expressed in a formal language, typically represented as a
dialogue act items (DAIs) (see Section 2) and the tracker updates probability
for each item. The dialogue state is a latent variable [20] and one needs to
label the conversations in order to train a dialogue state tracker using
supervised learning. It was shown that with a better dialogue state tracker,
conversation agents achieve better success rate in overall completion of the
their task [11].

…
Dial. state n: food:None, area:None, pricerange:None
System: What part of town do you have in mind?
User: West part of town.
Dial. state n+1: food:None, area:west, pricerange:None
System: What kind of food would you like?
User: Indian
Dial. state n+2: food:Indian, area:west, pricerange:None
System: India House is a nice place in the west of town serving tasty Indian
food.
…

Figure 1: Example of golden annotation of Dialogue Act Items (DAIs). The
dialogue act items comprise from act type (all examples have type inform) and
slots (food, area, pricerange) and their values (e.g. Indian, west, None).

This paper compares two different RNN architectures for dialogue state
tracking (see Section 3). We describe state-of-the art word-by-word dialogue
state tracker architectures and propose to use a new encoder-decoder
architecture for the DST task (see Section 4.2).

We focus only on the goal slot predictions because the other groups are
trivial to predict222The slots Requested and Method have accuracies 0.95 and
0.95 on the test set according to the state-of-the-art [19]..

We also experiment with re-splitting of the DSTC2 data because there are
considerable differences between the standard train and test datasets [7].
Since the training, development and test set data are distributed differently,
the resulting performance difference between training and test data is rather
high. Based on our experiments, we conclude that DSTC2 might suggest a too
pessimistic view of the state-of-the-art methods in dialogue state tracking
caused by the data distribution mismatch.

##  2 Dialogue state tracking on DSTC2 dataset

Dialogue state trackers maintain their believes beliefs about users’ goals by
updating probabilities of dialogue history representations. In the DSTC2
dataset, the history is captured by dialogue act items and their
probabilities. A Dialogue act item is a triple of the following form
$(actionType,slotName,slotValue)$.

The DSTC2 is a standard dataset for DST, and most of the state-of-the-art
systems in DST have reported their performance on this dataset [7]. The full
dataset is freely available since January 2014 and it contains 1612 dialogues
in the training set, 506 dialogues in the development set and 1117 dialogues
in the test set.333Available online at http://camdial.org/~mh521/dstc/. The
conversations are manually annotated at the turn level where the hidden
information state is expressed in form of $(actionType,slotName,slotValue)$
based on the domain ontology. The task of the domain is defined by a database
of restaurants and their properties444There are six columns in the database:
name, food, price_range, area, telephone, address.. The database and the
manually designed ontology that captures a restaurant domain are both
distributed with the dataset.

##  3 Models

Our models are all based on a RNN encoder [17]. The models update their hidden
states $h$ after processing each word similarly to the RNN encoder of Žilka
and Jurčíček [16]. The encoder takes as inputs the previous state $h_{t-1}$,
representing history for first $t-1$ words, and features $X_{t}$ for the
current word $w_{t}$. It outputs the current state $h_{t}$ representing the
whole dialogue history up to current word. We use a Gated Recurrent Unit (GRU)
cell [5] as the update function instead of a simple RNN cell because it does
not suffer from the vanishing gradient problem [10]. The model optimizes its
parameters including word embeddings [3] during training.

For each input token, our RNN encoder reads the word embedding of this token
along with several binary features. The binary features for each word are:

  * •

the speaker role, representing either user or system,

  * •

and also indicators describing whether the word is part of a named entity
representing a value from the database.

Since the DSTC2 database is a simple table with six columns, we introduce six
binary features firing if the word is a substring of named entity from the
given column. For example, the word indian will not only trigger the feature
for column $food$ and its value indian but also for column restaurant $name$
and its value indian heaven. The features make the data dense by abstracting
the meaning from the lexical level.

Our model variants differ only in the way they predict goal labels, i.e.,
$food$, $area$ and $pricerange$ from the RNN’s last encoded state.555Accuracy
measure with schedule 2 on slot $food$, $area$ and $pricerange$ about which
users can inform the system is a featured metric for DSTC2 challenge [7]. The
first model predicts the output slot labels independently by employing three
independent classifiers (see Section 3.1). The second model uses a decoder in
order to predict values one after the other from the $h_{T}$ (see Section
3.2).

The models were implemented using the TensorFlow [1] framework.

Figure 2: The joint label predictions using RNN from last hidden state
$h_{T}$. The $h_{T}$ represents the whole dialog history of $T$ words. The RNN
takes as input for each word $i$ an embedding and binary features concatenated
to vector $X_{i}$.

###  3.1 Independent classifiers for each label

The independent model (see Figure 3) consists of three models which predict
$food$, $area$ and $pricerange$ based on the last hidden state $h_{T}$
independently. The independent slot prediction that uses one classifier per
slot is straightforward to implement, but the model introduces an unrealistic
assumption of uncorrelated slot properties. In case of DSTC2 and the Cambridge
restaurant domain, it is hard to believe that, e.g., the slots $area$ and
$pricerange$ are not correlated.

We also experimented with a single classifier which predicts the labels
jointly (see Figure 2) but it suffers from data sparsity of the predicted
tuples, so we focused only on the independent label prediction and encoder-
decoder models.

Figure 3: The RNN encodes the word history into dialogue state $h_{T}$ and
predicts slot values independently.

###  3.2 Encoder-decoder framework

We cast the slot predictions problem as a sequence-to-sequence predictions
task and we use a encoder-decoder model with attention [2] to learn this
representation together with slot predictions (see Figure 4). To our
knowledge, we are the first who used this model for dialogue state tracking.
The model is successfully used in machine translation where it is able to
handle long sequences with good accuracy [2]. In DST, it captures correlation
between the decoded slots easily. By introducing the encoder-decoder
architecture, we aim to overcome the data sparsity problem and the incorrect
independence assumptions.

We employ the encoder RNN cell that captures the history of the dialogue which
is represented as a sequence of words from the user and the system. The words
are fed to the encoder as they appear in the dialogue - turn by turn - where
the user and the system responses switch regularly. The encoder updates its
internal state $h_{T}$ after each processed word. The RNN decoder model is
used when the system needs to generate its output, in our case it is at the
end of the user response. The decoder generates arbitrary length sequences of
words given the encoded state $h_{T}$ step by step. In each step, an output
word and a new hidden state $h_{T+1}$ is generated. The generation process is
finished when a special token End of Sequence (EOS) is decoded. This mechanism
allows the model to terminate the output sequence. The attention part of model
is used at decoding time for weighting the importance of the history.

The disadvantage of this model is its complexity. Firstly, the model is not
trivial to implement666We modified code from the TensorFlow ‘seq2seq‘ module..
Secondly, the decoding time is asymptotically quadratic in the length of the
decoded sequences, but our target sequences are always four tokens long
nevertheless.

Figure 4: Encoder decoder with attention predicts goals.

##  4 Experiments

The results are reported on the standard DSTC2 data split where we used 516
dialogues as a validation set for early stopping [14] and the remaining 1612
dialogues for training. We use 1-best Automatic Speech Recognition (ASR)
transcriptions of conversation history of the input and measure the joint slot
accuracy. The models are evaluated using the recommended measure accuracy [7]
with schedule 2 which skips the first turns where the believe tracker does not
track any values. In addition, our models are also evaluated on a randomly
split DSTC2 dataset (see Section 4.3).

For all our experiments, we train word embeddings of size 100 and use the
encoder state size of size 100, together with a dropout keep probability of
$0.7$ for both encoder inputs and outputs. These parameters are selected by a
grid search over the hyper-parameters on development data.

###  4.1 Training

The training procedure minimizes the cross-entropy loss function using the
Adam optimizer [12] with a batch size of 10. We train predicting the goal slot
values for each turn. We treat each dialogue turn as a separate training
example, feeding the whole dialogue history up to the current turn into the
encoder and predicting the slot labels for the current turn.

We use early stopping with patience [14], validating on the development set
after each epoch and stopping if the three top models does not change for four
epochs.

The predicted labels in DST task depend not only on the last turn, but on the
dialogue full history as well. Since the lengths of dialogue histories vary a
lot777The maximum dialogue history length is 487 words and 95% percentile is
205 words for the training set. and we batch our inputs, we separated the
dialogues into ten buckets accordingly to their lengths in order to provide a
computational speed-up. We reshuffle the data after each epoch only within
each bucket.

In informal experiments, we tried to speed-up the training by optimizing the
parameters only on the last turn888The prediction was conditioned on the full
history but we back-propagated the error only in words within the last turn.
but the performance dropped relatively by more than 40%.

###  4.2 Comparing models

Predicting the labels jointly is quite challenging because the distribution of
the labels is skewed as demonstrated in Figure 5. Some of the labels
combinations are very rare, and they occur only in the development and test
set so the joint model is not able to predict them. During first informal
experiments the joint model performed poorly arguably due to data sparsity of
slot triples. We further focus on model with independent classifiers and
encoder-decoder architecture.

Figure 5: The number occurrences of labels in form of $(food,area,pricerange)$
triples from the least to the most frequent.

The model with independent label prediction is a strong baseline which was
used, among others, in work of Žilka and Jurčíček [16]. The model suffers less
from dataset mismatch because it does not model the correlation between
predicted labels. This property can explain a smaller performance drop between
the test set from reshuffled data and the official test set in comparison to
encoder-decoder model.

Model | Dev set | Test set
---|---|---
Indep | 0.892 | 0.727
EncDec | 0.867 | 0.730
Vodolán et al. [15] | - | 0.745
Žilka and Jurčíček [16] | 0.69 | 0.72
Henderson et al. [6] | - | 0.737
DSTC2 stacking ensemble [7] | - | 0.789
Table 1: Accuracy on DSTC2 dataset. The first group contains our systems which
use ASR output as input, the second group lists systems using also ASR
hypothesis as input. The third group shows the results for ensemble model
using ASR output nd also live language understanding annotations.

Since the encoder-decoder architecture is very general and can predict
arbitrary output sequences, it also needs to learn how to predict only three
slot labels in the correct order. It turned out that the architecture learned
to predict quadruples with three slot values and the EOS symbol quickly, even
before seeing a half of the training data in the first epoch.999We could have
modified the decoder to always predict three symbols for our three slots, but
our experiments showed that the encoder-decoder architecture does not make
mistakes at predicting the order of the three slots and EOS symbol. At the end
of the first epoch, the system made no more mistakes on predicting slot values
in the incorrect order. The encoder-decoder system is competitive with state-
of-the art architectures and the time needed for learning the output structure
was surprisingly short.101010The best model weights were found after 18 to 23
epochs for all model architectures.

###  4.3 Data preparation experiments

Model | Dev set | Test set
---|---|---
Indep | 0.87 | 0.89
EncDec | 0.94 | 0.91
Table 2: Accuracy of our models on the re-split DSTC2 data.

The data for the DSTC2 test set were collected using a different spoken
dialogue system configuration than the data for the validation and the
training set.[7]. We intended to investigate the influence of the complexity
of the task, hence we merged all DSTC2 data together and created splits of
80%, 10% and 10% for the training, development and test sets. The results in
Table 2 show that the complexity of the task dropped significantly.

##  5 Related work

Since there are numerous systems which reported on the DSTC2 dataset, we
discuss only the systems which use RNNs. In general, the RNN systems achieved
excellent results.

Our system is related to the RNN tracker of Žilka and Jurčíček [16], which
reported near state-of-the art results on the DSTC2 dataset and introduced the
first incremental system which was able to update the dialogue state word-by-
word with such accuracy. In contrast to work of Žilka and Jurčíček [16], we
use no abstraction of slot values. Instead, we add the additional features as
described in Section 3. The first system which used a neural network for
dialogue state tracking [6] used a feed-forward network and more than 10
manually engineered features across different levels of abstraction of the
user input, including the outputs of the spoken language understanding
component (SLU). In our work, we focus on simplifying the architecture, hence
we used only features which were explicitly given by the dialogue history word
representation and the database.

The system of Henderson et al. [9] achieves the state-of-the-art results and,
similarly to our system, it predicts the dialogue state from words by
employing a RNN. On the other hand, their system heavily relies on the user
input abstraction. Another dialogue state tracker with LSTM was used in the
reinforcement setting but the authors also used information from the SLU
pipeline [13].

An interesting approach is presented in the work of Vodolán et al. [15], who
combine a rule-based and a machine learning based approach. The handcrafted
features are fed to an LSTM-based RNN which performs a dialog-state update.
However, unlike our work, their system requires SLU output on its input.

It is worth noting that there are first attempts to train an end-to-end
dialogue system even without explicitly modeling the dialogue state [4], which
further simplifies the architecture of a dialogue system. However, the
reported end-to-end model was evaluated only on artificial dataset and cannot
be compared to DSTC2 dataset directly.

##  6 Conclusion

We presented and compared two dialogue state tracking models which are based
on state-of-the-art architectures using recurrent neural networks. To our
knowledge, we are the first to use an encoder-decoder model for the dialogue
state tracking task, and we encourage others to do so because it is
competitive with the standard RNN model.111111The presented experiments are
published at https://github.com/oplatek/e2end/ under Apache license. Informal
experiments were conducted during the Statistical Dialogue Systems course at
Charles University (see https://github.com/oplatek/sds-tracker). The models
are comparable to the state-of-the-art models.

We evaluate the models on DSTC2 dataset containing task-oriented dialogues in
the restaurant domain. The models are trained using only ASR 1-best
transcriptions and task-specific lexical features defined by the task
database. We observe that dialogue state tracking on DSTC2 test set is
notoriously hard and that the task becomes substantially easier if the data is
reshuffled.

As future work, we plan to investigate the influence of the introduced
database features on models accuracy. To our knowledge there is no dataset
which can be used for evaluating incremental dialogue state trackers, so it
would be beneficial to collect the word-level annotations so one can evaluate
incremental DST models.

## References

  * [1] Martın Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, et al.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  _Software available from tensorflow. org_.
  * Bahdanau et al. [2014] Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio.  Neural machine translation by jointly learning to align and translate.  _arXiv preprint arXiv:1409.0473_ , 2014.
  * Bengio et al. [2003] Yoshua Bengio, Rejean Ducharme, and Pascal Vincent.  A Neural probabilistic language model.  _Journal of Machine Learning Research_ , 3:1137–1155, 2003\.
  * Bordes and Weston [2016] Antoine Bordes and Jason Weston.  Learning End-to-End Goal-Oriented Dialog.  _arXiv preprint arXiv:1605.07683_ , 2016.
  * Cho et al. [2014] Kyunghyun Cho, Bart van Merrienboer, Çaglar Gülçehre, Fethi Bougares, Holger Schwenk, and Yoshua Bengio.  Learning Phrase Representations using RNN Encoder-Decoder for Statistical Machine Translation.  _CoRR_ , abs/1406.1078, 2014.  URL http://arxiv.org/abs/1406.1078.
  * Henderson et al. [2013] Matthew Henderson, Blaise Thomson, and Steve Young.  Deep neural network approach for the dialog state tracking challenge.  _Proceedings of the SIGDIAL 2013 Conference_ , pages 467–471, 2013\.
  * Henderson et al. [2014a] Matthew Henderson, Blaise Thomson, and Jason Williams.  The second dialog state tracking challenge.  In _15th Annual Meeting of the Special Interest Group on Discourse and Dialogue_ , volume 263, 2014a.
  * Henderson et al. [2014b] Matthew Henderson, Blaise Thomson, and Jason D Williams.  The third dialog state tracking challenge.  In _Spoken Language Technology Workshop (SLT), 2014 IEEE_ , pages 324–329. IEEE, 2014b.
  * Henderson et al. [2014c] Matthew Henderson, Blaise Thomson, and Steve Young.  Word-based dialog state tracking with recurrent neural networks.  In _Proceedings of the 15th Annual Meeting of the Special Interest Group on Discourse and Dialogue (SIGDIAL)_ , pages 292–299, 2014c.
  * Hochreiter et al. [2001] Sepp Hochreiter, Yoshua Bengio, Paolo Frasconi, and Jürgen Schmidhuber.  Gradient flow in recurrent nets: the difficulty of learning long-term dependencies, 2001.
  * Jurčíček et al. [2012] Filip Jurčíček, Blaise Thomson, and Steve Young.  Reinforcement learning for parameter estimation in statistical spoken dialogue systems.  _Computer Speech & Language_, 26(3):168–192, 2012.
  * Kingma and Ba [2014] Diederik Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Lee and Kim [2016] Byung-Jun Lee and Kee-Eung Kim.  Dialog History Construction with Long-Short Term Memory for Robust Generative Dialog State Tracking.  _Dialogue & Discourse_, 7(3):47–64, 2016.
  * Prechelt [1998] Lutz Prechelt.  Early stopping-but when?  In _Neural Networks: Tricks of the trade_ , pages 55–69. Springer, 1998.
  * Vodolán et al. [2015] Miroslav Vodolán, Rudolf Kadlec, and Jan Kleindienst.  Hybrid Dialog State Tracker.  _CoRR_ , abs/1510.03710, 2015.  URL http://arxiv.org/abs/1510.03710.
  * Žilka and Jurčíček [2015] Lukás Žilka and Filip Jurčíček.  Incremental LSTM-based dialog state tracker.  _arXiv preprint arXiv:1507.03471_ , 2015.
  * Werbos [1990] Paul J Werbos.  Backpropagation through time: what it does and how to do it.  _Proceedings of the IEEE_ , 78(10):1550–1560, 1990.
  * Williams et al. [2013] Jason Williams, Antoine Raux, Deepak Ramachandran, and Alan Black.  The dialog state tracking challenge.  In _Proceedings of the SIGDIAL 2013 Conference_ , pages 404–413, 2013\.
  * Williams [2014] Jason D Williams.  Web-style ranking and SLU combination for dialog state tracking.  In _15th Annual Meeting of the Special Interest Group on Discourse and Dialogue_ , volume 282, 2014.
  * Young et al. [2010] Steve Young, Milica Gašić, Simon Keizer, François Mairesse, Jost Schatzmann, Blaise Thomson, and Kai Yu.  The hidden information state model: A practical framework for POMDP-based spoken dialogue management.  _Computer Speech & Language_, 24(2):150–174, 2010.