Discovering Bands from Graphs

TL;DR

This paper introduces a method to identify band structures in graphs by ordering vertices and partitioning the adjacency matrix, using optimization and regression techniques, to reveal underlying temporal or structural patterns.

Contribution

It formulates a novel optimization framework for discovering band structures in graphs and proposes efficient algorithms for vertex ordering and band detection.

Findings

Polynomial-time band discovery with isotonic regression

Effective heuristic for vertex ordering

Empirical success in practical scenarios

Abstract

Discovering the underlying structure of a given graph is one of the fundamental goals in graph mining. Given a graph, we can often order vertices in a way that neighboring vertices have a higher probability of being connected to each other. This implies that the edges form a band around the diagonal in the adjacency matrix. Such structure may rise for example if the graph was created over time: each vertex had an active time interval during which the vertex was connected with other active vertices. The goal of this paper is to model this phenomenon. To this end, we formulate an optimization problem: given a graph and an integer $K$, we want to order graph vertices and partition the ordered adjacency matrix into $K$ bands such that bands closer to the diagonal are more dense. We measure the goodness of a segmentation using the log-likelihood of a log-linear model, a flexible family of distributions containing many standard distributions. We divide the problem into two subproblems: finding the order and finding the bands. We show that discovering bands can be done in polynomial time with isotonic regression, and we also introduce a heuristic iterative approach. For discovering the order we use Fiedler order accompanied with a simple combinatorial refinement. We demonstrate empirically that our heuristic works well in practice.

Tables

Table 1. Table 1: Basic characteristics of the datasets, the number of edges in lattice, see Section  4.4 .

Name	$\left|V\right|$	$\left|E\right|$	$E\left(L\right)$
DblpCF	176	457	794
DblpCP	154	348	717
Fb107	1034	26 749	61 279
Fb1912	747	30 025	62 842
Paleo	139	4428	8737
Mammals	2183	2 378 193	4 752 044
Synthetic	250 000	279 223	2 124 683

Table 2. Table 2: Basic statistics from experiments using exact and approximate borders. The columns contain running times, the number of borders, the number of flip tie-breakers, the number of random tie-breakers, the number of refinement rounds, initial (negative) score based on a fiedler order, and final score after the refinement.

	Approximate borders								Exact borders
Name	time	brd	iter	rnd	ref	initial	final		time	initial	final
DblpCF	0.2s	55	235	48	3	945	908		.02s	945	905
DblpCP	0.4s	53	701	135	3	966	927		.05s	966	918
Fb107	12m	476	7217	676	7	61 734	60 444		20s	61 723	60 427
Fb1912	5m	375	7357	813	4	43 212	42 909		3.2s	43 212	42 930
Paleo	4s	201	423	51	4	$-8645$	$-8906$		.13s	$-8645$	$-8906$
Mammals	33m	2975	2000	40		19 798			2m	19 798	19 798
Synthetic	37m	625	2000	40		6 956 048

Equations

$$s\left(X\right)=\left|X\cap E\right|\log a\left(X\right)+\left|X\setminus E\right|\log (1-a\left(X\right)),\text{where}a\left(X\right)=\frac{\left|X\cap E\right|}{\left|X\right|},$$

$$p(x\mid r)=\exp (q(x)+Z(r)+rS(x)),$$

$$\begin{split}\sum_{i=1}^{K}\sum_{c\in C_{i}}\log p(D(c)\mid r_{i})&=\sum_{i=1}^{K}\big{(}{\left|C_{i}\right|}Z(r)+\sum_{c\in C_{i}}q(D(c))+\sum_{c\in C_{i}}r_{i}S(D(c))\big{)}\\ &=\sum_{c\in A\setminus B}q(D(c))+\sum_{i=1}^{K}\big{(}{\left|C_{i}\right|}Z(r)+r_{i}\sum_{c\in C_{i}}S(D(c))\big{)}\quad.\end{split}$$

$$s\left(C\mid r\right)=\left|C\right|(Z(r)+ra\left(C\right)),\text{where}a\left(C\right)=\frac{1}{\left|C\right|}\sum _{c\in C}D(c).$$

$$s\left(\mathcal{U}\right)=\sum _{i=1}^{K}s\left(U_i\setminus U_{i-1}\right).$$

$$\sum _i^{K}s\left(C_{[b_i,e_i]}\right).$$

$$s\left(O_{k-1}(j)\right)+s\left(C_{[j+1,i]}\right)$$

$$\sum _{(x,y)}(f(x,y)-D(x,y){)}^2.$$

$$\frac{d}{dc}\sum _{(x,y)\in Y}(f(x,y)-c-D(x,y){)}^2=2\sum _{(x,y)\in Y}-f(x,y)+D(x,y)<0,$$

$$right\left(i;T\right)=\left\{t_1,\dots ,t_j\right\},\text{where}j=\arg \underset{j>i}{\max }a\left(t_{i+1},\dots ,t_j\right).$$

$$F=\left\{(x,y)\in B\mid (x+1,y)\notin B\text{ and }(x,y+1)\notin B\right\}$$

$$(V\setminus U)\setminus (V\setminus X)=(V\cap X)\setminus U,$$

$$h(X;r,t)=s\left(X\setminus U_{i-1}\mid r\right)+s\left(U_{i+1}\setminus X\mid t\right).$$

$$g(l,\delta ;r,t)=l(Z(r)-Z(t)+(r-t)\delta ).$$

$$\begin{array}{rl}h(X;r,t)& =(\left|X\right|-\left|U_{i-1}\right|)Z(r)+r(cs\left(X\right)-cs\left(U_{i-1}\right))\\ & +(\left|U_{i+1}\right|-\left|X\right|)Z(t)+t(cs\left(U_{i+1}\right)-cs\left(X\right))\\ & =\left|X\right|(Z(r)-Z(t))+(r-t)cs\left(X\right)+\text{const},\end{array}$$

$$\begin{split}&h(Y;r,t)-h(X;r,t)\\ &\quad={\left|Y\right|}d+(r-t)\mathit{cs}\mathopen{}\left(Y\right)-{\left|X\right|}d-(r-t)\mathit{cs}\mathopen{}\left(X\right)\\ &\quad={\left|Z\right|}\big{(}d+(r-t)\frac{\mathit{cs}\mathopen{}\left(Z\right)}{{\left|Z\right|}}\big{)}=g({\left|Y\setminus X\right|},\mathit{a}\mathopen{}\left(Y\setminus X\right);r,t)\quad.\end{split}$$

$$r_{j-1}>r_j\text{ for }j=2,\dots ,K\text{and}s\left({\mathcal{U}}^{\ast }\mid r_1,\dots ,r_K\right)\ge s\left(\mathcal{U}\right)-ϵ.$$

$$\begin{array}{rl}2ϵ& \ge h(Z)-h(U_i)=g(n,\beta )=cg(m,\beta )\\ & =cg(m,\alpha )+cm(r_i-r_{i+1})(\beta -\alpha )\ge cg(m,\alpha )\\ & =c(h(U_i)-h(X))\ge c(y-ϵ-h(X)),\end{array}$$

Full text

The goal of this paper is to model this phenomenon. To this end, we formulate an optimization problem: given a graph and an integer $K$, we want to order graph vertices and partition the ordered adjacency matrix into $K$ bands such that bands closer to the diagonal are more dense. We measure the goodness of a segmentation using the log-likelihood of a log-linear model, a flexible family of distributions containing many standard distributions. We divide the problem into two subproblems: finding the order and finding the bands. We show that discovering bands can be done in polynomial time with isotonic regression, and we also introduce a heuristic iterative approach. For discovering the order we use Fiedler order accompanied with a simple combinatorial refinement. We demonstrate empirically that our heuristic works well in practice.

††journal: Data Mining and Knowledge Discovery

Categories and Subject Descriptors

H.2.8 [Database management]: Database applications–Data mining

1 Introduction

Consider a dataset given in Figure 1(a). This data contains 139 species discovered at 501 sites (Fortelius, 2005). As different species live in different eras, the dataset can be sorted111Here, we sorted the data using the Fiedler order, see Section 5. such that the data points form a band. Let us construct a similarity matrix between the sorted species, where the weight between two species is the number of sites. Since a large number of the species-pairs will have do not share a single site, it is beneficial to view the matrix as a weighted graph. We see from the graph, given in Figure 1(b), that most of the edges will be located close to the diagonal, forming a band.

The phenomenon of having edges near the diagonal is not uncommon. For example, assume that the graph was constructed over time and that each vertex had an active time interval during which it connected to other active vertices with a higher probability. In such case, we should be able to arrange the vertices such that the edges are concentrated around the diagonal. As another example consider a graph with a(n overlapping) clustering structure. If we can rearrange clusters such that only neighboring clusters have significant overlap, if any, then there exists a vertex order such that the edges are close to the diagonal.

The goal of this paper is to quantify this phenomenon, see Figure 1(c) for an example. More formally, we introduce the following optimization problem. Given a graph and an integer $K$, order vertices and segment the adjacency matrix into $K$ bands such that a) each band respects the vertex order—when drawn the boundary of the band must either move down or right, b) the edge density in inner bands is higher than the density in the outer bands, and c) segments are as homogeneous as possible, according to some score function. Note that these bands may have varying thickness as illustrated in Figure 1(c).

In order to score the bands, we model the graph as a mixture model, where each band is a Erdős-Rényi model (Erdos and Renyi, 1960). Our goal is to find bands that optimize the likelihood of this model. As an application, discovered bands can be used for determining communities for individual vertices: if $(u,v)$ belongs to the $k$th band, then we say that $v$ belongs to the $k$th community of $u$, small $k$s corresponding to the inner circles.

We break this optimization problem into two natural subproblems. The first problem is to find optimal $K$ bands given the order and the second problem is to find a good order. We approach the latter problem by using Fiedler order (Fiedler, 1975) accompanied with a simple greedy refinement heuristic.

