# Estimation of the last passage percolation constant in a charged complete
directed acyclic graph via perfect simulation

Sergey Foss, Takis Konstantopoulos, Bastien Mallein
and Sanjay Ramassamy

###### Abstract

Our object of study is the asymptotic growth of heaviest paths in a charged
(weighted with signed weights) complete directed acyclic graph. Edge charges
are i.i.d. random variables with common distribution $F$ supported on
$[-\infty,1]$ with essential supremum equal to $1$ (a charge of $-\infty$ is
understood as the absence of an edge). The asymptotic growth rate is a
constant that we denote by $C(F)$. Even in the simplest case where
$F=p\delta_{1}+(1-p)\delta_{-\infty}$, corresponding to the longest path in
the Barak-Erdős random graph, there is no closed-form expression for this
function, but good bounds do exist. In this paper we construct a Markovian
particle system that we call “Max Growth System” (MGS), and show how it is
related to the charged random graph. The MGS is a generalization of the
Infinite Bin Model that has been the object of study of a number of papers. We
then identify a random functional of the process that admits a stationary
version and whose expectation equals the unknown constant $C(F)$. Furthermore,
we construct an effective perfect simulation algorithm for this functional
which produces samples from the random functional.

Keywords and phrases. perfect simulation, coupling (from the past), random
graph, Markov process, stationarity, last passage percolation

AMS 2010 subject classification. Primary 82M31; secondary 60K15, 60G10, 05C80

##  1 Preliminaries

A Barak-Erdős random graph is a directed acyclic version of the standard
Erdős-Rényi graph [3]. We let $\mathbb{Z}^{+}$, the set of non-negative
integers, serve as the set of vertices. For each pair of vertices $i,j$ with
$i<j$, declare $(i,j)$ as an edge directed from $i$ to $j$ with probability
$p$, independently from any other pair. Then the maximum length $L_{n}$ of all
paths from vertex $0$ to $n$ satisfies a law of large numbers [23, 10]:
$\lim_{n\to\infty}L_{n}/n=C(p)$, a.s., where $C(p)$ is a certain
deterministic, increasing, analytic function of $0<p\leq 1$ [20, 21]. Owing to
the fact that such a graph appears as a model in various natural applications,
such as in computer systems [15, 16], in mathematical ecology [23, 24] and
others, information about $C(p)$ has been the subject of a number of papers
[23, 10, 7, 12, 20, 21, 11, 13].

Longest paths in Barak-Erdős graphs may be seen as a special case of the last
passage percolation model, which studies the growth rate of the length of the
longest path in a directed acyclic graph whose edges and/or vertices are
equipped with random weights. By contrast, first passage percolation is a
model studying shortest paths between two points in a graph (usually
undirected) whose edges and/or vertices are equipped with random weights. Both
models have mainly been studied in the case when the graphs are of the form
$\mathbb{Z}^{d}$, see e.g. [8, 1] and references therein. First passage
percolation has recently been investigated for Barak-Erdős graphs in the
sparse inhomogeneous setting [22].

We take interest in a generalization of Barak-Erdős graphs, considering the
last passage percolation problem on a complete directed acyclic graph, in
which each edge has a signed charge distributed according to an independent
copy of the random variable $w$ taking values in $\mathbb{R}\cup\\{-\infty\\}$
with a support bounded from above. We denote by $F$ the law of $w$, and let
$\\{w_{i,j},0\leq i<j\\}$, be a collection of i.i.d. copies of $w$. If $\pi$
is a path from $i$ to $j$, namely an increasing collection of vertices
$(i=i_{0},i_{1},\ldots,i_{\ell}=j)$ then its charge is defined as the sum of
the charges of its edges:
$w(\pi)=w_{i_{0},i_{1}}+\cdots+w_{i_{\ell-1},i_{\ell}}$, using the convention
that $-\infty+x=-\infty$ for all $x\in\mathbb{R}$. In other words, if a path
goes through an edge with charge $-\infty$, then the charge of the path is
$-\infty$. We define by convention the charge of a path consisting of a single
vertex as $0$.

We are concerned with the quantity

$W_{n}:=\sup\\{w(\pi):\,\pi\text{ is a path from $0$ to $n$}\\},$ (1)

the maximum charge of all paths between $0$ and $n$. Observe that
$(W_{n},n\geq 0)$ is a superadditive random sequence as direct computations
show that for all $n,m\geq 0$,

$W_{n+m}\geq W_{n}+W_{n,n+m},$

with $W_{n,n+m}=\sup\\{w(\pi):\,\pi\text{ is a path from $n$ to $n+m$}\\}$. As
$W_{n,n+m}$ is a copy of $W_{m}$ independent of $W_{n}$, by Kingman’s
subadditive ergodic theorem [17] we have

$\frac{W_{n}}{n}\to C(F),\quad\text{a.s. as $n\to\infty$}$ (2)

where $C(F)$ is a deterministic function of the law $F$. We refer to $C(F)$ as
the _last passage percolation constant_ of $F$. The objective of the present
article is to develop an approach to compute this constant through Monte Carlo
methods.

We denote by

$L=\inf\\{z\in\mathbb{R}:\mathbb{P}(w>z)=0\\}$

the _essential supremum_ of $F$ (this is the maximal point of the support of
the distribution $F$).

Consider momentarily the case $L\leq 0$. Then edge weights are nonpositive
a.s., and, since

$-W_{n}=-\inf\\{-w(\pi),\pi\text{ is a path from $0$ to $n$}\\},$

the problem is that of first passage percolation on the complete directed
graph. We claim that $C(F)=0$. Indeed, with $n\geq 2$, considering the 2-edge
path $(0,j,n)$,

$0\geq W_{n}\geq\sup_{1\leq j\leq n-1}(w_{0,j}+w_{j,n})\quad\text{a.s.},$

hence

$\displaystyle\mathbb{P}(W_{n}\geq 2(L-1))$
$\displaystyle\geq\mathbb{P}\left(w_{0,j}+w_{j,n}\geq 2(L-1)\text{ for all
}1\leq j\leq n-1\right)$ $\displaystyle\geq
1-(1-\mathbb{P}(w>L-1)^{2})^{n-1}.$

Thus, by the Borel-Cantelli lemma, $W_{n}/n\to 0$ in probability, and hence
$C(F)=0$, as claimed. It is not hard to see that $W_{n}$ itself converges
weakly to the random variable $\max(w,2L)$.

We only consider the case $L>0$ in the rest of the article. In this situation,
up to replacing $w$ by $w/L$, we will assume without loss of generality that
the essential supremum of $F$ is $1$. Hence, we work under the following
assumption for the distribution $F$ of $w$:

$\forall\varepsilon>0,\quad F([1-\varepsilon,1])>0\text{ and
}F((1,\infty))=0.$ (3)

The case $F=p\delta_{1}+(1-p)\delta_{-\infty}$ formally corresponds to a
Barak-Erdős graph as any edge with charge $-\infty$ can be ignored. Such a
graph was studied in [10] and a more general version of it in [7]. The
constant $C(p)$ mentioned earlier is, with an abuse of notation, the constant
$C(p\delta_{1}+(1-p)\delta_{-\infty})$. We know that $C(p)>0$ for all $p>0$
which means that eventually, any two vertices that are far apart enough are
connected by a path that has charge $>-\infty$. It was shown in [20, 21] that
the function $p\to C(p)$ is analytic on $(0,1]$ and a two-term asymptotic
expansion was given in the limit $p\to 0$ as well as the power series
expansion around $p=1$.

The case $F=p\delta_{1}+(1-p)\delta_{x}$, where $x\in(-\infty,1)$, was studied
in [13]. For this case, the quantity $C(F)$ was denoted by $C_{p}(x)$, a
differentiable function of $x\in(-\infty,1)\setminus I$, where $I$ is the
union of nonpositive rationals and of the inverses $1/n$, $n\geq 2$. Moreover
$C(p)$ is the decreasing limit of $C_{p}(x)$, as $x\to-\infty$. In the special
case when $x=0$, it was shown in [9] that $C_{p}(0)=1/\psi(1-p)$ where $\psi$
is a Ramanujan theta function.

Let $F$ be a distribution on $[-\infty,1]$ with essential supremum $1$.
Comparing $F$ with the distribution $p\delta_{1/2}+(1-p)\delta_{-\infty}$
where $p=F([\tfrac{1}{2},1])$, it is not hard to see that $C(F)>0$. The goal
of this paper is to _construct_ a random variable with expectation $C(F)$ that
can be perfectly simulated via an explicit algorithm. Perfect simulation of a
functional of a Markov chain in its “steady-state” is a technique that,
whenever applicable, avoids the bias introduced by standard MCMC (=Markov
Chain Monte Carlo) methods, in which one would approach $C(F)$ by a
realization of $W_{n}/n$ for $n$ large enough. The terminology and algorithm
were introduced in [25].

A survey can be found in [18]. Its relation to the so-called backwards-
coupling was studied in [14]. It belongs to the broader area of coupling
methods for stochastic recursions that may entirely lack the Markovian
property [4, 6, 10].

Our perfect simulation algorithm is based on the construction of a particle
system, that we call the Max Growth System (MGS) associated to the charged
complete directed graph. This particle system can be seen as an extension of
the Infinite Bin Model (IBM) [10, 20, 5, 21] arising in connection to the
Barak-Erdős graph. We mention en passant that the IBM is a particle system in
discrete time introduced in [10] but one which falls in a natural class of
similar particle systems, manifestations of which have appeared frequently in
the literature, e.g. in [2].

In Section 2, we first define the MGS with charge distribution $F$ and
describe some properties of its dynamics. In Section 3, we show that a certain
functional of the MGS is a Markov chain that admits a stationary version. In
Section 4, we pull the random variable mentioned above from the stationary
version, show that its expectation is $C(F)$ and describe a perfect simulation
algorithm. We conclude by suggesting further directions of research in Section
5.

##  2 The Max Growth System (MGS)

The Max Growth System is a particle system on $\mathbb{R}\cup\\{-\infty\\}$ in
which at every step a new atom is added to the process. This auxiliary
particle system is constructed in such a way that starting from a single
particle at position $0$, the $n$th particle in the system will be placed at
position $W_{n}$. The precise connection between the Max Growth System and the
last passage percolation model introduced in the previous section is given in
Lemma 2.

###  2.1 Deterministic dynamics of the MGS

We let $\mathcal{N}$ be the set of locally finite point measures on
$\mathbb{R}\cup\\{-\infty\\}$ with a finite maximal element, namely, measures
whose values are nonnegative integers and which are finite on every interval
of the form $[x,\infty)$. Another way to define $\mathcal{N}$ is as the set of
Radon measures $\nu$ on $\mathbb{R}\cup\\{-\infty\\}$ such that
$x\mapsto\nu([x,\infty))$ is a non-increasing function from $\mathbb{R}$ to
$\mathbb{Z}_{+}$. This will be the state space on which the MGS is defined.
Any such measure $\nu\in\mathcal{N}$ is specified by the nonincreasing
sequence $\nu_{1}\geq\nu_{2}\geq\cdots\geq-\infty$ of the locations of the
points (atoms) of $\nu$. This sequence may be finite or infinite. For example,
$\nu=2\delta_{0}+\delta_{-1.5}+3\delta_{-4}$ is equivalently represented by
the finite sequence $(0,0,-1.5,-4,-4,-4)$. We shall therefore think of any
$\nu\in\mathcal{N}$ either as a point measure $\nu=\sum_{k\geq
1}\delta_{\nu_{k}}$ or as a sequence $(\nu_{1},\nu_{2},\ldots)$. Note that the
zero measure $0$ is an element of $\mathcal{N}$ and corresponds to an empty
sequence of points. The total mass $\|\nu\|=\nu(\mathbb{R}\cup\\{-\infty\\})$
of $\nu$ is the number of its points (counted with multiplicity). We let
$\inf\nu:=\nu_{\|\nu\|}$ be the location of the last point of $\nu$ if
$\|\nu\|<\infty$. If $\|\nu\|=\infty$, we let $\inf\nu=-\infty$.

Let $w=(w_{1},w_{2},\ldots)$ be a sequence of elements of
$\mathbb{R}\cup\\{-\infty\\}$, such that $\sup_{k\geq 1}w_{k}\leq 1$. Let
$\mathcal{W}$ be the collection of such sequences. Given $\nu$ a non-zero
element of $\mathcal{N}$, define the quantity

$\mathfrak{m}(\nu,w):=\sup_{k\geq 1}(\nu_{k}+w_{k}),\quad\nu\neq 0,\quad
w\in\mathcal{W}.$

Here the supremum is taken either over all $k\geq 1$ if $\|\nu\|=\infty$ or
over $1\leq k\leq\|\nu\|$ if $\|\nu\|$ is finite. Observe that
$\mathfrak{m}(\delta_{0},w)=w_{1}$ for all $w\in\mathcal{W}$. The map
responsible for the dynamics of the MGS is defined by

$\Psi_{w}\nu:=\nu+\delta_{\mathfrak{m}(\nu,w)},$

that consists in adding at every step an atom in the process at a position
given by $\mathfrak{m}(\nu,w)$.

We will later employ a “coupling from the past” technique. To this end, it is
worth describing the MGS starting from an arbitrary point in time. Let
$(w(t),t\in\mathbb{Z})$ be a sequence of elements of $\mathcal{W}$, $\nu$ a
point measure in $\mathcal{N}$ and $T\in\mathbb{Z}$. The MGS starting from
$\nu$ at time $T$ is the process $(\nu(t),t\geq T)$ defined recursively by

$\nu(T)=\nu\quad\text{and}\quad\nu(t+1)=\Psi_{w(t+1)}\nu(t),\quad t\geq T.$

When $(w_{j}(t),j\in\mathbb{Z}^{+},t\in\mathbb{Z})$ is i.i.d. with law $F$, we
say that $(\nu(t),t\geq 0)$ is an MGS with charge distribution $F$. To
simplify notation, for all $s\leq t\in\mathbb{Z}$, we write

$\Psi_{w}^{s,t}=\Psi_{w(t)}\circ\Psi_{w(t-1)}\circ\cdots\circ\Psi_{w(s)},$

in which case we have $\nu(t)=\Psi_{w}^{T+1,t}\nu$ for all $t>T$.

To consider stationary versions of the MGS, we will sometimes need to work
with the particle system seen from the rightmost particle. We denote by
$\mathcal{N}_{0}$ the set of $\nu\in\mathcal{N}$ with $\nu_{1}=0$. For
$\nu\in\mathcal{N}$, we define its shift $\sigma\nu$ seen from the front by

$\int f(x)\,\mathrm{d}(\sigma\nu)(x):=\int f(x-\nu_{1})\,\mathrm{d}\nu,$

for all $\nu\in\mathcal{N}$ and all positive bounded functions $f$. Thus
$\sigma:\mathcal{N}\to\mathcal{N}$ and can be thought of as: “place the origin
at the position of the rightmost atom”. For example,
$\sigma(\delta_{a}+\delta_{b})=\delta_{0}+\delta_{-|a-b|}$. Observe that
$\sigma$ is a projection of $\mathcal{N}$ onto $\mathcal{N}_{0}$, which is
consistent with the definition of the MGS as

$\sigma\Psi_{w}=\sigma\Psi_{w}\sigma,$

for all sequences $w\in\mathcal{W}$. It is also worth mentioning that, for all
$\nu\in\mathcal{N}$, we have

$\mathfrak{m}(\sigma\nu,w)=\mathfrak{m}(\nu,w)-\nu_{1}.$ (4)

###  2.2 Decoupling properties of the MGS

The following lemma shows that if there is a large enough gap in between the
first and the second atom in the point measure $\nu$, and the sequence of
charges satisfies a “triangular” property, then the positions of the new
particles only depend on a finite number of charges.

######  Lemma 1 (Decoupling property).

Fix $\ell\in[0,1)$ and a positive integer $n$. Let $T\in\mathbb{Z}$ and
$(w(T+t),t\geq 1)$ be a sequence in $\mathcal{W}$. Let $\nu$ be a point
measure in $\mathcal{N}_{0}$ such that $\nu_{2}\leq-\ell$. We define the
sequences

$\nu(t)=\Psi_{w(t)}\nu(t-1)\quad\text{and}\quad\tilde{\nu}(t)=\Psi_{w(t)}\tilde{\nu}(t-1),\quad
t\geq T,$

with $\nu(T)=\nu$ and $\tilde{\nu}(T)=\delta_{0}$. For all $n\in\mathbb{N}$,
if

$\overline{w}(T;t):=\max\\{w_{1}(T+t),\ldots,w_{t}(T+t)\\}\geq
1-\ell\quad\text{ for all $1\leq t\leq n$,}$ (5)

then
$\mathfrak{m}(\nu(T+n-1),w(T+n))=\mathfrak{m}(\tilde{\nu}(T+n-1),w(T+n))$.

###### Proof.

It suffices to prove this statement for $T=0$. We prove, by induction, that

$\overline{w}(0;t)\geq 1-\ell\text{ for all }1\leq t\leq
n~{}\Rightarrow~{}{\nu}(n)_{|\mathbb{R}^{+}}=\tilde{\nu}(n)_{|\mathbb{R}^{+}}~{}\text{and}~{}\nu(n)(\mathbb{R}^{+})=n+1$
(6)

Assume first that $n=1$. In this case,

$\mathfrak{m}(\nu,w(1))=\max\big{(}w_{1}(1),\,\max_{j\geq
2}[\nu_{j}+w_{j}(1)]\big{)}.$

Since, by assumption, $\nu_{j}\leq\nu_{2}\leq-\ell$ for all $j\geq 2$, we have
$\nu_{j}+w_{j}(n)\leq-\ell+1$ for all $j\geq 2$. To prove (6) for $n=1$ we
must assume that $w_{1}(1)\geq 1-\ell$. But then $w_{1}(1)\geq\max_{j\geq
2}[\nu_{j}+w_{j}(1)]$ and so

$\mathfrak{m}(\nu,w(1))=w_{1}(1)=\mathfrak{m}(\tilde{\nu}(0),w(1)).$

Hence $\nu(1)=\nu+\delta_{w_{1}(1)}$ and, with $\tilde{\nu}(0)=\delta_{0}$,
$\tilde{\nu}(1)=\delta_{0}+\delta_{w_{1}(1)}$. Hence (6) holds for $n=1$.

Assume next that (6) holds for some $n\geq 2$. We prove that it also holds for
$n+1$. To do this, it suffices to assume that
$\nu(n)_{|\mathbb{R}^{+}}=\tilde{\nu}(n)_{|\mathbb{R}^{+}}$,
$\nu(n)(\mathbb{R}^{+})=n+1$, and $\overline{w}(0;t)\geq 1-\ell$ for all
$1\leq t\leq n+1$. In this case, we have

$\mathfrak{m}(\nu(n),w(n+1))=\max\left(\max_{j\leq
n+1}[\nu_{j}(n)+w_{j}(n+1)],\,\max_{j\geq n+2}[\nu_{j}(n)+w_{j}(n+1)]\right).$

But, for all $j\leq n+1$, $\nu_{j}(n)\geq 0$ and so

$\max_{j\leq n+1}[\nu_{j}(n)+w_{j}(n+1)]\geq 1-\ell.$

Taking into account the assumption $\nu_{2}\leq-\ell$, we have, for all $j\geq
n+2$, $\nu_{j}(n)\leq-\ell$ and so

$\max_{j\leq n+1}[\nu_{j}(n)+w_{j}(n+1)]\geq 1-\ell\geq\max_{j\geq
n+2}[\nu_{j}(n)+w_{j}(n+1)],$

which implies that

$\mathfrak{m}(\nu(n),w(n+1))=\mathfrak{m}(\tilde{\nu}(n),w(n+1))=\max_{j\leq
n+1}[\nu_{j}(n)+w_{j}(n+1)]\geq 1-\ell>0.$

The configuration $\nu(n+1)$ is thus obtained by adding a particle to $\nu(n)$
at a positive location. Since $\nu(n)=\tilde{\nu}(n)$ on $\mathbb{R}^{+}$ and
since the particle is added at the same location for both, we have
$\nu(n+1)=\tilde{\nu}(n+1)$ on $\mathbb{R}^{+}$. Clearly,
$\nu(n+1)(\mathbb{R}^{+})=\nu(n)(\mathbb{R}^{+})+1=n+2$, so (6) holds for
$n+1$. ∎

The above lemma allows us to describe a set of conditions on the sequences
$(w(t))$ so that the increments of $\nu$ and $\tilde{\nu}$ are algebraically
independent of $\nu$ and $\tilde{\nu}$.

######  Corollary 1.

Let $T\in\mathbb{Z}$, $n\in\mathbb{N}$, $\ell\in[0,1)$ and $(w(T+t),0\leq
t\leq n)$ a sequence such that

$w_{1}(T)\geq\ell\text{ and
}\min\\{\overline{w}(T;1),\overline{w}(T;2),\ldots,\overline{w}(T;n)\\}\geq
1-\ell.$ (7)

Let $\nu,\tilde{\nu}$ be two elements of $\mathcal{N}_{0}$ and define the
sequences

$\nu(t)=\Psi_{w(t)}\nu(t-1)\text{ and
}\tilde{\nu}(t)=\Psi_{w(t)}\tilde{\nu}(t-1),\quad t\geq T,$

with $\nu(T-1)=\nu$ and $\tilde{\nu}(T-1)=\tilde{\nu}$. Then
$\mathfrak{m}(\sigma\nu(t-1),w(t))=\mathfrak{m}(\sigma\tilde{\nu}(t-1),w(t))$
for all $T+1\leq t\leq T+n$.

In other words, the sequence
$\left(\mathfrak{m}(\sigma\nu(t-1),w(t)),~{}T+1\leq t\leq T+n\right)$ is
algebraically independent of $\nu(T-1)$ provided that $w$ satisfies (7).

###### Proof.

We observe that as $w_{1}(T)\geq\ell$, we have

$\mathfrak{m}(\nu(T-1),w(T))\geq\nu_{1}(T-1)+w_{1}(T)\geq\ell\quad\text{and}\quad\mathfrak{m}(\tilde{\nu}(T-1),w(T))\geq\ell.$

Therefore, the second largest atoms of $\sigma\Psi_{w(T)}\nu(T-1)$ and
$\sigma\Psi_{w(T)}\tilde{\nu}(T-1)$ are both smaller than $-\ell$, hence by
(4) we can apply Lemma 1, which completes the proof. ∎