Most of this work is devoted into solving the first subproblem which happens to have a polynomial solution. In fact, this problem resembles a monotonic segmentation problem (see Haiminen and Gionis, 2004), however, it is much more intricate and there is no obvious technique for solving such problem. We will show that for certain scores, we do not have to consider all possible segmentations. We introduce a concept of borders. Roughly speaking, a border divides the adjacency matrix into two parts such that the inner part has a higher average. We will show that the optimal segmentation can be constructed from the borders. This allows us to transform the original optimization problem into two separate problems. First we need to discover all the borders, secondly we need to select the optimal segmentation using borders as candidates. Surprisingly, the second subproblem turns out to be an instance of the sequence segmentation problem, and can be solved using a standard dynamic program given by Bellman (1961).

We consider two techniques for discovering borders. The first approach is based on isotonic regression (Spouge et al., 2003), and gives us an exact solution. Despite being a polynomial-time solution, this approach requires that the graph is stored in a full form. Hence, we also present an iterative heuristic technique that uses sparsity of the graph to its advantage.

The rest of the paper is organized as follows. We introduce preliminaries and formally state our problem in Section 2. We introduce the concept of borders in Section 3 and present algorithms for discovering borders in Section 4. We consider discovering orders in Section 5. We present related work in Section 6 and experimental evaluation in Section 7. Finally, we present our conclusions in Section 8.

2 Preliminaries and Problem Statement

In this section we present our notation and give the formal problem statement. We first introduce the optimization problem for graphs and then cast this problem into a more general setup.

2.1 Discovering Bands from Graphs

Our first task is to define formally what we mean by a band. In order to do this, assume an undirected graph $H=(V,F)$. If we are given an order $o$ on vertices, essentially a mapping ${o}:{1,\ldots,{\left|V\right|}}\to{V}$, we say that $H$ respects the order if the neighborhood of each vertex can be seen as a segment w.r.t. the order, that is, for every $v\in V$, there exist integers $s$ and $e$ such that $\left\{u\in V\mid(v,u)\in F\right\}\cup\left\{v\right\}=\left\{o(i)\mid s\leq i\leq e\right\}$. If we order the vertices according to $o$, then $H$ will have all its edges next to diagonal. Our goal is given a graph $G$, find $o$ and $H$ optimizing a certain score.

Let us now define the score that we wish to optimize. Assume that we are given a graph $G=(V,E)$. Let $X\subseteq V\times V$ be a subset of vertex pairs. Let us define

[TABLE]

which is a maximum log-likelihood of a Bernoulli variable. We can now formulate the optimization problem.

Problem 1 (2-band discovery)

Given a graph $G=(V,E)$ find an order $o$ on vertices and a graph $H$ respecting that order $o$ and maximizing $s\mathopen{}\left(E(H)\right)+s\mathopen{}\left((V\times V)\setminus E(H)\right)$ such that $\mathit{a}\mathopen{}\left(E(H)\right)\geq\mathit{a}\mathopen{}\left((V\times V)\setminus E(H)\right)$.

In other words, we are modelling edges in $G$ as a mixture of two Bernoulli variables. The last constraint requires that the density of $G$ in the edges of $H$, that is, next to diagonal, should be higher than in the non-edges of $H$. As mentioned in the introduction, we are interested in a more general setup where we can discover several bands. This gives us the following optimization problem.

Problem 2 ($K$-band discovery)

Given a graph $G=(V,E)$ and integer $K$, find an order $o$ on vertices and $K+1$ graphs $H_{0},\ldots H_{K}$ respecting order such that $\emptyset=E(H_{0})\subsetneq\cdots\subsetneq E(H_{K})=V\times V$, the density is decreasing, $\mathit{a}\mathopen{}\left(C_{i}\right)\geq\mathit{a}\mathopen{}\left(C_{i+1}\right)$, and the score $\sum_{i=1}^{K}s\mathopen{}\left(C_{i}\right)$ is maximized, where $C_{i}=E(H_{i})\setminus E(H_{i-1})$.

In order to approach this optimization problem we will split it into two subproblems. The first subproblem is to find the bands for a fixed order.

Problem 3 (ordered $K$-band discovery)

Given a graph $G=(V,E)$ an integer $K$ and an order $o$ on vertices, find $K+1$ graphs $H_{0},\ldots H_{K}$ respecting order such that $\emptyset=E(H_{0})\subsetneq\cdots\subsetneq E(H_{K})=V\times V$, the density is decreasing, $\mathit{a}\mathopen{}\left(C_{i}\right)\geq\mathit{a}\mathopen{}\left(C_{i+1}\right)$, and the score $\sum_{i=1}^{K}s\mathopen{}\left(C_{i}\right)$ is maximized, where $C_{i}=E(H_{i})\setminus E(H_{i-1})$.

The second subproblem is to find the actual order. Our main focus will be the first subproblem for which we develop a polynomial exact solution. We address discovering the order in Section 5.

2.2 Band Discovery as Monotonic 2D-segmentation

Graph $H$ in ordered band discovery has a special property: if we order $H$ based on the order $o$ and consider the upper-half of the adjacency matrix, then we see that all 1s are concentrated next to the diagonal. We will use this observation to cast band discovery into a more general segmentation problem.

In order to do so, assume that we are given a dataset of size $M\times N$. Define $A=\left\{(a,b)\mid 1\leq a\leq M,1\leq b\leq N\right\}$ to be the set of all entries. We say that $U\subseteq A$ is a corner if for every $(a,b)\in U$, and for every $x$ and $y$ such that $1\leq x\leq a$ and $1\leq y\leq b$, an entry $(x,y)$ is a member of $U$.

Given an integer $K$ and a corner $B$, we define a $K$-segmentation to be the set of $K+1$ corners, $U_{0},\ldots,U_{K}$ such that $U_{0}=B$, $U_{K}=A$, and $U_{i-1}\subseteq U_{i}$ for each $i=1,\ldots,K$. We will refer to the difference set $U_{i}\setminus U_{i-1}$ as a segment.

Our goal is to segment given data into $K$ segments such that this segmentation maximizes the likelihood of a log-linear model for each segment. By log-linear models, also known as exponential family, we mean models whose probability density function can be written as

[TABLE]

where $S$ is a function mapping each data point to a real number, $r$ is the parameter of the model, and $Z(r)$ is the normalization constant. Many standard distributions such as Bernoulli, binomial, Gaussian, and Poisson are log-linear distributions. Interestingly, one can show with a direct computation that using a Gaussian distribution with a fixed variance corresponds to minimizing $L_{2}$ error.

Before we define our score, let us demonstrate that we can safely assume that $S(x)=x$ and $q(x)=0$. In order to do that, let $C_{1},\ldots,C_{K}$ be $K$ subsets of $A\setminus B$ such that $C_{1}\cup\cdots\cup C_{K}=A\setminus B$ and $C_{i}\cap C_{j}=\emptyset$ for $i\neq j$. Assume that we have assigned to each $C_{i}$, a parameter $r_{i}$.

We will measure the goodness of a segmentation by the log-likelihood of the model,

[TABLE]

Note that the first term does not depend on $C_{i}$ or $r_{i}$. Consequently, we can ignore it and by doing so ignore $q(x)$ term. We can also safely assume that $S(x)=x$. Otherwise, we can transform data $D$ to a new dataset $D^{\prime}$ by setting $D^{\prime}(c)=S(D(c))$.

We can now formally define our score. Given a segment $C$ and a parameter $r$, we define the score $s\mathopen{}\left(C\mid r\right)$ as

[TABLE]

We also define $s\mathopen{}\left(C\right)=\sup_{r}s\mathopen{}\left(C\mid r\right)$ to be the score of the optimal model. Given a $K$-segmentation $\mathcal{U}=\left(U_{0},\ldots,U_{K}\right)$, we define the score $s\mathopen{}\left(\mathcal{U}\right)$ as

[TABLE]

We say that a $K$-segmentation $\mathcal{U}$ is monotonic if $\mathit{a}\mathopen{}\left(U_{i+1}\right)\leq\mathit{a}\mathopen{}\left(U_{i}\right)$ for each $1\leq i\leq K-1$ such that $U_{i}\neq\emptyset$. Our goal is to solve the following problem.

Problem 4 (2D-segmentation)

Given a dataset $D$, a corner $B$, a log-linear model, and an integer $K$, find a monotonic $K$-segmentation $\mathcal{U}$ maximizing $s\mathopen{}\left(\mathcal{U}\right)$.

We can now see that band discovery is in fact an instance of 2D-segmentation. The dataset $D$ is in fact the adjacency matrix of $G$, the corner $B$ is a diagonal, and the score model is Bernoulli variable. We should point out that solving the more general problem has its advantages. If we have weights on edges, we can use some other log-linear model, such as Poisson model, to score the segmentation. Moreover, we do not have to restrict ourselves to graphs, we can segment any given matrix. On the other hand, discovering ordered bands is essentially as difficult as solving 2D-segmentation, that is, the amount of additional work we need to do to solve the more general case is negligible.

From now on, we will ignore $B$ for the sake of simplicity, and assume that we want to segment the whole dataset. The algorithms that we present can be easily adjusted to handle the more general case when $B$ is given.

The next two sections are devoted to solving Problem 4. We discuss discovering the vertex order in Section 5.

3 Borders

We can easily show that there are ${N+M\choose M}$ different corners for a data of size $M\times N$. This implies that we cannot enumerate all possible corners in order to solve Problem 4. However, we can show that we do not have to consider all possible corners.

In order to do so, we say that a corner $U$ is a border if there are no corners $X$ and $Y$ such that $X\subsetneq U\subsetneq Y$ and $\mathit{a}\mathopen{}\left(Y\setminus U\right)\geq\mathit{a}\mathopen{}\left(U\setminus X\right)$. We denote all borders by $\mathit{brd}\mathopen{}\left(D\right)$. An example of a non-border is given in Figure 2(a).