###  2.3 The MGS derived from the charged complete directed graph

Consider the charged complete directed graph with i.i.d. edge charges
$\\{w_{i,j},0\leq i<j\\}$ of law $F$, a collection of i.i.d. random variables
in $\mathbb{R}\cup\\{-\infty\\}$ with common law $F$ satisfying assumption
(3). For all $n\in\mathbb{N}$, we write $W_{n}$ for the length of the longest
path between $0$ and $n$. We observe that $(W_{n},n\geq 1)$ can be coupled
with the MGS with charge distribution $F$.

######  Lemma 2.

Let $(\nu(t),t\geq 0)$ be an MGS with charge distribution $F$ such that
$\nu(0)=\delta_{0}$ and let $(W_{t},t\geq 0)$ as defined in (1). We have the
following equality in distribution:

$(\nu(t),t\geq 0)\overset{(d)}{=}\left(\sum_{j=0}^{t}\delta_{W_{j}},t\geq
0\right).$ (8)

###### Proof.

By definition, we have $\nu(0)=\delta_{0}=\delta_{W_{0}}$, using that the path
of length $0$ between $1$ and $1$ has mass $0$. Let $t_{0}\in\mathbb{Z}^{+}$
and assume that we can construct a coupling between $\nu$ and $W$ such that
$(\nu(t),t\leq t_{0})=(\sum_{j=0}^{t}\delta_{W_{j}},t\leq t_{0})$ a.s.
Conditionally on this coupling, let $(w_{j,t_{0}+1},j\in\mathbb{Z}^{+})$ and
$(w_{j}(t_{0}+1),j\in\mathbb{Z}^{+})$ be independent families of i.i.d. random
variables with law $F$. By (1), decomposing all paths $\pi$ ending at
$t_{0}+1$ according to their last step, we have

$W_{t_{0}+1}=\max\\{W_{j}+w_{j,t_{0}+1},0\leq j\leq t_{0}\\}\\\
\overset{(d)}{=}\max\\{\nu_{j}(t_{0})+w_{j}(t_{0}+1),1\leq j\leq
t_{0}+1\\}=\mathfrak{m}(\nu(t_{0}),w(t_{0}+1)),$

therefore $\nu(t_{0}+1)\overset{(d)}{=}\sum_{j=0}^{t_{0}+1}\delta_{W_{j}}$. As
a result, we can couple the two sequences of random variables in such a way
that the above equality holds almost surely.

Hence, by recursion, there exists a coupling between the MGS and the last
passage percolation problem such that (8) holds for all times. ∎

A noteworthy observation is that the increments of $W$ are the same as the
relative increments of the MGS. More precisely, defining

$M_{n}=\max_{0\leq k\leq n}W_{k}=\sup_{0\leq k\leq n}\left\\{w(\pi),\pi\text{
path from $0$ to $k$}\right\\},$ (9)

the increments of the sequence $(M_{n},n\geq 0)$ can be connected to the
relative increments of the MGS.

######  Corollary 2.

Under the foregoing assumptions,

$(M_{n}-M_{n-1},n\geq
1)\overset{(d)}{=}(\mathfrak{m}(\sigma\nu(n-1),w(n))^{+},n\geq 1).$

###### Proof.

In the proof above, we established a coupling between $(\nu(n),n\geq 0)$ and
$(W_{k},k\geq 0)$. Under this coupling, for $n\in\mathbb{Z}^{+}$, we have

$W_{n}=\mathfrak{m}(\nu(n-1),w(n))\quad\text{and}\quad M_{n}=\nu_{1}(n).$

As a result, under this coupling, we have

$M_{n}-M_{n-1}=(W_{n}-M_{n-1})^{+}=(\mathfrak{m}(\nu(n-1),w(n))-\nu_{1}(n-1))^{+}=\mathfrak{m}(\sigma\nu(n-1),w(n))^{+}.\qed$

##  3 Stationarity via coupling

We recall that our aim is to compute the quantity $C(F)$ defined by

$C(F):=\lim_{n\to\infty}\frac{W_{n}}{n}\quad\text{a.s.}$

As $F$ has a finite essential supremum, it holds that
$\int_{0}^{\infty}xF(\mathrm{d}x)<\infty$. Therefore, by [12], it is known
that

$C(F)=\lim_{n\to\infty}\frac{M_{n}}{n}\quad\text{a.s. and in }L^{1}.$

Thus, if $F$ satisfies (3), we have

$C(F)=\lim_{n\to\infty}\frac{\mathbb{E}(M_{n})}{n}.$

Using Corollary 2, we remark that for all $n\in\mathbb{N}$,

$\frac{M_{n}}{n}=\frac{1}{n}\sum_{j=1}^{n}(M_{j}-M_{j-1})=\frac{1}{n}\sum_{j=1}^{n}\mathfrak{m}(\sigma\nu(j-1),w(j))^{+}.$

We show in this section that $\left(\mathfrak{m}(\sigma\nu(n-1),w(n)),n\geq
1\right)$ admits a stationary version (where we recall that $\sigma\nu$ is the
point measure shifted so that its rightmost element is at position $0$). Since
the process $\left(\mathfrak{m}(\sigma\nu(n-1),w(n)),n\geq 1\right)$ is not
Markovian, the term “stationary version” should be used with caution. For us,
it means that it couples with a stationary process in finite time, as in the
statement of Theorem 1 below. Then, letting $\overline{\mathfrak{m}}$ be the
limit in distribution of $\mathfrak{m}(\sigma\nu(n-1),w(n))$ as $n\to\infty$,
we have

$C(F)=\lim_{n\to\infty}\frac{1}{n}\sum_{j=1}^{n}\mathbb{E}\left(\mathfrak{m}(\sigma\nu(j-1),w(j))^{+}\right)=\mathbb{E}(\overline{\mathfrak{m}}^{+}),$

as $\mathfrak{m}(\sigma\nu,w)^{+}\in[0,1]$ a.s. In the next section, we
introduce the perfect simulation algorithm, which consists in giving a
realization of $\overline{\mathfrak{m}}$ without constructing the limit in
distribution of $\sigma\nu(n)$ as $n\to\infty$. In the special case of Barak-
Erdős graphs, a simpler case of a perfect simulation algorithm was explained
in [10].

######  Theorem 1.

Suppose that $F$ is a distribution satisfying (3). Let $\ell\in[0,1)$ be such
that $p:=F([1-\ell,1])\in(0,1]$. Let $w=(w_{i},i\in\mathbb{N})$ be i.i.d.
random variables with law $F$ and $(w(t),t\in\mathbb{Z})$ i.i.d. copies of
$w$. Given $\nu(0)\in\mathcal{N}$, we define the MGS by

$\nu(t+1)=\Phi_{w(t+1)}\nu(t),\quad\text{ for all $t\geq 0$}.$

There exists a stationary process
$(\overline{\mathfrak{m}}(t),t\in\mathbb{Z})$ such that

$\mathfrak{m}(\sigma\nu(t-1),w(t))=\overline{\mathfrak{m}}(t)\quad\text{a.s.
for $t$ large enough}.$

In particular $\mathbb{E}(\overline{\mathfrak{m}}(0)^{+})=C(F)$.

###### Proof.

For $T\in\mathbb{Z}$ and $t\in\mathbb{N}$, we recall the notation
$\overline{w}(T;t)=\max\\{w_{1}(T+t),\ldots,w_{t}(T+t)\\}$ from Lemma 1. We
introduce the event

$R_{k}:=\bigcap\limits_{j=1}^{\infty}\\{w_{1}(k)\geq\ell,\overline{w}(k;j)\geq
1-\ell\\}.$

It is clear from its definition that $(R_{k},k\in\mathbb{Z})$ is a stationary
sequence of events with

$\mathbb{P}(R_{k})=\mathbb{P}(R_{0})=F([\ell,1])\prod_{j=1}^{\infty}(1-(1-p)^{j})>0.$

Consider the stationary random set $J:=\\{k\in\mathbb{Z}:R_{k}\text{
holds}\\}$. Since $\mathbb{P}(R_{k})>0$, we have, by ergodicity (more
specifically by the Poincaré recurrence theorem), $\inf J=-\infty$ and $\sup
J=\infty$ a.s. We enumerate the elements of $J$ by

$\cdots<T_{-1}<T_{0}\leq 0<T_{1}<T_{2}<\cdots$

We define

$\tilde{\nu}(t):=\sigma\sum_{i\in\mathbb{Z}}\mathbbm{1}_{\\{T_{i}<t\leq
T_{i+1}\\}}\Psi^{T_{i},t}_{w}\delta_{0},\quad t\in\mathbb{Z}.$

It is clear from its definition that $(\tilde{\nu}(t),t\in\mathbb{Z})$ is
stationary, as $(w(t),t\in\mathbb{Z})$ is a stationary sequence, and
$(\tilde{\nu}(T_{i}+1),i\in\mathbb{Z})$ are i.i.d. elements of
$\mathcal{N}_{0}$. Next, we define

$\overline{\mathfrak{m}}(t)=\mathfrak{m}(\tilde{\nu}(t-1),w(t)),\quad
t\in\mathbb{Z},$

which is again a stationary sequence.

By Corollary 1, we observe that for all $t\geq T_{1}+1$, the quantity
$\mathfrak{m}(\sigma\nu(t-1),w(t))$ does not algebraically depend on
$\nu(T_{1}-1)$. Hence, we have
$\mathfrak{m}(\sigma\nu(t-1),w(t))=\overline{\mathfrak{m}}(t)$, using that
$\overline{\mathfrak{m}}(t)$ is the same quantity for the MGS started from
$\delta_{0}$ at time $T_{1}-1$. As $T_{1}<\infty$ a.s. this completes the
proof of the first part of the theorem.

Next, using that

$C(F)=\lim_{n\to\infty}\frac{\mathbb{E}(M_{n})}{n}=\lim_{n\to\infty}\frac{1}{n}\sum_{j=1}^{n}\mathbb{E}(\mathfrak{m}(\sigma\nu(j-1),w(j))^{+})\quad\text{a.s.},$

and using the eventual equality between $\mathfrak{m}(\sigma\nu(t-1),w(t))$
and $\overline{\mathfrak{m}}(t)$, we have

$C(F)=\lim_{n\to\infty}\frac{1}{n}\sum_{j=1}^{n}\mathbb{E}(\overline{\mathfrak{m}}(j)^{+})=\mathbb{E}(\overline{\mathfrak{m}}(0)^{+}),$

by stationarity and ergodicity of the sequence. ∎

######  Remark 1.

The random times $T_{i}$ split the process into independent and identically
distributed pieces (thereby making the process strictly regenerative) yielding
a number of limiting results including a (functional) central limit theorem.
In terms of the last passage percolation model, the $T_{i}$ are the locations
of points through which every longest path must pass. Thus the stationary last
passage percolation model admits bi-infinite longest paths, and any longest
path in a finite graph grown from a single initial vertex will eventually
coalesce with some bi-infinite longest path.

##  4 Perfect simulation

The formula (2) for $C(F)$ suggests a straightforward method for estimating
$C(F)$: starting from $\nu(0)=0$, generate iteratively
$\nu(1),\nu(2),\ldots,\nu(n)$, and take $\nu_{1}(n)$ for an estimation of
$C(F)$. This standard (so-called MCMC) method introduces a bias. Indeed,
$\mathbb{E}(\nu_{1}(n))/n$ is not equal to $C(F)$, but merely converges to
that constant.

To eliminate this bias, we produce an algorithm that constructs the variable
$\overline{\mathfrak{m}}(0)$, whose distribution is unknown. Then, by standard
Monte Carlo method, an unbiased estimation of $C(F)$ can be constructed. This
is done in this case by using the construction described in the proof of
Theorem 1.

This algorithm is a development of a similar construction for functionals of
stochastic recursions in [10] and is based on the ideas of so-called “backward
coupling”, see [14]. It is close in spirit to the coupling-from-the-past
method for Markov chains [25] and to the perfect simulation construction for
processes with “long memory” [6]. Note that the algorithm from [25] is
applicable to either finite Markov chains or ordered monotone Markov chains
possessing a unique minimal state and a unique maximal state, so it cannot be
applied to our case.

######  Theorem 2 (Perfect simulation).

Define

$T^{*}:=\sup\\{t\leq-1:\,{w}_{1}(t)\geq\ell,\min_{1\leq
j\leq|t|}\overline{w}(t;j)\geq 1-\ell\\}.$

Then $|T^{*}|<\infty$ a.s., and

$\overline{\mathfrak{m}}(0)=\mathfrak{m}\left(\sigma\Psi^{T^{*},-1}_{w}\delta_{0};w(0)\right)\text{
a.s.}$

###### Proof.

We recall that $(T_{-j},j\in\mathbb{N})$ are the negative elements of the
random set $J$, with $T_{-1}>-\infty$. We remark that

$w_{1}(T_{-1})\geq\ell,\quad\overline{w}(T_{-1};j)\geq 1-\ell\text{ for all
$j>0$},$

therefore $T_{*}\geq T_{-1}$, proving its finiteness.

Moreover, since

$w_{1}(T^{*})\geq\ell,\min_{1\leq j\leq|t|}\overline{w}(T^{*};j)\geq 1-\ell,$

by Corollary 1, the quantity
$\mathfrak{m}\left(\sigma\Psi^{T^{*},-1}_{w}\nu;w(0)\right)$ does not
algebraically depend on the value of $\nu\in\mathcal{N}_{0}$. As a result, it
is equal to $\overline{\mathfrak{m}}(0)$, defined as
$\mathfrak{m}\left(\sigma\Psi^{T^{*},-1}_{w}\delta_{0};w(0)\right)$. ∎

######  Remark 2.

If the essential supremum $L$ of $F$ is infinite, then the perfect simulation
algorithm we defined cannot apply. Indeed, in this situation, Lemma 1 does not
apply and we could not find an event depending on a finite number of charges
such that an analogue of this lemma would hold. When $L=\infty$, even if the
tail of $F$ decays fast enough, we would still need to look at infinitely many
values of $w_{j}(1)$ to increment just the first time step of the MGS, making
it impossible to hope for a perfect simulation algorithm which ends in finite
time for any starting configuration.

Fix $t=0$ and $J=1$;

Generate the variable $w_{1}(0)$;

Fix Stopping = False;

while  _Stopping = False_ do

while  _$\max_{1\leq j\leq J}w_{j}(t) <1-\ell$_ do

Increase $J$ by $1$;

Generate the variable $w_{J}(t)$;

while  _$J >1$_ do

Decrease $J$ by $1$ and $t$ by $1$;

Generate $w_{1}(t),\ldots w_{J}(t)$;

while  _$\max_{1\leq j\leq J}w_{j}(t) <1-\ell$_ do

Increase $J$ by $1$;

Generate the variable $w_{J}(t)$;

Decrease $t$ by $1$;

Generate $w_{1}(t)$;

Fix Stopping = $\\{w_{1}(t)\geq\ell\\}$;

Fix $\nu=\delta_{0}$;

for  _$s$ from $t+1$ to $-1$_ do

Generate the variables $w_{1}(s),\ldots,w_{||\nu||}(s)$ ;

Set $\mathfrak{m}=\max\\{\nu_{j}+w_{j}(s)\text{ for $1\leq j\leq\|\nu\|$}\\}$;

Add $\delta_{\mathfrak{m}}$ to $\nu$;

Set $\mathfrak{m}=\max\\{\nu_{j}+w_{j}(0)\text{ for }1\leq j\leq\|\nu\|\\}$;

Return: $\mathfrak{m}-\nu_{1}$;

Algorithm 1 Construction of a variable of law $\overline{\mathfrak{m}}(0)$.

### The perfect simulation algorithm

We now describe more precisely the perfect simulation algorithm. Let $F$ be a
probability distribution satisfying (3), we fix $\ell\in[0,1)$ such that
$F([1-\ell,1])\in(0,1)$. The algorithm requires the construction of an array
of i.i.d. random variables with common distribution $F$ until the random
variable $T^{*}$ can be constructed.

To construct $T^{*}$ as well as $\overline{\mathfrak{m}}(0)$ from the sequence
$\\{w_{j}(t),j\in\mathbb{N},t\in\mathbb{Z}\\}$, one only needs to consider
a.s. finitely many elements of this set, as $\\{T^{*}=t\\}$ is a measurable
function of

$\\{w_{1}(t)\\}\cup\\{w_{j}(t+k),1\leq j\leq k\leq|t|\\}$

and $\overline{\mathfrak{m}}(0)$ is a measurable function of

$\\{w_{1}(T^{*})\\}\cup\\{w_{j}(T^{*}+k),1\leq j\leq k\leq|T^{*}|\\}.$

Therefore, we can explore triangular arrays of the form

$\\{w_{1}(t)\\}\cup\\{w_{j}(t+k),1\leq j\leq k\leq|t|\\},$

progressively decreasing $t$ until time $T^{*}$ is detected. Once this random
variable is known, we construct the random variable
$\overline{\mathfrak{m}}(0)$ using the procedure described in Theorem 1 from
the previously discovered random variables. A possible implementation is
described in Algorithm 1. We show a graphical representation of a run of
Algorithm 1 in Figure 1.

Figure 1: Illustration of the execution of Algorithm 1 on an example, in the
case where $\ell<1-\ell$. The variables sampled until the Boolean variable
Stopping becomes True are pictured by black/white squares and disks. One
searches for the first time $T^{*}$ such that every line of index $T^{*}+1\leq
t\leq 0$ has at least one black disk between columns $1$ and $t-T^{*}$ and
such that there is a black square in position $(T^{*},1)$. The full triangular
array of variables used in the construction of $\nu$ is enclosed by a red
boundary.

We observe that this algorithm has a complexity of $(T^{*})^{2}$, as it is the
number of steps needed to generate the variable $\overline{\mathfrak{m}}(0)$.
It is worth noting that $-T^{*}$ can be constructed as the first hitting time
of $0$ of the Markov chain $(X_{n})$ with initial state

$X_{0}=\min\\{j\geq 1,w_{j}(0)\geq 1-\ell\\}$

and with transition probabilities defined for all $j\geq 2$ and $i\geq j$ by

$P(j,j-1)=1-(1-p)^{j-1}\text{ and }P(j,i)=p(1-p)^{i-1}$

where $p=\mathbb{P}(w_{1}(0)\geq 1-\ell)$, with

$P(1,0)=\mathbb{P}(w_{1}(0)\geq\ell),\quad\mathbb{P}(1,1)=\mathbb{P}(1-\ell\leq
w_{1}(0)<\ell),\quad P(1,j)=p(1-p)^{j-1}\text{ for }j\geq 2.$

The quantity $X_{n}$ corresponds to the value of the variable $J$ at the end
of the period when $t=-n$ in Algorithm 1. In the example shown in Figure 1, we
have

$(X_{0},X_{-1},X_{-2},X_{-3},X_{-4},X_{-5},X_{-6})=(4,3,2,1,2,1,0).$

Note that $T^{*}$ has exponential tails.

The choice of the parameter $\ell$ may have an important effect on the
behaviour of the average complexity $\mathbb{E}((T^{*})^{2})$ of the
algorithm. We plotted $\ell\mapsto\mathbb{E}((T^{*})^{2})$ in Figure 2, when
the charge distribution is given by $F(\mathrm{d}x)=\mathbbm{1}_{\\{x\leq
1\\}}e^{x-1}\mathrm{d}x$. Additionally, as $p\to 0$, the quantity
$\mathbb{E}((T^{*})^{2})$ grows to $\infty$. We estimated
$\mathbb{E}((T^{*})^{2})$ for $F=p\delta_{1}+(1-p)\delta_{-\infty}$ and
plotted this quantity as a function of $p$ in Figure 3.

Figure 2: Dependency in the parameter $\ell$ of the complexity of Algorithm 1
with a charge distribution $F(\mathrm{d}x)=\mathbbm{1}_{\\{x\leq
1\\}}e^{x-1}\mathrm{d}x$. The figure was obtained with a Monte Carlo
simulation of $N=10^{4}$ copies of $T^{*}$ for $100$ different values of
$\ell$. For this charge distribution, the Monte Carlo simulations give
$C(F)=0.4432\pm 0.0006$.

We observe in Figure 2 that different choices of the value $\ell$ can have a
dramatic impact on the efficiency of Algorithm 1. Choosing a value $\ell$ too
small has the effect of making the first appearance of a triangular event too
late. On the other hand, if $\ell$ is too big then with high probability, one
will have $w_{1}(T)\leq\ell$, and thus the first “successful” triangular event
will appear much later. For the distribution $F$ we chose, it appears that an
optimal choice of $\ell$ seems to be around $\ell=0.7$, which balances between
these two extremes.

Figure 3: Dependency in the parameter $p$ of the complexity of Algorithm 1
applied to the detection of the longest path in the Barak-Erdős graph with
parameter $p$. Figure obtained through Monte Carlo simulation of $N=10^{5}$
copies of $T^{*}$ for 120 different values of $p$.

We observe in Figure 3 that if $F$ puts a large mass on the negative half-
line, the complexity of Algorithm 1 can become quite large. The function
$p\mapsto\mathbb{E}(T^{*})^{2}$ grows at least exponentially in $1/p$ as $p\to
0$ in the Barak-Erdős graph, but we were not able to obtain a good estimate of
this rate of increase.

##  5 Further directions of research

In this article we considered last passage percolation on the directed
complete graph, which has a total order on its vertex set. One extension of
this would be to construct a perfect simulation algorithm for so-called
directed slab graphs [7] where the set of vertices is only partially ordered.
Another possible extension would be to add i.i.d. vertex weights with a
distribution that has a finite essential supremum. In both cases, as well as
in the setting considered in this paper, one should be able to obtain a
perfect simulation algorithm if one replaces the i.i.d. weights by more
general stochastic recursions with stationary drivers, as was considered in
[10].

As discussed in the previous section, the complexity of our perfect simulation
algorithm may dramatically vary with $\ell$. In the case of
$F(\mathrm{d}x)=\mathbbm{1}_{\\{x\leq 1\\}}e^{x-1}\mathrm{d}x$ presented in
Figure 2, there seems to be a unique optimal choice for $\hat{\ell}$ around
$0.7$. It would be interesting to find some classes of distributions $F$ for
which one has good bounds on the optimal value $\hat{\ell}$.