The next key proposition shows that we can safely ignore all corners that are not borders.

Proposition 1

Let $\mathcal{U}$ be a $K$-segmentation. Then there exists a $K$-segmentation $\mathcal{V}=\left(V_{0},\ldots,V_{K}\right)$ such that $V_{i}\in\mathit{brd}\mathopen{}\left(D\right)$ for every $i=0,\ldots,K$ and $s\mathopen{}\left(\mathcal{V}\right)\geq s\mathopen{}\left(\mathcal{U}\right)$.

We present the proof for Proposition 1 in Appendix A.

The following proposition states that the borders form a chain, a key property, which we also illustrate in Figure 2(b).

Proposition 2

Let $D$ be a dataset. Let $U,\ V\in\mathit{brd}\mathopen{}\left(D\right)$ be two borders. Then either $U\subseteq V$ or $V\subseteq U$.

Proof

Assume that the proposition does not hold, then both sets $U\setminus V$ and $U\setminus V$ are non-empty. Assume that $\mathit{a}\mathopen{}\left(U\setminus V\right)\geq\mathit{a}\mathopen{}\left(V\setminus U\right)$, otherwise swap $U$ and $V$. Let $X=U\cap V$ and $Y=U\cup V$. We have $\mathit{a}\mathopen{}\left(Y\setminus V\right)=\mathit{a}\mathopen{}\left(U\setminus V\right)\geq\mathit{a}\mathopen{}\left(V\setminus U\right)=\mathit{a}\mathopen{}\left(V\setminus X\right)$. Consequently, $V$ is not a border, which is a contradiction. ∎

Proposition 2 has several key implications. Assume that we have a dataset $D$ of size $N\times M$. Then the number of borders is bounded by $NM+1$. From now on we will assume that $\mathit{brd}\mathopen{}\left(D\right)=U_{0},\ldots,U_{L}$ is ordered from smallest to largest. In order to store this list efficiently, we define $C_{i}=U_{i}\setminus U_{i-1}$ for $i=1,\ldots,L$. Instead of storing $U_{i}$ we simply store $C_{i}$. Since $C_{i}$ are disjoint and the union of $C_{i}$ is equal to $A$, the set of all entries, we can store the borders in $O(NM)$ space.

Propositions 1–2 allow us to divide Problem 4 into two subproblems. The first problem is to find all the borders.

Problem 5

Given a dataset $D$, compute $\mathit{brd}\mathopen{}\left(D\right)$.

After we have discovered the borders, we can now use them to find the optimal segmentation. We will restate the problem in a slightly different manner, using directly segments instead of corners. In order to do so, given a list of segments $C_{1},\ldots,C_{L}$, let us define $C_{[a,b]}$ to mean $\bigcup_{i=a}^{b}C_{i}$. Note that if we set $C_{i}=U_{i}\setminus U_{i-1}$, then it follows that $C_{[a,b]}=U_{b}\setminus U_{a-1}$. This implies that we can reformulate the optimization problem as follows.

Problem 6

Given a sequence of $L$ segments $C_{1},\ldots,C_{L}$ and an integer $K$, find $K$ intervals, $[b_{i},e_{i}]$, such that $b_{1}=1$, $e_{K}=L$, and $b_{i}=e_{i-1}+1$, for $i=2,\ldots,K$, optimizing

[TABLE]

Note that we have dropped the monotonicity condition from the definition of the problem. We will see later in Corollary 1 that if $C_{i}$ are constructed from borders, that is, $C_{i}=U_{i}\setminus U_{i-1}$, then $\mathit{a}\mathopen{}\left(C_{i+1}\right)<\mathit{a}\mathopen{}\left(C_{i}\right)$. Thus monotonicity will automatically be guaranteed.

Problem 6 is in fact an instance of a classic sequence segmentation problem, where the goal is to split a sequence of length $L$ into $K$ homogeneous segments. This can be solved by a dynamic program in $O(L^{2}K)$ time and $O(KL)$ space (Bellman, 1961). To see this let us write $O_{k}(i)$ to be an optimal $k$-segmentation for a sequence $C_{1},\ldots,C_{i}$. Then $O_{k}(i)$ is equal to $O_{k-1}(j)$ augmented with $[j+1,i]$, and $j<i$ is selected such that the score

[TABLE]

is maximized. We can iteratively compute this by first computing $O_{1}(i)$ for each $i=1,\ldots,L$, and use the above identity to $O_{k}(i)$ from $O_{k-1}(i)$ until we reach $K$ segments.

To summarize, we discover bands in 3 steps: {enumerate*}

Order the dataset, if one is not given.

Compute borders $\mathit{brd}\mathopen{}\left(D\right)$ (Problem 5).

Segment borders to obtain $K$-segmentation (Problem 6).

We have already shown that Problem 6 can be solved by using the classic segmentation technique. In the next two sections we discuss how to obtain the borders, either exactly or approximately. Finally, in Section 5 we discuss how to obtain the order. As discovering the order seems to be computationally intractable we resort to spectral heuristics.

4 Discovering borders

In this section we focus on discovering borders (Problem 5). Namely, we provide a polynomial-time dynamic program that discovers borders correctly. In addition we provide a heuristic for the cases when the exact discovery is too time-consuming.

4.1 Maximal and Minimal Corners

In order to define the discovery algorithm, we need to introduce the notion of maximal and minimal corners. We will show that these notions are closely related to borders.

Let $U$ be a corner. We define $\mathit{left}\mathopen{}\left(U\right)$ to be the minimal corner $V\subsetneq U$ such that $\mathit{a}\mathopen{}\left(U\setminus V\right)$ is the smallest possible. If there are several candidates, we choose the smallest one.​222We can easily show that this choice is unique. We also define $\mathit{right}\mathopen{}\left(U\right)$ to be the maximal corner $V\supsetneq U$ such that $\mathit{a}\mathopen{}\left(V\setminus U\right)$ is the largest possible. If there are several candidates, we choose the largest one.

We can use maximal and minimal corners to describe borders.

Proposition 3

Let $U$ and $V$ be two consecutive borders. Then $U=\mathit{left}\mathopen{}\left(V\right)$ and $V=\mathit{right}\mathopen{}\left(U\right)$.

We prove this proposition in Appendix B.

This proposition has an important corollary that shows why we can ignore the monotonicity condition in Problem 6. Indeed segments between borders will have automatically decreasing average.

Corollary 1

Let $U,V,W\in\mathit{brd}\mathopen{}\left(D\right)$ be three consecutive borders, $U\subsetneq V\subsetneq W$. Then $\mathit{a}\mathopen{}\left(V\setminus U\right)>\mathit{a}\mathopen{}\left(W\setminus V\right)$.

Proof

Proposition 3 implies that $\mathit{right}\mathopen{}\left(U\right)=V$. Since $V\subsetneq W$, Lemma 1 (given in Appendix B) implies that $\mathit{a}\mathopen{}\left(V\setminus U\right)>\mathit{a}\mathopen{}\left(W\setminus V\right)$. ∎

4.2 Exact discovery

In this section we present an algorithm for computing the borders.

Discovering borders exactly is in fact an instance of isotonic regression for a grid. In this regression we are given a grid, and a set of values on each grid entry. The goal is to find another set of values such that they decrease as we move towards the one selected corner and $L_{2}$ error is minimized as we move towards that corner.

Problem 7

Let $D$ be a dataset of size $M\times N$. Find a function $f$ such that $f(x,y)\geq f(x,y+1)$ and $f(x,y)\geq f(x+1,y)$ minimizing the cost

[TABLE]

Once these these values are discovered, we can reconstruct the borders using the following proposition.

Proposition 4

Let $D$ be a dataset and let $f$ be the solution to the grid isotonic regression. Let $B$ be a border. Then there is $\sigma$ such that $B=\left\{(x,y)\mid f(x,y)\geq\sigma\right\}$.

Proof

Let $m=\min_{(x,y)\in B}f(x,y)$. Define $Y=\left\{(x,y)\notin B\mid f(x,y)\geq m\right\}$. We need to show that $Y=\emptyset$. Assume otherwise. Note that since $f$ is monotonic, $B\cup Y$ is a corner. Since $B$ is a border, Proposition 3 implies that $\mathit{a}\mathopen{}\left(Y\right)<m$. An entry, if such exists outside $B\cup Y$ must be lower than $m$. Hence, we can decrease the values of $f$ in $Y$ by a small amount, say $\epsilon$, without violating the monotonic constraint. Let us consider the effect. In order to do this, consider the following derivative,

[TABLE]

where the inequality follows from the fact that $\mathit{a}\mathopen{}\left(Y\right)<m\leq f(x,y)$ for any $(x,y)\in Y$. Hence, we can decrease the values of $f$ in $Y$ by a small amount such that the monotonicity still holds and the score is decreased. This contradicts the fact that $f$ is the optimal solution. ∎

This proposition gives us a simple approach to discover borders by varying $\sigma$.

Finding the solution for grid isotonic regression can be done in $O((NM)^{2})$ time and $O(NM)$ space by an algorithm of Spouge et al. (2003). However, the time complexity of the algorithm is overly pessimistic. The algorithm runs in $O((NM)L)$, where $L$ is the number of borders. The number of borders is $NM$, at worst, but in practice it is much smaller. The algorithm of Spouge et al. (2003) is a conquer-and-divide algorithm. The running time $O((NM)L)$ is based on the pessimistic assumption that each division is imbalanced, that is, only one point is separated. If these divisions are (nearly) balanced, then in practice the running time will be closer to $O((NM)\log L)$. Additional speed-ups are possible if instead of considering the full matrix of size $NM$ we first compute the smallest corner containing the whole data, and apply the algorithm to the corner. The points that are left outside constitute a border with [math] average. This trick may speed-up the search significantly but it is highly vulnerable to noise as one non-zero point is enough to counter this optimization.