Yet another research direction would be the estimation of the constants
appearing for last passage percolation on a 2-dimensional version of the
Barak-Erdős directed graph on the set $\mathbb{N}\times\mathbb{N}$ and whose
edges are as follows: if
$u=(u_{1},u_{2}),v=(v_{1},v_{2})\in\mathbb{N}\times\mathbb{N}$ are two
vertices such that $i_{1}\leq j_{1}$, $i_{2}\leq j_{2}$, then declare the pair
$(u,v)$ as an edge directed from $u$ to $v$ with probability $p$,
independently over all such pairs. Then maximum length $L_{n}$ of all paths
from $(1,1)$ to $(n,n^{a})$, for a certain $a>0$, rescaled appropriately,
converges weakly [19] to a random variable having a Tracy-Widom distribution
depending on two parameters. The estimation of these parameters is an open
problem.

## Acknowledgements

We thank the referees for suggestions to improve the exposition.

SF was partially supported by the RFBR collaborative grant 19-51-15001 and TK,
BM and SR were partially supported by the CNRS PRC collaborative grant
CNRS-193-382 with the common title “Asymptotic and analytic properties of
stochastic ordered graphs and infinite bin models”.

## References

  * [1] A. Auffinger, M. Damron and J. Hanson. 50 years of first-passage percolation. University Lecture Series, American Mathematical Society, Providence, RI, 68, v+161 pp, 2017.
  * [2] D. Aldous and J. Pitman. The asymptotic speed and shape of a particle system. In: Probability, statistics and analysis, London Math. Soc. Lecture Note Ser. 79, 1–23, 1983.
  * [3] A.B. Barak and P. Erdős. On the maximal number of strongly independent vertices in a random acyclic directed graph. SIAM J. Algebraic Discrete Methods 5, 508–514, 1984.
  * [4] A.A. Borovkov and S.G. Foss. Stochastically recursive sequences and their generalizations. Siberian Adv. Math., 2(1):16–81, 1992.
  * [5] K. Chernysh and S. Ramassamy. Coupling any number of balls in the infinite-bin model. J. Appl. Probab. 54 no. 2, 540–549, 2017.
  * [6] F. Comets, R. Fernández and P.A. Ferrari. Processes with long memory: regenerative construction and perfect simulation. Ann. Appl. Probab. 12, 921–943, 2002.
  * [7] D. Denisov, S. Foss and T. Konstantopoulos. Limit theorems for a random directed slab graph. Ann. Appl. Probab., 22, 702–733, 2012.
  * [8] D. Dauvergne, M. Nica and B. Virág. RSK in last passage percolation: a unified approach. Probab. Surv., 19, 65–112, 2022.
  * [9] K. Dutta. On limit constants in last passage percolation in transitive tournaments. arXiv:2005.09922 [math.PR]
  * [10] S. Foss and T. Konstantopoulos. Extended renovation theory and limit theorems for stochastic ordered graphs. Markov Process. Related Fields 9, no. 3, 413–468, 2003.
  * [11] S. Foss and T. Konstantopoulos. Limiting properties of random graph models with vertex and edge weights. J. Stat. Phys. 173, 626–643, 2018.
  * [12] S. Foss, J. Martin and P. Schmitt. Long-range last-passage percolation on the line. Ann. Appl. Probab. 24, no. 1, 198–234, 2014.
  * [13] S. Foss, T. Konstantopoulos and A. Pyatkin. Probabilistic and analytical properties of the last passage percolation constant in a weighted random directed graph. Ann. Appl. Probab., to appear. arXiv:2006.01727v3 [math.PR]
  * [14] S. Foss and R.L. Tweedie. Perfect simulation and backward coupling. Comm. Statist. Stochastic Models 14, no. 1–2, 187–203, 1998.
  * [15] E. Gelenbe, R. Nelson, T. Philips and A. Tantawi. An approximation of the processing time for a random graph model of parallel computation. Proceedings of 1986 ACM Fall Joint Computer Conference, Los Alamitos, CA, 691–697, 1986.
  * [16] M. Isopi and C.M. Newman. Speed of parallel processing for random task graphs. Comm. Pure Appl. Math. 47, 261–276, 1994.
  * [17] J. F. C. Kingman.  Subadditive ergodic theory.  Ann. Probab., 1:883–909, 1973.
  * [18] W.S. Kendall. Notes on perfect simulation. In: Markov chain Monte Carlo, Lect. Notes Ser. Inst. Math. Sci. Natl. Univ. Singap. 7, 93–146, 2005.
  * [19] T. Konstantopoulos and K. Trinajstić. Convergence to the Tracy-Widom distribution for longest paths in a directed random graph. ALEA, Lat. Am. J. Probab. Math. Stat. 10, no. 2, 711–730, 2013.
  * [20] B. Mallein and S. Ramassamy. Barak-Erdős graphs and the infinite-bin model. Ann. Inst. Henri Poincaré Probab. Stat. 57, no. 4, 1940–1967, 2021.
  * [21] B. Mallein and S. Ramassamy. Two-sided infinite-bin models and analyticity for Barak-Erdős graphs. Bernoulli 25, no. 4B, 3479–3495, 2019.
  * [22] B. Mallein and P. Tesemnikov. On the length of the shortest path in a sparse Barak-Erdős graph. Statist. Probab. Lett. 190, Paper No. 109634, 2022.
  * [23] C.M. Newman. Chain lengths in certain random directed graphs. Random Structures Algorithms 3, no. 3, 243–253, 1992.
  * [24] C.M. Newman and J.E. Cohen. A stochastic theory of community food webs IV: theory of food chains in large webs. Proc. Roy. Soc. Lond. B: Biol. Sci. 228, 355–377, 1986.
  * [25] J.G. Propp and D.B. Wilson. Exact sampling with coupled Markov chains and applications to statistical mechanics. Random Structures Algorithms 9, 223–252, 1996.

Sergey Foss,
Sobolev Institute of Mathematics, Siberian Branch of the Russian Academy of
Sciences, Novosibirsk, Russia
Novosibirsk State University, Novosibirsk, Russia
School of Mathematical and Computer Sciences, Heriot-Watt University,
Edinburgh, UK
E-mail address<EMAIL_ADDRESS>
Takis Konstantopoulos,
Department of Mathematical Sciences, University of Liverpool, Liverpool, UK
E-mail address<EMAIL_ADDRESS>
Bastien Mallein,
Université Sorbonne Paris Nord, LAGA, UMR 7539, F- 93430, Villetaneuse, France
E-mail address<EMAIL_ADDRESS>
Sanjay Ramassamy,
Université Paris-Saclay, CNRS, CEA, Institut de physique théorique, 91191 Gif-
sur-Yvette, France
E-mail address<EMAIL_ADDRESS>


# FreaAI: Automated extraction of data slices to test machine learning models

Samuel Ackerman, Orna Raz, Marcel Zalmanovici

(IBM Research, Haifa, Israel
)

###### Abstract

Machine learning (ML) solutions are prevalent. However, many challenges exist
in making these solutions business-grade. One major challenge is to ensure
that the ML solution provides its expected business value. In order to do
that, one has to bridge the gap between the way ML model performance is
measured and the solution requirements. In previous work (Barash et al. 2019)
we demonstrated the effectiveness of utilizing feature models in bridging this
gap. Whereas ML performance metrics, such as the accuracy or F1-score of a
classifier, typically measure the average ML performance, feature models shed
light on explainable data slices that are too far from that average, and
therefore might indicate unsatisfied requirements. For example, the overall
accuracy of a bank text terms classifier may be very high, say $98\%\pm 2\%$,
yet it might perform poorly for terms that include short descriptions and
originate from commercial accounts. A business requirement, which may be
implicit in the training data, may be to perform well regardless of the type
of account and length of the description. Therefore, the under-performing data
slice that includes short descriptions and commercial accounts suggests
poorly-met requirements. In this paper we show the feasibility of
automatically extracting feature models that result in explainable data slices
over which the ML solution under-performs. Our novel technique, IBM FreaAI aka
FreaAI, extracts such slices from structured ML test data or any other labeled
data. We demonstrate that FreaAI can automatically produce explainable and
statistically-significant data slices over seven open datasets.

## Introduction

Software systems that utilize ML are becoming more and more prevalent. The
quality of ML-based solutions greatly depends on the data used to train a
solution. Common approaches for assessing the quality of trained ML solutions
(we will refer to to these as ML) use the average performance of the ML
according to some metric, such as accuracy or F1-score for a classifier and
$R^{2}$ or root mean square error (RMSE) for regression, over unseen yet
representative data, such as that of a test set. Of course, even if these
average measures are satisfactory, there are likely to exist data records for
which the ML performance is far below the average. It is important to
understand this behavior as these far-away records might be indicative of
unmet requirements. For example, the third row in Table 8 indicates that the
ML model trained over loan approval data (the Adult dataset) under-performs
for people in the ages of 33 to 64. While the overall accuracy is roughly 85%,
the accuracy for people in this age group is roughly 80%. Probably this fails
to meet a requirement of providing similar ML performance regardless of age.
Another example in that table is under-performance of roughly 67% accuracy for
people who work 40 to 43 hours per week and belong to workclass category 5. It
is probably required that the solution be insensitive to the number of hours a
person works and their workclass. These are examples of ‘data slices’. The
first is a single feature data slice and the second is a 2-way feature
interaction data slice.

In previous work (Barash et al. 2019) we showed that data slices, defined by a
Combinatorial Testing (CT) feature model, result in identifying far-away
records in a way that a user can understand and act upon. In that work the
feature models were defined manually, and included various abstractions of the
ML input data as well as additional metadata. In this work we concentrate on
automatically creating feature models for one simple yet useful type of
abstraction—that of binning input data feature values. The next section
provides background about feature models.

While searching over single categorical feature values to find interesting
slices is straightforward, searching over continuous-valued features and over
interactions of multiple features is challenging. An exhaustive search is
infeasible as the search space is exponential. Moreover, reporting such
records or data slices should be done such that the user can take corrective
actions. This means that the user needs to understand what defines a data
slice. FreaAI is our novel technology that implements a set of heuristics to
efficiently suggest data slices that are, by construction, explainable,
correct, and statistically significant. ‘Explainable’ means that data slices
are described solely as ranges or single values and their unions. FreaAI is
named after the Norse goddess Frea (spelling varies) who, among other traits,
can foresee the future, but does not share this knowledge with humans. Our
technology suggests data slices that are correct, yet the results may be
incomplete in that not all possible under-performing slices may be found.

The contribution of this paper is a set of heuristics for automatically
defining feature models, such that each slice in that model is indicative of
ML under-performance, is statistically significant, and can be easily
understood by a human. We validate these heuristics by analyzing seven open
datasets. FreaAI generates data slices that contain significantly more mis-
predicted data records than would be expected by randomly drawing the same
number of data records.

We next list related work and provide background from our previous work about
feature models and data slices. The methodology of heuristically and
automatically suggesting data slices and the definition of our requirements
from data slices are then presented. The experimental results of running our
technology on seven open data-sets follow, and we conclude with a discussion
and brief summary of the contributions of this work.

## Related work and background

Several works capture challenges and best practices in developing ML
solutions, e.g., (Sculley et al. 2014; Breck et al. 2016; Zinkevich 2018), or
suggest organizational work-flows for building successful AI projects (Ng
2018).

We summarize the most relevant background and results from our previous work
(Barash et al. 2019). The work we report on in here is motivated by that work.
In a nutshell, we suggested the use of CT modeling methodology for validating
ML solutions. Building on that work, here we report on FreaAI, a novel
technology for automatically extracting feature models, which are an instance
of CT models, out of ML data, implementing an abstraction of automated binning
according to ML performance.

Following are more details on our previous work. We demonstrated that the
Combinatorial modeling methodology is a viable approach for gaining actionable
insights on the quality of ML solutions. Applying the approach results in
identifying coverage gaps indicating uncovered requirements. The approach also
identifies under-performing slices, suggesting areas where the current ML
model cannot be relied upon. When used in the development stage of the ML
solution life-cycle, the approach can also suggest additional training data or
the inclusion of multiple ML model variants as part of the overall ML
solution. The methodology may also assist in directing the labeling effort to
areas that would best cover the business requirements. That work utilized the
modeling process used by CT (Burroughs, Jain, and Erickson 1994; Cohen,
Snyder, and Rothermel 2006; Dalal et al. 1999; Grindal et al. 2006; Kuhn,
Kacker, and Lei 2013; Williams 2000; Wojciak and Tzoref-Brill 2014) to detect
weak business areas of the ML solution and strengthen them.

CT is a well-known testing technique to achieve a small and effective test
plan. The systematic combinatorial modeling methodology used by CT provides
insight on the business requirements and ties them to testing artifacts. We
utilized the combinatorial model as a means to detect gaps and weaknesses of
the ML solution. Business requirements often do not appear as features in the
ML solution training data. Rather, they can be defined as abstractions over
these features or may utilize additional information that is not always
available for learning. The CT model parameters may then include abstractions
as simple as dividing values into bins, such as defining ranges over age as
young, middle-aged and elderly, as well as more complex abstractions such as
defining mathematical functions over the relation between multiple features.

Using the mapping between CT model parameters and ML model features, we
defined data slices over the training data as follows. We observed that a
$t$-way interaction coverage requirement induces a division of the training
data into slices. Each slice corresponds to a specific value combination of
size $t$ and consists of all data records that map to that specific value
combination. Naturally, the number of data slices is equal to the number of
$t$-way value combinations in the CT model. Note that the slices are not
disjoint, since each data record maps to multiple value combinations of size
$t$, reflecting different business requirements that it captures. Data slices
over credit data may include records that belong to young women (a 1-way
interaction slice) and records that belong to young women with only high-
school level education (a 2-way interaction slice). We performed coverage gap
analysis of the training data associated with each data slice. Slices with no
corresponding training data expose business areas in which the ML solution
will likely under-perform. For slices with sufficient data or support, we
computed the ML metric over each slice. We expected the result for a slice to
be within the confidence interval of the ML metric over the entire data. If
that was not the case, we highlighted the slice as indicative of business
requirements that are not well addressed by the ML solution.

The work we report on here provides a set of heuristics, FreaAI, for applying
the above CT modeling methodology on the ML solution test data in order to
better validate it and verify that it addresses its requirements. As even the
simple abstractions were shown to be useful, FreaAI was designed to
automatically generate simple abstractions such as binning of continuous
features and identification of categorical feature values and combinations of
feature. The resulting data slices are guaranteed to: (1) under-perform
relative to the overall ML solution performance, (2) be explainable, (3) be
statistically significant. Details on the heuristics and how they guarantee
these properties follow in the Methodology Section.

## Methodology

FreaAI automatically analyzes features and feature interactions in order to
define feature models in which each of the data slices under-performs. It
requires only the test data and the ML prediction for each record, without
requiring knowledge of the ML model, for instance.

### Single feature analysis

FreaAI starts by analyzing single features. The rationale is that the simpler
the reason for an issue, the easier it is to understand its cause and, as a
result, easier to come up with a fix. FreaAI heuristically decides which
features should be treated as categorical and which as continuous and does so
automatically based on the type of feature data and number of unique values.
Alternatively, the user can supply this information. Looking for under-
performing data slices over single categorical features is simple. FreaAI
computes the ML performance of data slices that are defined by each of the
feature possible values (categories). Finding data slices over single
continuous features is challenging, as we need a way to find interesting
ranges. We do that in various ways, explained in the following paragraphs.

### Feature interactions analysis

FreaAI also analyzes feature interactions. It does so in two ways: (1) by
slicing on a category or range found when analyzing a single feature and
repeating the process for a single feature on that subset, or (2) by using a
method which can divide the search space on multiple axis at once—in our case,
a decision tree (DT) on n features. From CT we know that most software
problems involve a small number of variables. Empirically, we notice that this
seems to be the case for data features used for ML as well; single features
and interactions of two or three features already uncover many problems. An
inherent problem with increasing the interaction level is that the higher the
interaction level, the less likely it is to have sufficient support.

### Heuristics for feature and interaction analysis

We experimented with multiple feature-analysis techniques. These included
clustering methods, building decision trees, and Highest Posterior Density
(HPD), also known as High Density Regions (HDR) (Hyndman 1996) to extract
feature models with under-performing slices. We shortly describe each of the
techniques. We found decision trees and HPD to be the most effective,
therefore we report only on the results of applying these heuristics in the
Experimental Results section.

#### Clustering as a feature analysis heuristic

We applied clustering techniques to group the test data. We utilized multiple
clustering algorithms that are inherently different in concept, such as
k-means and DBScan. We then tested whether any of the resulting clusters
exhibited ML under-performance. We also varied the clustering metric and the
parameters to get different clusters. While this method worked well for
identifying under-performing slices, it was hard to automatically point out
the offending feature(s) in a way that was user-actionable. We could only mark
a (small) subset of features as being more dominant in the under-performing
cluster than in the rest of the test set, suggesting one or more of them could
be responsible.

#### HPD as a feature analysis heuristic

HPD is a nonparametric method that computes the shortest possible interval or
union of non-overlapping intervals containing a given proportion of records in
a univariate and numeric sample; it essentially gives the ‘best possible’
empirical confidence interval for a sample. An illustration from the original
paper is given in Figure 1.

Figure 1: Illustration of five different intervals containing 75% of the
density from an example (bimodal) empirical sample (source: (Hyndman 1996)).
The first one is the HDR (or HPD) since of the five, it is the shortest
(highest-density); in this case, the HPD is given by a union of two non-
overlapping intervals due to the bimodality of the distribution.

Since HPDs can only be used on univariate numeric data (unlike decision trees,
which do not have such restrictions), we use it in our heuristic only to
identify interesting slices for each numeric feature alone (i.e., no higher-
order interactions). So the slices are explainable and to avoid potential
overfitting due to multi-modality of the distribution, we furthermore restrict
the HPD to consist of a single interval, rather than a union. For the sample
of each relevant feature in the test dataset, we begin by finding the initial
HPD, and iteratively shrink that area by some epsilon (typically 0.05)
checking whether the performance on the new sample increases or decreases
relative to the previous one. If it decreases, then we have found a smaller
range which under-performs. Alternatively, if it increases, then at least one
of the discarded ranges (one on either side) under-performs. This process is
repeated after dropping the highest-density range from the feature until we
remain with too little data, e.g. 10% of the data records.

Figure 2 depicts one iteration of the HPD heuristic over one feature,
‘MFCCs_22’ from the Anuran (family) dataset. First, the area under the
probability density curve, up to and including the vertical ranges, is
identified by HPD as containing 90% of the probability density. Our HPD
heuristic iteration shrinks that area by 5% by removing the vertical ranges
and re-checks the performance of the updated range. It then proceeds as
described above.

Figure 2: Shows the areas discarded when changing HPD to include 85% [-0.04553
to 0.26878 with a performance of 0.9858] instead of 90% [-0.07738 to 0.27257
with 0.9840] density. The two marked bars are the areas FreaAI checks when the
performance of the lower density is higher. This data belongs to feature
‘MFCCs_22’ from Anuran (family).

#### Decision Trees as a feature analysis heuristic

The decision trees heuristic works as follows. For each subset of one, two or
three features the heuristic defines a decision tree classification problem
using these features as input. The output is always defined by a binary target
of True/False. ‘True’ means that the original ML model predicted the original
output successfully and ‘False’ means that it failed in its prediction. Notice
that this makes each classification problem binary, regardless of the number
of original labels. The decision tree heuristic then fits a decision tree
model to each of the above training data, but ignores the performance of the
resulting model. Instead, it walks the generated tree looking for splits where
the purity of ‘False’ is high enough for the slice defined by it to be under-
performing and keeps only those.

#### Minimal support

The slicing heuristics consider as candidates only slices that have at least
minimum support. Minimum support is heuristically calculated as the maximum
between 2 and 0.5% of mis-predicted instances in the test dataset. The
rationale is that in a dataset with a high metric performance there will be
few mis-predictions, therefore a slice with a small number of mis-predictions
is important, whereas for a model whose performance is relatively low there
will be many mis-predictions requiring a slice to have more mis-predictions to
be interesting. In addition, by taking a percentage we automatically adjust to
different test dataset sizes, e.g., Table 1 shows that the Statlog dataset has
300 records, while the Adult dataset has over 14,000 records. The percentage
is user-configurable. However, in our experiments the empirical values that we
chose were satisfactory.

### Defining data slices requirements

We define the following requirements from the data slices reported by FreaAI,
regardless of the heuristic that generated a data slice. A data slice must be

  1. 1.

correct,

  2. 2.

statistically significant, and

  3. 3.

explainable.

We use a methodology that abides to these requirements and generates slices
which are correct-by-construction. The slices are correct because they are
subsets of the test set for which we run the same performance metric as on the
original dataset. They are statistically significant because we take only
slices where the performance is considerably (4% by default) lower than the
expected lower-bound of the confidence interval of the performance metric.
FreaAI further reports only on data slices that have both minimal support and
a sufficiently low p-value, as computed by the hypergeometric distribution
test, described in the Experimental Results section. The slices are also self-
explainable since they are generated from an easily defined range on a small
number of features. Moreover, since we focus on structured datasets with
engineered features, the ML model developer should be able to draw useful
conclusions. Table 8 lists a few examples of data slices. These are
explainable by the way they are defined. For example, the first data slice
indicates a problem in the Adult ML model performance when the ‘relationship’
category is $10$. This is an example of a 1-way or single feature data slice.
Another data slice indicates a problem when the feature ‘hours_per_week’ has
values between $40$ and $43$ and the feature ‘workclass’ category is $5$. This
is an example of a 2-way or feature interactions data slice.

## Experimental results

Our experiments demonstrate the ability of our heuristics, as implemented by
FreaAI, to find under-performing data slices that are correct, statistically
significant, and explainable.

### Input

We tested FreaAI on open-source data sets, six taken from UCI and one from
ProPublica. The datasets are Adult, Avila, Anuran (Family, Genus, Species
targets), statlog (German credit data) from UCI and COMPAS Recidivism Risk
Score from ProPublica.

Table 1 lists some general information about the datasets and the performance
of the models trained on them. We made no attempt at creating the best
possible model for each dataset as it is not the goal of this work. We trained
a random forest model for each dataset and got ML performance that was similar
to the reported performance of each in the literature.