4.3 Heuristic discovery

The computational complexity of isotonic regression may become too high in practice, especially due to the $O(NM)$ term. Hence, in this section we present a practical heuristic approach. The idea here is to sort individual entries of $D$ into a sequence. Once we have this sequence we can use it to compute candidates for borders. We can then use these candidates to rearrange the entries again, and repeat the procedure until convergence.

To make this more formal, assume that we are given a dataset $D$ of size $M\times N$. Let $T=\left(t_{1}=(x_{1},y_{1}),\ldots,t_{NM}=(x_{NM},y_{NM})\right)$ be a sequence of all entries of $D$. We say that $T$ is a monotonic entry order if for any $k$, $1\leq k\leq NM$, the set $t_{1},\ldots,t_{k}$ is a corner.

Given a monotonic entry order $T=t_{1},\ldots,t_{NM}$ and an integer $i$, we define

[TABLE]

If there are ties, we select the largest index.

Given a monotonic entry order $T$, consider the following process. Set $U_{0}=\emptyset$ and then iteratively compute $U_{i}=\mathit{right}\mathopen{}\left({\left|U_{i-1}\right|};T\right)$ until we reach $U_{L}$ containing all the entries. We will write $\mathit{brd}\mathopen{}\left(T\right)=U_{0},\ldots,U_{L}$.

We say that a monotonic entry order $T=t_{1},\ldots,t_{NM}$ is compatible with $\mathit{brd}\mathopen{}\left(D\right)$ if for each border $U\in\mathit{brd}\mathopen{}\left(D\right)$ there exists an index $k$ such that $t_{1},\ldots,t_{k}=U$. Such order always exists.

Assume that we are given an order $T$ compatible with $\mathit{brd}\mathopen{}\left(D\right)$. Select and fix a border $U_{i}\in\mathit{brd}\mathopen{}\left(D\right)$. Let $k={\left|U_{i}\right|}$. We must have $U_{i+1}=\mathit{right}\mathopen{}\left(k;T\right)$. This implies that if we are given an order $T$ that is compatible with $\mathit{brd}\mathopen{}\left(D\right)$, then $\mathit{brd}\mathopen{}\left(T\right)=\mathit{brd}\mathopen{}\left(D\right)$.

A naive approach to compute an individual $\mathit{right}\mathopen{}\left(i;T\right)$ requires $O(NM)$ time, however, we can compute $\mathit{right}\mathopen{}\left(i;T\right)$ for each $i$ simultaneously in $O(NM)$ total time using the approach given by Calders et al. (2007), where the goal of the authors was to discover the densest interval, essentially $\mathit{right}\mathopen{}\left(i;T\right)$, from a stream $T$ in an amortized constant time for each $i$. This algorithm is actually a variation of PAVA algorithm for solving isotonic regression for a total order, see (Ayer et al., 1955), though the goal and the output of the algorithms are different. Connection to isotonic regression is natural as we saw previously that the exact borders can be discovered by solving the grid isotonic regression in $O(N^{2}M^{2})$ time. Since we no longer consider a grid but a monotonic entry order, we can use a more efficient algorithm whose copmputational complexity is $O(NM)$.

Let us now consider the optimization problem from an another angle. Assume that we do not know $\mathit{brd}\mathopen{}\left(D\right)$ but instead we know only the average of values in segments, that is, for each entry $(i,j)$, we know the average, say $w(i,j)$, of the segment that contains $(i,j)$. We can construct the order compatible with $\mathit{brd}\mathopen{}\left(D\right)$ from the weights. Corollary 1 implies that weights should decrease. Hence, if we build a monotonic entry order by greedily selecting the entry with the largest weight while at the same time making sure that the order is indeed monotonic, we end up with an order that is compatible with $\mathit{brd}\mathopen{}\left(D\right)$. We present the pseudo-code for this approach in Algorithm 1.

We can now construct our algorithm for discovering approximate borders. Given a monotonic entry order $T=t_{1},\ldots,t_{NM}$, we can compute the borders $\mathit{brd}\mathopen{}\left(T\right)=U_{0},\ldots,U_{L}$. Assume that we are given an entry $p$. Let $U_{k}\in\mathit{brd}\mathopen{}\left(T\right)$ be the border such that $p\in U_{k}\setminus U_{k-1}$. Let us define $w_{S}(p;T)=\mathit{a}\mathopen{}\left(U_{k}\setminus U_{k-1}\right)$ to be the density of the segment containing $p$. Once these weights are computed, we can use them to compute a new order using FindOrder. Computing weights can be done in $O(NM)$ time while computing a new order can be done in $O(NM\log\min(N,M))$ time.

Using just $w_{S}$ is problematic in practice. The reason for this is that during FindOrder there will be often several entries in $H$ that belong to the same segment, and so will have the same values of $w_{S}$. Hence, we need a weight function that that will break these ties. Breaking ties properly is important since it is possible to design a weight $w$ that $T=\textsc{FindOrder}(w(\cdot;T))$ for any $T$. In other words, iterating FindOrder and computing weights will never improve the current order.

Given an order $T$, we say that a weight function $w$ is a tie-breaker if $w_{S}(p;T)<w_{S}(q;T)$ implies $w(p;T)<w(q;T)$, where $p$ and $q$ are entries in $T$. Before considering any specific tie-breakers, let us first show that using a tie-breaker $w$ leads to a convergence.

Proposition 5

Let $w$ be a tie-breaker. Let $T^{0}$ be any monotonic entry order and define $T^{i+1}=\textsc{FindOrder}\mathopen{}\left(w(\cdot;T^{i})\right)$. Then there exists $k$ such that $\mathit{brd}\mathopen{}\left(T^{i}\right)=\mathit{brd}\mathopen{}\left(T^{k}\right)$ for any $i\geq k$.

Proof

Fix $m$ and let $U=T^{m}$ and $V=v_{1},\ldots,v_{NM}=\textsc{FindOrder}\mathopen{}\left(w(\cdot;U)\right)$. Let $B_{0},\ldots,B_{L}=\mathit{brd}\mathopen{}\left(U\right)$ and define $C_{0},\ldots,C_{K}=\mathit{brd}\mathopen{}\left(V\right)$. Define a vector $\alpha$ of length $2L$ such that $\alpha_{2i-1}=\mathit{a}\mathopen{}\left(B_{i}\setminus B_{i-1}\right)$ and $\alpha_{2i}={\left|B_{i}\right|}$. Define similarly $\beta$ using $\mathit{brd}\mathopen{}\left(V\right)$.

Assume two entries $p$ and $q$ such that $p\in B_{i}$ and $q\notin B_{i}$. This means that $w_{S}(p;U)>w_{S}(q;U)$. Since $w$ is a tie-breaker, this means that $p$ will occur earlier than $q$ in $V$. In other words, $\{v_{1},\ldots,v_{{\left|B_{j}\right|}}\}=B_{j}$ as sets, for any $j=1,\ldots,L$.

If $\mathit{brd}\mathopen{}\left(U\right)=\mathit{brd}\mathopen{}\left(V\right)$, then it follows immediately that $\alpha=\beta$.

Assume that $\mathit{brd}\mathopen{}\left(U\right)\neq\mathit{brd}\mathopen{}\left(V\right)$ and let $j$ be the largest index such that $B_{i}=C_{i}$ for any $i\leq j$. This implies that $\alpha_{i}=\beta_{i}$, for $i\leq 2j$. Since $\{v_{{\left|B_{j}\right|}+1},\ldots,v_{{\left|B_{j+1}\right|}}\}=B_{j+1}\setminus B_{j}$ and $C_{j}=B_{j}$, we know that $C_{j+1}\setminus C_{j}$ will be as dense as $B_{j+1}\setminus B_{j}$, that is, it follows that either $\mathit{a}\mathopen{}\left(C_{j+1}\setminus C_{j}\right)>\mathit{a}\mathopen{}\left(B_{j+1}\setminus B_{j}\right)$, or $\mathit{a}\mathopen{}\left(C_{j+1}\setminus C_{j}\right)=\mathit{a}\mathopen{}\left(B_{j+1}\setminus B_{j}\right)$ and ${\left|C_{j+1}\right|}>{\left|B_{j+1}\right|}$. That is, either $\alpha_{2j+1}<\beta_{2j+1}$, or $\alpha_{2j+1}=\beta_{2j+1}$ and $\alpha_{2j+2}<\beta_{2j+2}$. Consequently, $\beta$ is larger than $\alpha$ with respect to the lexicographical order.

We have shown that if $\mathit{brd}\mathopen{}\left(T^{i}\right)\neq\mathit{brd}\mathopen{}\left(T^{i+1}\right)$, then for any $j>i$, $\mathit{brd}\mathopen{}\left(T^{i}\right)\neq\mathit{brd}\mathopen{}\left(T^{j}\right)$. Since there are only finite number of possible borders, there exists an index $k$ such that $\mathit{brd}\mathopen{}\left(T^{i}\right)=\mathit{brd}\mathopen{}\left(T^{k}\right)$ for any $i\geq k$. ∎

We consider two tie-breakers. The first tie-breaker, $w_{R}$ breaks the ties in a random order. Formally, we define $w_{R}(p;T)=(w_{S}(p;T),r)$ to be a vector of length 2, where $T$ is an order, $p$ is an entry in $T$, and $r$ is a random number between [math] and $1$. The weights are compared in lexicographical order. This immediately implies that $w_{R}$ is a tie-breaker. Our second tie-breaker tries to flip the original order as much as possible. Formally, we define $w_{F}(p;T)=(w_{S}(p;T),i)$, where $T$ is an order, $p$ is an entry in $T$, and $i$ is the index of $p$ in $T$. Obviously, $w_{F}$ is a tie-breaker, and it will favor the entries that appear later in $T$.

Proposition 5 states that the iteration will converge. However, we need means to detect this convergence. This is difficult since while the borders themselves will converge, the actual orders may cycle. Fortunately, we can show that $w_{F}$ has a cycle of at most 2, we prove this proposition in Appendix C.

Proposition 6

Let $T^{0}$ be any monotonic entry order and define iteratively $T^{i+1}=\textsc{FindOrder}\mathopen{}\left(w_{F}(\cdot;T^{i})\right)$. There exists $m$ such that $T^{m+2}=T^{m}$.

We can easily detect convergence for $w_{F}$ by remembering the second last order and compare it to the current one. If the orders are the same, then we have converged.

Unfortunately, we cannot detect convergence easily for $w_{R}$. Hence we adopt the following hybrid approach. We begin with an random monotonic order, and apply $w_{F}$ until convergence. Once converged, we apply $w_{R}$ once, and repeat applying $w_{F}$ until we converge again. We repeat this until the borders have not changed after 20 applications of $w_{R}$. This is summarized in Algorithm 2.

4.4 Using sparsity to speed-up the discovery

So far we have operated with ordinary matrices, and consequently the running time for our heuristic discovery is $O(N^{2}\log N)$ when applied to a graph with $N$ vertices. In practice, graphs are sparse and we can use this sparsity to our advantage.

Assume that we are given a graph $G$ with positive weights and a corner $B$. Let us define

[TABLE]

to be the set of frontier points, that is, points in $B$ from which you cannot advance away from the diagonal, see Figure 3(a). Now, Proposition 3 implies that if $B$ is a border, then $F$ must be a subset of edges. Otherwise, we can always remove a non-edge frontier point and by doing so increase the density of $B$. Note also that given $F$, we can always recover $B$ by taking the smallest corner containing $F$.

This suggests that instead of ordering all entries of the adjacency matrix it is enough to order just the edges. This, however, poses two complications. First, we need to be able keep the edges in a monotonic order. While this was easy when dealing with the full matrix since visiting one entry revealed only 2 adjacent entries, at maximum, visiting an edge may reveal any number of unvisited edges. Secondly, in order to compute the density of $B$, whenever we visit an edge we need to compute the number of non-edges we visited implicitly in order to reach that edge, see Figure 3(c).

In order to solve the first problem, we define a lattice structure on the edges. More formally, we define a directed acyclic graph $L$ such that the vertices correspond to the edges of $G$ and the edges in $E(L)$ are formed as follows: vertex $(x_{1},y_{1})$ is connected to $(x_{2},y_{2})$ in $L$ if $x_{1}\leq x_{2}$ and $y_{1}\leq y_{2}$ and there is no third vertex $(x_{3},y_{3})$ such that $x_{1}\leq x_{3}\leq x_{2}$ and $y_{1}\leq y_{3}\leq y_{2}$, see Figure 3(b), for example. In order to guarantee that the edges are visited in a monotonic order, we simply visit the edges in a topological order of $L$. This modification of FindOrder can be done to run in $O({\left|E(L)\right|}+{\left|E(G)\right|}\log{\left|V(G)\right|})$ time. Moreover, the lattice can be constructed in $O({\left|E(L)\right|}\log{\left|V(G)\right|})$ time.

In order to solve the second problem, that is, to compute the number of non-edges encapsulated by visiting a new entry, see Figure 3(c), first note that we need this information when computing the borders from the monotonic order. Consequently, at this point we have an order at our disposal. We enumerate the entries in the order and during this enumeration we compute and update the frontier set of the border corresponding to the entries visited so far. This update can be done efficiently by keeping the entries, which are edges in $G$ and have a form, say $(x,y)$, in a red-black tree indexed by $x$. Adding a new entry into a border will delete at most $k$ entries from the frontier set, where $k$ is the number of parents in the lattice. There can be at most ${\left|V(G)\right|}$ frontier entries, at any time. Hence, the running time for this enumeration is $O({\left|E(L)\right|}\log{\left|V(G)\right|})$. Assume now that we have computed the frontier set for $i-1$ entries, and we need to compute the encapsulated non-edges for $i$th entry, say $(x,y)$. We do this by finding the entry from the frontier set $(u,v)$ such that $u$ is the largest possible value and $u\leq x$. We then use this entry as a starting point and iterate to the following entries w.r.t. the red-black tree. At each entry we compute the number of non-edges captured between two adjacent nodes of the frontier set and $(x,y)$, see Figure 3(c) for illustration. We need to visit only $O(k)$ entries, where $k$ is the number of parents of $i$th entry in $L$. Hence, computing the areas for all entries can be done in $O({\left|E(L)\right|}\log{\left|V(G)\right|})$ time.

5 Discovering Order

Our main focus so far was to compute the segmentation from already ordered data. This order may be given, for example, if each vertex has a time stamp. If the order is not known, we will compute the order using the Fiedler vector (Fiedler, 1975), a popular technique for ordering matrices and graphs. Fiedler order has a tendency of pushing the edges towards the diagonal. In fact, we can show that if there exists an order of vertices of the graph $G$ such that the edges respect the order, that is, all edges are around diagonal, then Fiedler vector will find this order (Fiedler, 1975).

We use Fiedler order as an initial order and introduce a simple heuristic refinement. Assume that we have computed a segmentation using the order. Let $B$ one of the segments and let $(x,y)$ be a frontier, as defined in Section 4.4, of $B$. Find the smallest vertex $x^{\prime}$ such that for any $u$, $x^{\prime}\leq u\leq x$, it holds that for any $v$, $(x,v)$ and $(u,v)$ belong to the same segment. Similarly, find the largest vertex $y^{\prime}$ w.r.t. $y$, see Figure 4(a). It immediately follows that a permutation of vertices between $x^{\prime}$ and $x$, and $y$ and $y^{\prime}$ cannot decrease the score since none of the entries will leave its segment. Assume that there is a non-edge entry $(u,v)$ such that $x^{\prime}\leq u\leq x$ and $y\leq v\leq y^{\prime}$. Then if we swap $x$ and $u$, and $v$ and $y$, we essentially decrease the area of the segment since none of the entries will leave the segment but at the same time $(u,v)$ cannot be a frontier since it is non-edge, see Figure 4(b). If there are several non-edges, we select the non-edge, say $(u,v)$, minimizing ${\left|(u,z)\in E(G)\mid y\leq z\leq y^{\prime}\right|}+{\left|(z,v)\in E(G)\mid x^{\prime}\leq z\leq x\right|}$, that is, we select $u$ and $v$ with the smallest number of edges.

We will do these swaps until no swaps are possible. For the sake of efficiency, we do these swaps in a batch style, that is, for each frontier point $(x,y)$ we first compute $x^{\prime}$ and $y^{\prime}$, and then perform swap. During these swaps we make sure that once an interval $[x^{\prime},x]$ or $[y,y^{\prime}]$ is used for a swap, it will not be used again in the same batch. Once no swaps are possible, we recompute the segmentation, and repeat the refinement. This will eventually converge since the score will always increase after the refinement.

6 Related Work

To our knowledge our proposed optimization problem is novel. However, there are similar tasks. A common plotting technique is a contour plot, where a line follows a constant value. Our optimization technique is inherently different as we are considering all the points inside a segment and not just the data points within the vicinity of the boundary.

A related optimization problem was introduced by Mannila and Terzi (2007), where the authors introduced a concept of nestedness for binary data. A binary data is nested if for all pairs of rows, one row is either a superset or a subset of another. If data is fully nested, we can order it in a way that all 1s occur in the top-left corner. The authors considered minimizing number of 0s that needs to be replaced such data becomes fully nested. Authors also consider partitioning columns into $K$ parts such that each part is almost fully segmented. While both approaches model the same phenomenon, there are significant differences. We can handle more general datasets such as counting data and real-valued data, and our scoring function is based on the log-likelihood. We are interested in finding $K$-segmentations, whereas Mannila and Terzi (2007) focus on 2-segmentation. On the other hand, we assume that our data is preordered whereas no such assumption is made by Mannila and Terzi (2007).

The discovered bands give rise for each vertex a set of expanding communities. Discovering similar structures have been suggested. For example, Alvarez-Hamelin et al. (2006) order nodes into a tree by deleting iteratively low-degree vertices. In another example, Tatti and Gionis (2013) suggested discovering nested communities given a set of seed nodes.

The spectral approach to discover an order falls into a larger category of approaches used for solving the seriation problem. In seriation, a typical goal is, given a similarity matrix between objects, to organize the objects that minimize some stress function (Hahsler et al., 2008). In the ideal case, the discovered order will rearrange the similarity matrix into a Robinson form, that is, the values of the matrix will get smaller as we move away from the diagonal. In practice, such a permutation rarely exists and hence the stress function reflects how far away we are from the Robinson form. Fiedler (1975) demonstrated that the spectral method, that is, ordering using the Fiedler vector, will find the Robinson form if such exists. If such an order does not exist, then the spectral order can be viewed as an algorithm for placing objects on a straight line with the goal of minimizing pair-wise distance weighted by the similarity matrix. Note that the goal here is to find the locations on a line, the order of the objects comes as a by-product. Optimizing stress functions that are directly based on the order of the similarity matrix is typically a NP-hard problem. In such cases the problem typically resembles a travelling salesman problem. Consequently, heuristics to discover a good order are variants of TSP solvers (Johnson et al., 2004).

Ordering data rows and columns and using this order to discover an underlying structure has been suggested. Gionis et al. (2004) and Tatti and Vreeken (2012) suggested discovering tile hierarchies from ordered binary datasets. Both work used Bernoulli models, a log-linear model, to score the hierarchies.