Table 1: Information about the datasets that we tested, listing the number of test records, the model performance measured by accuracy and the low/high bounds of the confidence interval for the accuracy Dataset | # records | Acc | Low CI | High CI
---|---|---|---|---
Adult | 14653 | 0.852 | 0.845 | 0.857
Avila | 6169 | 0.979 | 0.976 | 0.983
Anuran (family) | 2159 | 0.980 | 0.974 | 0.987
Anuran (genus) | 2159 | 0.967 | 0.959 | 0.974
Anuran (species) | 2159 | 0.969 | 0.962 | 0.977
Statlog | 300 | 0.763 | 0.708 | 0.821
ProPublica | 1852 | 0.742 | 0.718 | 0.759

In our experiments, FreaAI gets as input a test dataset. This is data that was
not used for training the ML model, is considered representative of the
underlying data distribution, and contains both the ground truth value of the
target feature as well as the ML model output, either a numeric or class
prediction as appropriate.

### Finding under-performing data slices

FreaAI implements various slicing heuristics to find candidate under-
performing data slices, as discussed in the Methodology Section. FreaAI
computes the performance by comparing the ground truth and the ML model
prediction that are provided as part of the input data, according to the same
ML performance metric used to assess the ML model performance. For simplicity,
in all the experiments that we report on here we use accuracy as the ML
performance metric and treat all features as numeric. As the Methodology
Section explains, we report on the data slices resulting from applying the HPD
and the DT slicing heuristics.

FreaAI’s heuristics consider as candidates only slices that have at least
minimum support, as the Methodology section describes. The candidate slices
are further filtered according to their statistical significance and FreaAI
reports only on statistically significant slices. The next subsection provides
details about the hypergeometric distribution test that FreaAI applies to
achieve that.

Tables 4–7 show that regardless of the heuristic used to calculate the slices
(HPD or DT) or the level of feature interaction (1-way or 2-way), the support
varies greatly for all datasets, with the occasional exception of the Statlog
dataset which is the smallest dataset (see dataset statistics in Table 1).
This is to be expected, as our goal is to find groups of mis-predicted records
for which we can automatically provide a readable explanation. Table 8
provides examples from the Adult dataset of under-performing data slices.

It is highly encouraging that our heuristics consistently find under-
performing data slices. Moreover, the reported slices often have substantial
support, increasing the usefulness of the report to the user.

Table 2: Summary of the number of under-performing slices for 1 and 2-way feature interaction and their total. ‘cand’ indicates the candidates found the HPD heuristic. ‘rep’ indicates how many candidate slices were actually reported out of the ‘cand’ list after applying the hypergeometric distribution test with p-value $<$ 0.05. | # 1-way | # 2-way | # total
---|---|---|---
Dataset | cand | rep | cand | rep | cand | rep
Adult | 26 | 25 | 82 | 82 | 108 | 107
Avila | 38 | 29 | 17 | 17 | 55 | 46
Anuran (family) | 80 | 58 | 386 | 371 | 466 | 429
Anuran (genus) | 105 | 61 | 723 | 638 | 828 | 699
Anuran (species) | 91 | 64 | 653 | 569 | 744 | 633
Statlog | 24 | 7 | 56 | 25 | 80 | 32
ProPublica | 12 | 4 | 13 | 7 | 25 | 11
Table 3: Summary of the number of under-performing slices for 1 and 2-way feature interaction and their total. ‘cand’ indicates the candidates found the DT (decision trees) heuristic. ‘rep’ indicates how many candidate slices were actually reported out of the ‘cand’ list after applying the hypergeometric distribution test with p-value $<$ 0.05. | # 1-way | # 2-way | # total
---|---|---|---
Dataset | cand | rep | cand | rep | cand | rep
Adult | 36 | 35 | 285 | 271 | 321 | 306
Avila | 120 | 84 | 928 | 689 | 1048 | 773
Anuran (family) | 148 | 113 | 1539 | 1265 | 1687 | 1378
Anuran (genus) | 209 | 142 | 2123 | 1498 | 2332 | 1640
Anuran (species) | 205 | 134 | 1881 | 1336 | 2086 | 1470
Statlog | 29 | 8 | 610 | 177 | 639 | 185
ProPublica | 1 | 1 | 24 | 22 | 25 | 23
Table 4: Summary of the support for 1-way HPD slices using a stricter p-value in the hypergeometric distribution test of 0.01 as threshold. The first column lists the number of reported slices. Then for each dataset the minimal (‘MIN’), average (‘AVG’), maximal (‘MAX’) and standard deviation (‘STD’) for the reported slices. Dataset | # | MIN | AVG | MAX | STD
---|---|---|---|---|---
Adult | 19 | 108 | 2601.5 | 9888 | 2931
Avila | 13 | 6 | 121.1 | 305 | 77.6
Anuran (family) | 32 | 2 | 109.3 | 197 | 51.9
Anuran (genus) | 36 | 3 | 161.1 | 351 | 70.2
Anuran (species) | 30 | 3 | 201.6 | 460 | 111.3
Statlog | 2 | 3 | 78.5 | 80 | 1.5
ProPublica | 3 | 24 | 485.7 | 942 | 341.9
Table 5: Summary of the support for 2-way HPD slices using a stricter p-value of 0.01 as a threshold. Dataset | # | MIN | AVG | MAX | STD
---|---|---|---|---|---
Adult | 81 | 108 | 1278.7 | 8749 | 1741.6
Avila | 10 | 6 | 28.2 | 137 | 37
Anuran (family) | 280 | 2 | 12.92 | 92 | 13.8
Anuran (genus) | 351 | 3 | 23 | 114 | 22.7
Anuran (species) | 391 | 3 | 31.1 | 231 | 37.5
Statlog | 11 | 3 | 11.82 | 29 | 7.1
ProPublica | 3 | 24 | 143 | 285 | 102.8
Table 6: Summary of the support for 1-way DT slices using a stricter p-value of 0.01 as a threshold. Dataset | # | MIN | AVG | MAX | STD
---|---|---|---|---|---
Adult | 22 | 108 | 2301.5 | 9941 | 2869
Avila | 44 | 6 | 55.8 | 236 | 51.4
Anuran (family) | 70 | 2 | 54.6 | 262 | 57.9
Anuran (genus) | 72 | 3 | 74 | 348 | 78.4
Anuran (species) | 71 | 3 | 73.4 | 451 | 100.1
Statlog | 4 | 3 | 80 | 124 | 30.1
ProPublica | 1 | 24 | 519 | 519 | 0
Table 7: Summary of the support for 2-way DT slices using a stricter p-value of 0.01 as a threshold. Dataset | # | MIN | AVG | MAX | STD
---|---|---|---|---|---
Adult | 202 | 108 | 2682.5 | 10736 | 3086.8
Avila | 368 | 6 | 75.4 | 366 | 66.9
Anuran (family) | 870 | 2 | 71 | 366 | 70.3
Anuran (genus) | 934 | 3 | 92.3 | 555 | 90.3
Anuran (species) | 818 | 3 | 104.1 | 582 | 111
Statlog | 86 | 3 | 87 | 176 | 41.2
ProPublica | 17 | 24 | 372.9 | 826 | 220.4
Table 8: Examples of under-performing slices taken from both the HPD and the DT heuristics on the Adult dataset. The average performance for Adult is roughly 0.85. The first 4 rows are examples of 1-way under-performing slices and the last 3 of 2-way slices. The first two rows represent categorical features, the next two continuous features, as decided by heuristic, and the following three come from the union of continuous feature ranges. Note that the last row represents a subset of the range in the third row, i.e. a slice on a secondary feature. ATTR | VALUE | SUP | PERF | p-val
---|---|---|---|---
relationship | 5 | 692 | 0.699 | 1.2E-25
education | 10 | 174 | 0.707 | 7.4E-07
age | 33-64 | 8538 | 0.797 | 1E-113
capital_gain | 3942-4934 | 122 | 0.688 | 3.8E-06
(hours_per_week, | 40-43, | 141 | 0.673 | 8.2E-08
workclass) | 5 |  |  |
(education_num, | (15, 16), | 122 | 0.680 | 1.4E-06
relationship) | (1, 2, 3, 4) |  |  |
(age, | 33-64, | 3019 | 0.749 | 1.3E-62
hours_per_week) | 41-99 |  |  |

### Slice significance testing

Our combinatorial procedure, along with the HPD and DT heuristics, return data
slices of combinations of 1–3 features, that is subsets of the test set
records, for which the model performance is worse than the performance on the
test set overall. The performance on these slices may not be significantly
worse, however. That is, in the case of a classifier, the proportion of
correctly-classified observations in these subsets may not be statistically
significantly lower than the proportion in the test set overall. Due to the
number of potential combinations, which grows polynomially with the number of
features, we wish to only return to the user slices with statistically-
significantly worse performance than the ML model overall performance.

To assess this, let the test dataset consist of $N$ records, of which
$K\in\\{0,\dots,N\\}$ are correctly classified. Consider a slice of
$n\in\\{1,\dots,N\\}$ records, of which $k\in\\{0,\dots,n\\}$ are correctly
classified. The model under-performance (accuracy $k/n$) on the slice is
considered significant if $k$ is low relative to $n\frac{K}{N}$, the expected
number of correctly-classified observation in a randomly drawn subset the same
size $n$ from the test set, which has overall accuracy $K/N$.

#### The hypergeometric distribution test

We use the discrete statistical hypergeometric distribution, which in a
general setup models the probability of a certain number of binary outcomes
($k$) from a random draw without replacement of a certain size ($n$) from a
finite population (of size $N$). In our work, we are concerned with how many
correct classifications (a binary outcome) are expected in a random
observation subset of a given size ($n$) from the test set (size $N$) with a
total number $K$ of correct classifications. This likelihood is the
probability mass function of the hypergeometric distribution with parameters
$N,\>n,\>K$ which are observed in the slice, that is, $\text{Pr}(X=x\mid
N,n,K)$. The lower-tailed p-value, which evaluates the significance of our
result, is simply $\sum_{x=0}^{k}\text{Pr}(X=x\mid N,n,K)$, below the observed
correct classification count $k$ in the slice.

Figure 3 illustrates the hypergeometric calculation for the statlog on the
slice of observations with feature value ‘credithistory’=5 (not shown in this
paper). The slice contains $n=21$ observations, of which $k=14$ (66.7%) are
correctly classified; the dataset overall has $N=300$ observations, of which
$K=230$ (76.7%) are correct. The hypergeometric distribution, which models the
probability of $x$ correct classifications out of $n$, in this case has
support on $x=\\{0,\dots,21\\}$. The dashed line (at 16.1, or 76.7% of 21)
indicates the number of correct classifications that would be expected for a
random draw of size $n=21$, that is, the $k$ giving same accuracy as on the
test set overall.

Since we are only dealing with cases where the performance is less than
expected ($k$ observed $<n\frac{K}{N}$), the p-value (shown in red in Figure
3) is the area to the left of the observed value $k=14$. The performance is
not below the expected value significantly enough, however, as the p-value
$\approx 0.193$. In order to achieve a significant p-value (below 0.05), we
would need $k\leq 12$ (accuracy at most 57.1%), indicated by the dotted line.
As discussed below, among slices with the same accuracy $k/n$, the
hypergeometric p-value gives more significance to larger slices ($n$ larger).
Thus, a slice where $k=28$ and $n=42$—both doubled from before, thus leaving
the slice accuracy, and the overall $K$ and $N$ unchanged—would have a p-value
of $\approx 0.04$ and would be significant. To have the same significance, in
terms of the p-value, as the actual slice, this larger slice would have about
29 correct classifications (69%), that is, more than twice the observed
$k=14$. This means that our original slice of size $n=21$ and accuracy 66.7%
is about as significantly as bad as a slice twice as large but with a slightly
better accuracy of 69%.

Figure 3: Hypergeometric distribution of correct classifications and p-value
($\approx 0.193$) for a slice of size $n=21$ based on the ‘credithistory’
feature of the statlog dataset.

#### Potential issues in computing significant slices

Our slicing heuristics require the slices to consist of all the records with
the set of given categorical feature values, or of all records with numeric
feature values in a given range, or to be a higher-order interaction of such
single-feature slices. By virtue of this restriction, not all possible random
subsets of size $n$ of records would qualify as a slice, if there were omitted
observations that would occupy the same range along one of the features. The
ML performance on an unconstrained random sample would be on average the same
as that on the test set overall. It is possible that, this constraint on the
relevant subsets inherently results in ML model performances on the slice that
are typically more different from the overall performance than the typical
random subset is. This may be investigated in further research. Another issue
we will investigate is how to integrate multiple hypothesis testing procedures
on the slices. This may need to account for dependence between the model’s
performance on difference slices due to overlap in the features included, to
ensure that slices we identify as significant are so in fact. However, we
believe that slices that are filtered out as insignificant using this
criterion are definitely not worth investigating further, so this
hypergeometric p-value can be considered an initial sanity check on the
results.

#### Selecting slices with more support

The procedure of filtering candidate slices based on their low p-values tends
to select the larger slices among the candidates as those with worse model
performance than the dataset overall. This is a good characteristic of our
algorithm, since larger slices tend to indicate more fundamental data issues.
Small slices are more likely to represent artefacts of overfitting or simply
less-interesting areas of model weakness.

The tendency to select larger slices is a direct result of the hypergeometric
distribution, specifically the fact it is based on sampling without
replacement from a finite sample (in contrast to the binomial distribution).
Given a fixed dataset, the variance of the proportion of correct
classifications (that is, of $X/n$ when $X$ is a hypergeometric-distributed
random variable) decreases as the sample size ($n$, here the slice size)
increases. Therefore, for two slices with the same performance (i.e.,
proportion of correct classifications), the larger will be a more significant
result; the performance can even be slightly better in the larger slice while
still being more significant than the smaller one. For instance, in Table 8 on
the Adult dataset, the slice based on ‘capital_gain’ has worse model accuracy
(0.688) than that based on the ‘age’ feature (0.797), but has a much smaller
slice size ($n$=122 vs 8,538). Nevertheless, despite the fact that accuracy on
the age slice is better, the larger size makes it a more significant instance
of model under-performance, with a much lower p-value, though in this instance
both slices have low p-values. A illustrative example of this property is
given in the discussion of Figure 3.

In addition to the variance aspect, which would hold in sampling with
replacement, the same level of under-performance is also more significant in a
larger slice when sampling without replacement because the slice size is
nearer to the finite sample of correctly-classified observations. Thus, the
statistical rationale of the p-value reinforces our logical sense that larger
diagnostic slices of data are more significant and helpful to the user.

## Discussion

We report on initial promising results of automatically suggesting under-
performing data slices.

In our experiments we applied FreaAI on ML test data. Data slices in general
and especially ones automatically generated such as by FreaAI are useful
throughout the ML solution development life-cycle. We are working on utilizing
data slices from the very early stages of ML solution development, starting
from the design of a solution that provides value. We believe that our
technology has other useful use cases. One use case may be to identify
potential features or feature combinations for de-biasing. Bias is a prejudice
that is considered to be unfair, such as refusing a loan based on the gender
or race of a person. Even though FreaAI works on existing data which may
already be biased and therefore the ML performance on it might be within
bounds, FreaAI may occasionally identify data slices that are indicative of
bias. Table 9 provides anecdotal evidence of that.

Table 9: Anecdotal relation between FreaAI analysis and features identified as requiring de-biasing. The first two rows are taken from statlog in which age$<$25 was identified as a biasing factor and the next 3 rows from ProPublica where number of priors ($>1$) biases the model. The last two rows show this captured by the DT method which aggregated several values, while the other rows where detected by the HPD method. ATTR | VALUE | SUP | PERF | p-value
---|---|---|---|---
age | 21 | 4 | 0.250 | 0.038
age | 22 | 7 | 0.428 | 0.045
Number_of_Priors | 4 | 119 | 0.588 | 7E-05
Number_of_Priors | 5 | 79 | 0.658 | 0.024
Number_of_Priors | 3 | 132 | 0.659 | 0.007
age | (19, 20, .., 29) | 124 | 0.670 | 1E-04
Number_of_Priors | (2, 3, 4, 5) | 519 | 0.660 | 9E-07

Another use case is to improve the overall ML solution as a result of
analyzing the data slices. One may add training data as characterized by data
slices or add alternative logic for inputs that belong to specific slices.
Computing data slices performance on the training data may assist in doing
that. Table 10 provides an example. Slices that have low performance even on
the training set and have low support may be indicative of a need for more
training data with the slice’s characteristics. The first row is an example of
that. Slices that have low performance on the training data but high support,
such as the second row, may be indicative of a need for alternative logic for
input belonging to such slices.

Table 10: Running FreaAI on the training dataset we sometimes encounter slices which under-perform (relative to the train CI) even in this dataset. The higher the model performance the fewer such slices are found. The first two rows in the table are from the Avila dataset and the last two from Anuran (Family). The LOW CI column is the lower bound of the training dataset confidence interval ATTR | VALUE | SUP | PERF | p-value | LOW CI
---|---|---|---|---|---
F4 | -0.633–0.634 | 12 | 0.960 | 6E-14 | 0.9992
F3 | -2.172–0.18 | 172 | 0.988 | 0 | 0.9992
MFCCs_4 | 0.783-1 | 45 | 0.977 | 0 | 0.9986
MFCCs_9 | -0.587–0.311 | 98 | 0.989 | 3E-05 | 0.9986

FreaAI can automatically find only issues for which data exists. We highly
recommend adding manually-defined validation features to address the challenge
of identifying missing requirements.

There is a challenge of combining and reporting on multiple data slices. We
are experimenting with ways to optimize the different factors, such as a
Pareto front, to take into account support vs performance, size vs range of a
slice and more. We suspect that the use case for which the slices are being
analyzed also affects the kind of slices that are more interesting. For
example, if our aim is to identify features that should be de-biased, or if it
is to direct the labeling effort for getting more training data, we may expect
the set of interesting slices suggested by the analysis in these cases to
differ from one another.

Using multiple slicing heuristics (e.g., HPD or DT) may result in overlapping
slices. This creates a challenge of whether and how to combine the overlapping
slices.

We experimented with extracting data slices of up to 3-way feature
interactions. It is algorithmically possible to extract slices of any n-way
interaction (up to the number of data features). However, we suspect this is
not useful because the smaller the dataset is, the more fragmentation the
interactions result in and the smaller the likelihood that data exists.

We experimented with multiple heuristics to suggest data slices, and more
methods exist. We are exploring ways to rank the resulting slices.
Architecturally, to simplify the design, we could use the same method for
identifying both 1-way and higher-order interaction slices. This could be done
with the DT heuristic, for example. However, it may be the case that different
slicing heuristics capture different aspects of unmet requirements. We plan to
experiment with prioritizing slices generated with various techniques and
develop metrics, including ones based on user feedback, to better estimate the
variety of unmet requirements captured.

## Conclusions

We address the challenge of finding under-performing data slices to validate
the ML performance of ML solutions. We developed automated slicing heuristics
and implemented them in FreaAI, such that the resulting slices are correct,
statistically significant and explainable. We experimented with seven open
datasets, demonstrating that FreaAI is able to consistently find under-
performing data slices. Moreover, the reported slices often have substantial
support, increasing the usefulness of the report to the user.

## References

  * Barash et al. (2019) Barash, G.; Farchi, E.; Jayaraman, I.; Raz, O.; Tzoref-Brill, R.; and Zalmanovici, M. 2019.  Bridging the gap between ML solutions and their business requirements using feature interactions.  In _Proceedings of the 2019 27th ACM Joint Meeting on European Software Engineering Conference and Symposium on the Foundations of Software Engineering_ , 1048–1058. ACM.
  * Breck et al. (2016) Breck, E.; Cai, S.; Nielsen, E.; Salib, M.; and Sculley, D. 2016.  What’s your ML test score? A rubric for ML production systems.  In _Reliable Machine Learning in the Wild - NIPS 2016 Workshop (2016)_.
  * Burroughs, Jain, and Erickson (1994) Burroughs, K.; Jain, A.; and Erickson, R. 1994.  Improved quality of protocol testing through techniques of experimental design.  In _SUPERCOMM/ICC_ , 745–752.
  * Cohen, Snyder, and Rothermel (2006) Cohen, M. B.; Snyder, J.; and Rothermel, G. 2006.  Testing across configurations: implications for combinatorial testing.  _SIGSOFT Softw. Eng. Notes_ 31(6): 1–9.
  * Dalal et al. (1999) Dalal, S. R.; Jain, A.; Karunanithi, N.; Leaton, J. M.; Lott, C. M.; Patton, G. C.; and Horowitz, B. M. 1999.  Model-Based Testing in Practice.  In _ICSE_ , 285–294.
  * Grindal et al. (2006) Grindal, M.; Lindström, B.; Offutt, J.; and Andler, S. F. 2006.  An evaluation of combination strategies for test case selection.  _Empirical Softw. Eng._ 11(4): 583–611.
  * Hyndman (1996) Hyndman, R. J. 1996.  Computing and Graphing Highest Density Regions.  _The American Statistician_ 50(2): 120–126.  doi:10.1080/00031305.1996.10474359.  URL https://www.tandfonline.com/doi/abs/10.1080/00031305.1996.10474359.
  * Kuhn, Kacker, and Lei (2013) Kuhn, D. R.; Kacker, R. N.; and Lei, Y. 2013.  _Introduction to Combinatorial Testing_.  Chapman & Hall/CRC.
  * Ng (2018) Ng, A. 2018.  AI Transformation Playbook How to lead your company into the AI era.  Landing AI.  URL “url–https://landing.ai/ai-transformation-playbook/?utm˙source=MLYList&utm˙medium=TextLink&utm˙campaign=Playbook˝.
  * Sculley et al. (2014) Sculley, D.; Holt, G.; Golovin, D.; Davydov, E.; Phillips, T.; Ebner, D.; Chaudhary, V.; and Young, M. 2014.  Machine Learning: The High Interest Credit Card of Technical Debt.  In _SE4ML: Software Engineering for Machine Learning (NIPS 2014 Workshop)_.
  * Williams (2000) Williams, A. W. 2000.  Determination of Test Configurations for Pair-Wise Interaction Coverage.  In _TestCom_ , 59–74.
  * Wojciak and Tzoref-Brill (2014) Wojciak, P.; and Tzoref-Brill, R. 2014.  System Level Combinatorial Testing in Practice – The Concurrent Maintenance Case Study.  In _ICST_ , 103–112.
  * Zinkevich (2018) Zinkevich, M. 2018.  Rules of Machine Learning: Best Practices for ML Engineering.  Google.  URL “url–https://developers.google.com/machine-learning/rules-of-ml/?utm˙source=google-ai&utm˙medium=card-image&utm˙campaign=training-hub&utm˙content=ml-rules˝.