In our experiments, the number of borders is small. Hence, we are able to solve Problem 6 exactly using a dynamic program. This program requires quadratic time. If the number of borders becomes too large, that is, close to the number of entries, using an exact solver becomes impractical. In order to overcome this problem we can use heuristic approaches, such as top-down approaches where we select greedily a new change-point (see Shatkay and Zdonik (1996); Bernaola-Galván et al. (1996); Douglas and Peucker (1973); Lavrenko et al. (2000), for example) or bottom-up approaches where at the beginning each point is a segment, and points are combined in a greedy fashion (see Palpanas et al. (2004), for example). Yet another option is a randomized heuristic was suggested by Himberg et al. (2001), where we start from a random segmentation and optimize the segment boundaries. These approaches, although fast, are heuristics and have no theoretical guarantees of the approximation quality. A divide-and-segment approach, an approximation algorithm with theoretical guarantees on the approximation quality was given by Terzi and Tsaparas (2006).

7 Experiments

Setup:

In our experiments we used 6 real-world datasets and one synthetic data. The first two graphs, DblpCP and DblpCF, are ego-networks of Christos Papadimitriou and Christos Faloutsos, that is, the graphs contain the collaborators obtained from DBLP,​333http://www.informatik.uni-trier.de/~ley/db/ two researchers are connected if they have a joint paper. The other two datasets, Fb107 and Fb1912, were the two largest ego-networks taken from the Facebook dataset obtained from SNAP repository.​444http://snap.stanford.edu/snap/ The Mammals presence data consists of presence records of European mammals within geographical areas of $50\times 50$ kilometers (Mitchell-Jones et al., 1999).​555The full version of the Mammals dataset is available for research purposes from the Societas Europaea Mammalogica at http://www.european-mammals.org We computed a similarity matrix between different locations, say $i$ and $j$, by considering the observed joint number of different mammals, normalized by $c_{i}c_{j}$, where $c_{i}$ and $c_{j}$ are the number of mammals observed at $i$ and $j$, respectively. Dataset Paleo666NOW public release 030717 available from Fortelius (2005). contains information of species fossils found in specific paleontological sites in Europe Fortelius (2005). We constructed a similarity matrix between two fossils, say $a$ and $b$, by computing the number of sites in which both $a$ and $b$ have been discovered. We also created a synthetic dataset with $250\,000$ vertices and $279\,223$ edges where edges had a higher probability of being closer to the diagonal. The basic characteristics of the datasets are given in Table 1.

We ordered the vertices of each graph (except synthetic) according to Fiedler’s vector using algorithm given by Atkins et al. (1999). Since Atkins’ algorithm may produce several orders, this happens when the graph has disconnected components, we pick one order at random. Once the graph is ordered we apply our heuristic method given in Section 4.3 along with the refinement step described in Section 5. We also applied the algorithm for discovering for exact borders described in Section 4.2. We used Bernoulli model for unweighted graphs, Poisson model for Paleo data and $L_{2}$ error for Mammals. We set the number of bands to $3$ in DBLP and Paleo datasets and $4$ in Facebook and Mammals datasets. The statistics of the datasets and the experiments are given in Table 1. Due to the slow convergence of heuristic method in Mammals and Synthetic, we limited the number of iterations to $2000$ and did not apply refinement. We also noticed that after each random tie-breaker the score gets a (relatively small) bump. Hence, we force random tie-breaker after each $50$th iteration in these two datasets.

Results:

Let us first focus on computational complexity. We see that the exact algorithm works very well in practice, it is faster than the heuristic discovery, and it has the benefit of producing the exact result. The main reason for this is that for medium-sized datasets the running time for the exact algorithm is very competitive and the slowness of the heuristic algorithm is due to high number of iterations it requires to converge. This dynamic changes when the size of the data increases. For large datasets, the fact that we need to store the graph in the full form will slow the algorithm down (or simply run out of memory). This is the case with the Synthetic data where we could not run the exact algorithm due to the memory limitations.

Unlike the exact method, the heuristic discovery of bands can use sparsity to its advantage (see Section 4.4). More specifically, one iteration is linear w.r.t. the number of edges in the lattice constructed from the edges of the original graph. While it is possible in theory that the number of lattice edges is significantly higher than the number of edges in the original graph, third column in Table 1 demonstrates that in our experiments the number of lattice edges is about 2–3 times the number of edges in the original graph. Due to this property the running times stay reasonably small, minutes at worst.

The number of discovered borders, 5th column, is small, which in turns implies that the final segmentation is fast, and indeed we spend most of the time computing the borders.

Let us now study convergence of the heuristic in more detail. Table 2 shows that the heuristic converges close to the exact values w.r.t. the score function. The number of iterations needed in a single refinement round, given in $3$th and $4$th columns of Table 1, grows as the graph gets larger. The number of random tie-breaker iterations is about $10$ percent. In order to study the convergence in more details, we plotted $s\mathopen{}\left(\mathcal{B}\right)$ in Figure 5, where $\mathcal{B}$ is the current set of approximate borders, as a function of rounds left to convergence. Since the initial border set comes from a random monotonic order, we repeated this experiment 10 times. Note that $s\mathopen{}\left(\mathcal{B}\right)$ is an upper bound for score of the final segmentation. From the results we see that initial orders have a bad score, and there is a high variance in the initial score and the convergence time. However, the score stabilizes quickly, close to its final value, and the final score has little dependence of the starting point. Majority of the rounds is spent in cosmetic improvement. This hints that if the convergence time is a factor, we may stop the border discovery early and still get a good score. This is demonstrated with the Synthetic dataset. The convergence of the score is given in Figure 5 and again we see the same behaviour where the initial order has a bad score but quickly obtains a stable score. Running the heuristic took $37$ minutes and we were not able to run the exact algorithm due to the memory limitations.

Let us next look at the refinement. In our experiments, the number of refinement iterations is low, around 10. Moreover, the improvement in score is modest, around $1$%. This is also the case when we apply to the datasets that are not ordered with the fiedler vector. This suggests that the refinement procedure should not be used alone but instead it should be used only to improve the already good order. We also tried discovering bands using random orders. Naturally, the score for discovered bands is weaker. The gain of using spectral orders depends how strongly banded is the data. For example, Paleo contains a strong banded structure, and the score using a random order was significantly worse than the spectral order. Finally, using the exact border discovery with the refinement may produce a worse score than using the approximate border discovery, as shown by the scores for Fb1912.

Finally, let us look at the discovered bands given in Figures 6–8. The first band in the DBLP datasets is relatively small, typically discovering small clusters that are due to joint publication with the author of the ego-network. Facebook graphs contain a strong clustering structure which is seen in the discovered bands. However, we also see the overlap effect: in Fb1912 the clusters overlap. In Figure 8 we provide small snipets from DBLP graphs along with the authors. For example, in DblpCP, the inner band contains authors that have collaborate with each other but these authors do not form a clique. In Paleo, the bands correspond to the weights of the edges: the inner band has edges with higher weights while the outer bands have less edges. In Mammals, the spectral order—which was computed without the knowledge of the geological location—corresponds roughly to the latitude: one extreme being Scandinavia while on the other extreme being Greece and Crete. The order fails to capture Spain and Portugal suggesting that the data cannot be explained fully by a single order. The discovered bands, especially in the northern areas, suggest a continuous change in fauna as we move along the latitude.

8 Conclusion

In this paper we introduced an optimization problem in order to model the concentration of edges next to the diagonal. More specifically, given a graph and an integer $K$, our goal is to order the vertices and find $K$ bands around the diagonal, inner bands being more dense. As a measure of goodness we use the likelihood of a log-linear model, a family containing many standard distributions.

We divide the problem into two subproblems. The first problem is to find a good order while the second problem is to discover the bands. We introduced two solvers for the latter problem. The first solver is exact and based on isotonic regression, the second solver is a heuristic approach that sorts the entries into a linear order and applies one-dimensional isotonic regression. By doing so, the solver can exploit sparsity of the graph. Both approaches complement each other. If we can store the graph as a full adjacency matrix, then it is beneficial to use the exact algorithm. The algorithm runs in $O(N^{4})$ time, at worst, but is closer to $O(N^{2}\log P)$, where $N$ is the number of vertices and $P$ is the number of borders. On the other hand, if the number of vertices is large but the graph is sparse it is better to apply the heuristic approach with a limited number of iterations. A single iteration runs in $O({\left|E(L)\right|})$ time, where $L$ is the lattice describing the relations between the nodes (described in Section 4.4). At worst, ${\left|E(L)\right|}$ is $O(N^{2})$ but in practice it is smaller for sparse graphs.

We discover the order mainly by using Fiedler order. It would be fruitful to develop algorithms that induce an order while looking for the optimal segmentation simultaneously. The work done by Mannila and Terzi (2007) with nested datasets hints that this problem is NP-hard. However, it may be possible to develop an efficient heuristic approach.

Another open question, which we will leave to a future work, is the problem of choosing $K$, the number of segments. This question is not only specific to our problem but also occurs in many classic problems such as clustering or sequence segmentation. In our experiments computing the final segmentation from the set of borders is cheap. Hence, instead of studying a single segmentation we can compute several segmentations simultaneously and present all of them to the user. Alternatively, since our score is essentially a log-likelihood, we can design an MDL approach to select a segmentation with a good score.

Appendix A Proof of Proposition 1

In order to prove this proposition, we first need to establish a series of lemmas. We will use the following fact. Let $B$ be a set of entries, and let $C\subsetneq B$ be a non-empty subset of $B$ such that $\mathit{a}\mathopen{}\left(C\right)<\mathit{a}\mathopen{}\left(B\right)$, then $\mathit{a}\mathopen{}\left(B\setminus C\right)>\mathit{a}\mathopen{}\left(B\right)$.

The first lemma states that if you cut a segment from a maximal corner, then that segment will be more dense. Conversely, any segment adjacent to a maximal corner will be more sparse.

Lemma 1

Let $U$ be a corner. Let $V=\mathit{right}\mathopen{}\left(U\right)$ be a corner. Let $X\supseteq U$ be a corner. If $V\setminus X\neq\emptyset$, then $\mathit{a}\mathopen{}\left(V\setminus X\right)\geq\mathit{a}\mathopen{}\left(V\setminus U\right)$. If $X\setminus V\neq\emptyset$, then $\mathit{a}\mathopen{}\left(X\setminus V\right)<\mathit{a}\mathopen{}\left(V\setminus U\right)$.

Proof

Assume that $\mathit{a}\mathopen{}\left(V\setminus X\right)<\mathit{a}\mathopen{}\left(V\setminus U\right)$. Then, since

[TABLE]

$\mathit{a}\mathopen{}\left((V\cap X)\setminus U\right)>\mathit{a}\mathopen{}\left(V\setminus U\right)$, which contradicts the definition of $V$. The second case holds because otherwise $\mathit{a}\mathopen{}\left((V\cup X)\setminus U\right)\geq\mathit{a}\mathopen{}\left(V\setminus U\right)$, which contradicts the definition of $V$. ∎

The next lemma essentially states that if a corner cuts a maximal corner, then it cannot be a border.

Lemma 2

Let $U$ be a corner. Write $V=\mathit{right}\mathopen{}\left(U\right)$. Let $X\supsetneq U$ be a corner such that $V\setminus X\neq\emptyset$, then $\mathit{a}\mathopen{}\left((V\cup X)\setminus X\right)\geq\mathit{a}\mathopen{}\left(X\setminus U\right)$. Consequently, $X$ is not a border.

Proof

By definition of $V$, $\mathit{a}\mathopen{}\left(V\setminus U\right)\geq\mathit{a}\mathopen{}\left(X\setminus U\right)$.

Lemma 1 implies that $\mathit{a}\mathopen{}\left(V\setminus X\right)\geq\mathit{a}\mathopen{}\left(V\setminus U\right)$. As $(V\cup X)\setminus X=V\setminus X$, we can combine the two inequalities and prove the first claim. $X$ cannot be a border since $U\subsetneq X\subsetneq V\cup X$. ∎

Lemmas 1–2 can be modified to hold for minimal corners. The proofs for these lemmas are similar to the proofs of Lemmas 1–2.

Lemma 3

Let $U$ be a corner. Let $V=\mathit{left}\mathopen{}\left(U\right)$ be a corner. Let $X\subsetneq U$ be a corner. If $X\setminus V\neq\emptyset$, then $\mathit{a}\mathopen{}\left(X\setminus V\right)\leq\mathit{a}\mathopen{}\left(U\setminus V\right)$. If $V\setminus X\neq\emptyset$, then $\mathit{a}\mathopen{}\left(V\setminus X\right)>\mathit{a}\mathopen{}\left(U\setminus V\right)$.

Lemma 4

Let $U$ be a corner. Write $V=\mathit{left}\mathopen{}\left(U\right)$. Let $X\subsetneq U$ be a corner such that $X\setminus V\neq\emptyset$, then $\mathit{a}\mathopen{}\left(X\setminus(V\cap X)\right)\leq\mathit{a}\mathopen{}\left(U\setminus X\right)$. Consequently, $X$ is not a border.

Assume that we are given two $K$-segmentations, $\mathcal{U}=U_{0},\ldots,U_{K}$ and $\mathcal{V}=V_{0},\ldots,V_{K}$. We write $\mathcal{U}\precneqq\mathcal{V}$ if there exists $i$ such that $U_{j}=V_{j}$ for $j<i$, and $U_{i}\subsetneq V_{i}$. The idea behind the proof of Proposition 1 is that we can replace a non-border in a segmentation, say $\mathcal{U}$, either with its minimal or maximal corner such that either the score will increase or that the resulting segmentation, say $\mathcal{V}$, is smaller wrt. partial order, $\mathcal{V}\precneqq\mathcal{U}$.

Next lemma makes sure that we can always pick a non-border in a segmentation that such that its maximal and minimal corners form a chain with the other corners in the segmentation.

Lemma 5

Let $U_{0},\ldots,U_{K}$ be a $K$-segmentation. Assume that, say, $U_{i}$ is not a border. Then there exists $U_{j}$ with $0<j<K$, and two corners $X$ and $Y$ such that $U_{j-1}\subseteq X\subsetneq U_{j}\subsetneq Y\subseteq U_{j+1}$ and $\mathit{a}\mathopen{}\left(Y\setminus U_{j}\right)\geq\mathit{a}\mathopen{}\left(U_{j}\setminus X\right)$.

Proof

Since $U_{i}$ is not a border, there exist $X$ and $Y$ such that $X\subsetneq U_{i}\subsetneq Y$ and $\mathit{a}\mathopen{}\left(Y\setminus U_{i}\right)\geq\mathit{a}\mathopen{}\left(U_{i}\setminus X\right)$. We need to show that $U_{i-1}\subseteq X$ and $Y\subseteq U_{i+1}$, or possibly modify $X$, $Y$, and $i$ such that inclusions hold.

We can safely assume that $Y=\mathit{right}\mathopen{}\left(U_{i}\right)$.

Assume that $Y\nsubseteq U_{i+1}$, that is, $Y\setminus U_{i+1}\neq\emptyset$. Lemma 2 implies that $\mathit{a}\mathopen{}\left(Y\setminus U_{i+1}\right)\geq\mathit{a}\mathopen{}\left(U_{i+1}\setminus U_{i}\right)$. Hence, we can redefine $X=U_{i}$ and $Y=\mathit{right}\mathopen{}\left(U_{i+1}\right)$, and then increase $i$ by one and repeat the argument. This process will eventually stop since we increase $i$ every step. When we finally stop, note that both inclusions will hold, and we have proved the lemma.

If we start with the case where $Y\subseteq U_{i}$ and $U_{i-1}\nsubseteq X$, we proceed to modify $X$, $Y$, and $i$ in the opposite direction. ∎

Our next goal is Corollary 2 which will make sure that we can always make sure that the segmentation is monotonic. For that we need the following two lemmas, that follow immediately from the fact that our scoring function is a log-linear model.

Lemma 6

$s\mathopen{}\left(X\mid u\right)$* is a concave function of $u$. Let $U$ and $V$ be two segments such that $\mathit{a}\mathopen{}\left(U\right)\leq\mathit{a}\mathopen{}\left(V\right)$. Let $r=\arg\max_{u}s\mathopen{}\left(U\mid u\right)$ and $t=\arg\max_{u}s\mathopen{}\left(V\mid u\right)$. Then $r\leq t$.*

Proof

Let $X$ be a segment. A straightforward calculation shows that $\partial s\mathopen{}\left(X\mid u\right)/\partial u={\left|X\right|}(\mathit{a}\mathopen{}\left(X\right)-\operatorname{E}_{u}\mathopen{}\left[x\right])$, where $\operatorname{E}_{u}\mathopen{}\left[x\right]=\sum_{x}xp(x\mid u)$ is the mean of the log-linear model.

We also have $\partial\operatorname{E}_{u}\mathopen{}\left[x\right]/\partial u=\operatorname{Var}_{u}\mathopen{}\left[x\right]\geq 0$. This immediately proves that $s\mathopen{}\left(X\mid u\right)$ is a concave function. The optimal (possibly infinite) value is reached when $\mathit{a}\mathopen{}\left(X\right)=\operatorname{E}_{u}\mathopen{}\left[x\right]$. Since $\operatorname{E}_{u}\mathopen{}\left[x\right]$ is a monotone function of $u$, optimal value for $s\mathopen{}\left(U\mid u\right)$ will be smaller or equal than $s\mathopen{}\left(V\mid u\right)$. ∎

Lemma 7

Let $U$ and $V$ be two segments such that $\mathit{a}\mathopen{}\left(U\right)<\mathit{a}\mathopen{}\left(V\right)$. Let $r$ and $t$ be two parameters, $r\geq t$. Then there exists $u$, $r\geq u\geq t$ such that $s\mathopen{}\left(U\mid r\right)+s\mathopen{}\left(V\mid t\right)\leq s\mathopen{}\left(U\mid u\right)+s\mathopen{}\left(V\mid u\right)$.

Proof

Let $r^{*}$ be such that $s\mathopen{}\left(U\mid r^{*}\right)$ is optimal, and define $t^{*}$ similarly. Note that $r^{*}$ or $t^{*}$ can be infinite. Lemma 6 implies that $r^{*}\leq t^{*}$.

Assume that $r\leq t^{*}$. This implies that $t\leq r\leq t^{*}$. Since the score function is concave, $s\mathopen{}\left(V\mid r\right)\geq s\mathopen{}\left(V\mid t\right)$. Set $u=r$ to prove the lemma.

Assume that $r>t^{*}$. Set $u=\max\mathopen{}\left(t^{*},t\right)$. Since $r\geq u\geq r^{*}$, due to concavity of $s\mathopen{}\left(U\mid u\right)\geq s\mathopen{}\left(U\mid r\right)$ and by definition $s\mathopen{}\left(V\mid u\right)\geq s\mathopen{}\left(V\mid t\right)$. This proves the lemma. ∎

Corollary 2

Let $\mathcal{U}$ be a $K$-segmentation. Let $r_{1},\ldots,r_{K}$ be $K$ parameters such that $r_{i}\geq r_{i-1}$. There exists a monotone $K$-segmentation $\mathcal{V}\preceq\mathcal{U}$ such that $s\mathopen{}\left(\mathcal{V}\right)\geq s\mathopen{}\left(\mathcal{U}\mid r_{1},\ldots,r_{K}\right)$.

Proof

Assume that $\mathcal{U}$ is not monotone, that is, there exists $U_{i}$ and $U_{i+1}$ such that $\mathit{a}\mathopen{}\left(U_{i}\setminus U_{i-1}\right)<\mathit{a}\mathopen{}\left(U_{i+1}\setminus U_{i}\right)$. Lemma 7 implies that we can replace $r_{i}$ and $r_{i+1}$ with a common parameter, say $u$, and not decrease the score. We can now join the $i$th and $i+1$th segments into one segment with a parameter $u$ and add a new empty corner at the beginning of the segmentation (with an infinitely large parameter). This modification makes the segmentation smaller w.r.t our order. We repeat this step until convergence. When converged, we have obtained a monotone $K$-segmentation whose score is at least as good as the original segmentation. ∎

We will now show that we can modify a segmentation containing a non-border.

Proposition 7

Let $U_{0},\ldots,U_{K}=\mathcal{U}$ be a monotone $K$-segmentation. Assume that there exist $U_{i}$, $X$, and $Z$ such that $U_{i-1}\subseteq X\subsetneq U_{i}\subsetneq Z\subseteq U_{i+1}$ and $\mathit{a}\mathopen{}\left(Z\setminus U_{i}\right)\geq\mathit{a}\mathopen{}\left(U_{i}\setminus X\right)$.

Then there exists a monotone $K$-segmentation $\mathcal{V}$ such that either $\mathcal{V}\precneqq\mathcal{U}$ and $s\mathopen{}\left(\mathcal{V}\right)\geq s\mathopen{}\left(\mathcal{U}\right)$ or $s\mathopen{}\left(\mathcal{V}\right)>s\mathopen{}\left(\mathcal{U}\right)$.

In order to prove the proposition we will introduce some helpful notation. First, given two parameters $r$ and $t$, we define

[TABLE]

We also define

[TABLE]

This function is essentially the difference between two scores.

Lemma 8

Let $U_{i-1}\subseteq X\subsetneq Y\subseteq U_{i+1}$. We have $h(Y;r,t)-h(X;r,t)=g({\left|Y\setminus X\right|},\mathit{a}\mathopen{}\left(Y\setminus X\right);r,t)$.

Proof

Let us define $\mathit{cs}\mathopen{}\left(Z\right)=\sum_{z\in Z}D(z)$. Note that

[TABLE]

where const does not depend on $X$. Let us write $d=Z(r)-Z(t)$ and $Z=Y\setminus X$. This allows us to write

[TABLE]

This completes the proof. ∎

Proof (of Proposition 7)

Replace $U_{i}$ with $Z$ and let $\mathcal{U}^{*}$ be the resulting segmentation. Similarly, replace $U_{i}$ with $X$ and let $\mathcal{U}^{\prime}$ be the resulting segmentation. Define $y=\sup_{r,t}h(U_{i};r,t)$.

Fix $\epsilon>0$. Lemma 6 implies that there exist $r_{1},\ldots,r_{K}$ s.t.

[TABLE]

From now on we will write $h(X)$ to mean $h(X;r_{i},r_{i+1})$ and $g(k,\delta)$ to mean $g(k,\delta;r_{i},r_{i+1})$. Note that we must have $h(U_{i})\geq y-\epsilon$.

Assume that $h(Z)>y+\epsilon$. Then $s\mathopen{}\left(\mathcal{U}^{*}\mid r_{1},\ldots,r_{K}\right)>s\mathopen{}\left(\mathcal{U}\mid r_{1},\ldots,r_{K}\right)+\epsilon\geq s\mathopen{}\left(\mathcal{U}\right)$. Corollary 2 applied to $\mathcal{U}^{*}$ shows that there is a monotone segmentation with a score better than $s\mathopen{}\left(\mathcal{U}\right)$.

Assume that $h(Z)\leq y+\epsilon$. We must have $h(U_{i})+\epsilon\geq y\geq h(Z)-\epsilon$ or, equivalently, $2\epsilon\geq h(Z)-h(U_{i})$.

Define $\beta=\mathit{a}\mathopen{}\left(Z\setminus U_{i}\right)$ and $\alpha=\mathit{a}\mathopen{}\left(U_{i}\setminus X\right)$, $n={\left|Z\setminus U_{i}\right|}$, $m={\left|U_{i}\setminus X\right|}$. Define $c=n/m$. We now have

[TABLE]

which implies $y-h(X)\leq\epsilon(1+2c^{-1})\leq\epsilon(1+2{\left|U_{K}\right|})$. Corollary 2 now implies that there exists a monotone segmentation $\mathcal{V}$ with $\mathcal{V}\precneqq\mathcal{U}$ such that $y-s\mathopen{}\left(\mathcal{V}\right)\leq\epsilon(1+2{\left|U_{K}\right|})$. Since this holds for any $\epsilon>0$, we have proved the proposition. ∎

Proof (of Proposition 1)

Assume that $U_{j}$ is not a border, then Lemma 5 implies that there exist $U_{i}$, $X$, and $Y$ such that the conditions in Proposition 7 are satisfied. Apply Proposition 7 to obtain a new monotone segmentation, $\mathcal{V}$. Reapply the step to $\mathcal{V}$ until $\mathcal{V}$ consists only of borders. This procedure terminates since at each step we either increase score or move segmentation is moved to the left w.r.t the partial order $\prec$. There are finite number of segmentations and no segmentation is visited twice, hence we converge to a segmentation consisting only of borders. ∎

Appendix B Proof of Proposition 3

Proof (of Proposition 3)

Define $Z=\mathit{right}\mathopen{}\left(U\right)$. Let us first prove that $Z$ is a border.

Let $X=\mathit{left}\mathopen{}\left(Z\right)$. If $U\setminus X\neq\emptyset$, then Lemma 4 implies that $U$ is not a border. Hence $U\subseteq X$. Lemma 1 implies that $\mathit{a}\mathopen{}\left(Z\setminus X\right)\geq\mathit{a}\mathopen{}\left(Z\setminus U\right)$.

Let $Y\supsetneq Z$ be a corner, then Lemma 1 implies that $\mathit{a}\mathopen{}\left(Y\setminus Z\right)<\mathit{a}\mathopen{}\left(Z\setminus U\right)\leq\mathit{a}\mathopen{}\left(Z\setminus X\right)$. By definition, $Z\setminus X$ has the smallest possible average. Consequently, $Z$ is a border.

Assume that $Z\setminus V\neq\emptyset$, then Lemma 2 implies that $V$ is not a border, which is a contradiction. Hence $Z\subseteq V$. Since $V$ is the border next to $U$, we must have $Z=V$.

The proof in other direction is similar. ∎

Appendix C Proof of Proposition 6

Proof

Let $k$ be as in Proposition 5. Fix $U_{j}\in\mathit{brd}\mathopen{}\left(T^{k}\right)$ with $j>0$. Let us write $V^{i}$ to be the portion of $T^{i}$ that corresponds to the entries in $C=U_{j}\setminus U_{j-1}$. Since $\mathit{brd}\mathopen{}\left(T^{i}\right)=\mathit{brd}\mathopen{}\left(T^{k}\right)$ for $i>k$, it is enough to prove the result by showing that there exists $m$ such that $V^{m}=V^{m+2}$. We will prove the result by induction over the size of $C$.

To that end, consider a DAG $G$ where the nodes are the entries in $C$. Two distinct nodes $(a,b)$, $(c,d)$ are connected if and only if $a\leq c$ and $b\leq d$. Let $p$ be a sink in $G$. Define $V^{k}_{p}$ to be equal to $V^{k}$ without the entry $p$ and let $V^{i}_{p}$, for $i>k$, be the order obtained from $V^{i-1}_{p}$ by simulating FindOrder with $w_{F}$. Since $p$ is a sink, it does not block any entries as FindOrder updates the order. This implies that $V^{i}_{p}$ is equal to $V^{i}$ with $p$ deleted for any $i\geq k$.

By the induction assumption there exists $m_{p}$ such that $V^{m_{p}}_{p}=V^{m_{p}+2}_{p}$ for each sink $p$. Define $m=\max m_{p}$. Note that $V^{m}_{p}=V^{m+2}_{p}$ for each sink $p$.

Assume now that we have only one sink, say $p$. Then $p$ will always be last entry in $V^{i}$ and it follows immediately that $V^{m}=V^{m+2}$. We can safely assume that we have more than one sink.

Assume that we have exactly two sinks, say $p$ and $q$. Assume that $V^{m}_{p}$ does not end on $q$. Then it must be that $V^{m}_{p}$ ends on a parent of $p$ and $p$ is the last entry in $V^{m}$. Since $V^{m}_{p}=V^{m+2}_{p}$, the proposition follows. Similar argument holds for $V^{m}_{q}$, $V^{m+1}_{p}$, and $V^{m+1}_{q}$. Hence, we can safely assume that $V^{m}_{p}$, $V^{m}_{q}$, $V^{m+1}_{p}$, and $V^{m+1}_{q}$ all have $p$ or $q$ as their last entry. This implies that $V^{m}$, $V^{m+1}$, and $V^{m+2}$. have $p$ and $q$ as their last entries. Assume that $p$ occurs before $q$ in $V^{m}$. There will be a point during $\textsc{FindOrder}(w_{F}(\cdot;T^{m}))$ when $p$ and $q$ are in the heap $H$ at the same time, otherwise there will be an entry between $p$ and $q$ in $V^{m+1}$. This implies that $q$ occurs before $p$ in $V^{m+1}$. We apply the same argument to $V^{m+1}$ to conclude that $p$ occurs before $q$ in $V^{m+2}$. This implies that $V^{m+1}=V^{m+2}$.

Assume that we have more than two sinks. Let $p$, $q$, and $r$ be three sinks. We can deduce from $V^{m}_{r}$ whether $p$ occurs before $q$, or vice versa. Since we can do this for any sink triplet, we can deduce the order of sinks in $V^{m}$. The last sink, say $p$, will be the last in $V^{m}$. By definition of $m$, the order will be the same in $V^{m+2}$ and the same sink will be also the last in $V^{m+2}$. Since $V^{m}_{p}=V^{m+2}_{p}$, the proposition follows. ∎