1

A Generalized Framework for Kullback-Leibler
Markov Aggregation

Abstract—This paper proposes an information-theoretic cost
function for aggregating a Markov chain via a (possibly stochastic) mapping. The cost function is motivated by two objectives:
1) The process obtained by observing the Markov chain through
the mapping should be close to a Markov chain, and 2) the
aggregated Markov chain should retain as much of the temporal
dependence structure of the original Markov chain as possible.
We discuss properties of this parameterized cost function and
show that it contains the cost functions previously proposed by
Deng et al., Xu et al., and Geiger et al. as special cases. We
moreover discuss these special cases providing a better understanding and highlighting potential shortcomings: For example,
the cost function proposed by Geiger et al. is tightly connected
to approximate probabilistic bisimulation, but leads to trivial
solutions if optimized without regularization. We furthermore
propose a simple heuristic to optimize our cost function for
deterministic aggregations and illustrate its performance on a
set of synthetic examples.
Index Terms—Markov chain, lumpability, predictability, bisimulation, model reduction

I. I NTRODUCTION
Markov aggregation is the task of representing a Markov
chain with a large alphabet by a Markov chain with a smaller
alphabet, thus reducing model complexity while at the same
time retaining the computationally and analytically desirable
Markov property (see Fig. 1). Such a model reduction is
necessary if the original Markov chain is too large to admit
simulation, estimating model parameters from data, or control
(in the case of Markov decision processes). These situations
occur often in computational chemistry (where aggregation is
called coarse-graining, e.g., [1]), natural language processing,
and the simulation and control of large systems (giving rise
to the notion of bisimulation, e.g., [2]). Additionally, Markov
aggregation can be used as a tool in exploratory data analysis,
either to discover groups of “similar” states of a stochastic
process or to cluster data points, cf. [3], [4].
Information-theoretic cost functions were proposed for
Markov aggregation in [5]–[8]. Specifically, the authors of [5]
proposed a cost function linked to the predictability of the
aggregated Markov chain. Such an approach is justified if the
original model is nearly completely decomposable, i.e., if there
is a partition of the alphabet such that transitions within each
element of the partition occur quickly and randomly, while
transitions between elements of the partition occur only rarely.
Building on this work, the authors of [7] proposed a cost
function linked to lumpability, i.e., to the phenomenon where
a function of a Markov chain is Markov. Such an approach

X ∼ Mar(X , P)

Observation pY |X

arXiv:1709.05907v1 [] 18 Sep 2017

Rana Ali Amjad, Clemens Blöchl, Bernhard C. Geiger
Institute for Communications Engineering, Technical University of Munich, Germany
ranaali.amjad@tum.de, bernhard.geiger@tum.de

Ag

gr

eg
a

tio
n

D̄(Y||Ỹ)
Y

Ỹ ∼ Mar(Y, Q)

Fig. 1. Illustration of the aggregation problem: A stationary first-order Markov
chain X is given. We are interested in finding a conditional distribution pY |X
and an aggregation of X, i.e., a Markov chain Ỹ on Y. The conditional
distribution pY |X defines a stationary process Y, a noisy observation of X.
Y might not be Markov of any order, but can be approximated by a Markov
chain Ỹ.

is justified whenever there are groups of states with similar
probabilistic properties (in a well-defined sense). Both [5]
and [7] focus on deterministic aggregations, i.e., every state
of the original alphabet is mapped to exactly one state of the
reduced alphabet. Moreover, the authors of both references
arrive at their cost functions by lifting the aggregated Markov
chain to the original alphabet. The authors of [6] present an
information-theoretic cost function for stochastic aggregations,
but they do not justify their choice by an operational characterization (such as predictability or lumpability). Instead, they
arrive at their cost function via the composite of the original
and the aggregated Markov chain.
In this paper, we extend the works [5]–[8] as follows:
1) We present a two-step approach to Markov aggregation (Section III): Observing the original Markov chain
through a (stochastic or deterministic) channel, and
then approximating this (typically non-Markov) process
as a Markov chain (see Fig. 1). This approach has
already been taken by [7], albeit only for deterministic
aggregations.
2) Using this two-step approach, we propose a parameterized, information-theoretic cost function for Markov
aggregation (Section IV). We arrive at this cost function
neither via lifting nor via the composite model, but via
requiring specific operational qualities of the process
observed through the channel: It should be close to a
Markov chain and it should retain the temporal dependence structure of the original Markov chain.

2

3) We show that our cost function contains the cost functions of [5]–[8] as special cases (Section V). We also
discuss previous algorithmic approaches to the Markov
aggregation problem.
4) We propose a simple, low-complexity heuristic to minimize our generalized cost function for deterministic
aggregations (Section VI).
5) As a side result, we justify the cost function proposed
in [7] by showing a tight connection to approximate
probabilistic bisimulation (Section III-A).
We illustrate our cost function for various examples in
Section VII. Specifically, we investigate the aggregation of
quasi-lumpable and nearly completely decomposable Markov
chains, and we look at a toy example from natural language
processing. We also take up the approach of [3], [4] to perform
clustering via Markov aggregation. In future work, we shall
extend our efforts to Markov decision processes, and provide
a theory for lifting stochastic aggregations as indicated in [6,
Remark 3].
II. N OTATION AND D EFINITIONS

1

(1)

The Markov chain is time-homogeneous if, the right-hand side
of (1) does not depend on n, i.e., if
pZn |Zn−1 (zn |zn−1 ) = pZ2 |Z1 (zn |zn−1 ) := Pzn−1 →zn .

z∈Z

H(Z2 |Z1 ) :=

X

H(Z2 |Z1 = z)pZ1 (z)

I(Z1 ; Z2 ) := H(Z2 ) − H(Z2 |Z1 ).

(3c)

The entropy rate and redundancy rate of a stationary stochastic
process Z (not necessarily Markov) are
H(Z1n )
= lim H(Zn |Z1n−1 )
n→∞
n→∞
n
R̄(Z) := lim I(Zn ; Z1n−1 ) = H(Z) − H̄(Z).

H̄(Z) := lim

n→∞

(3d)
(3e)

The Kullback-Leibler divergence rate (KLDR) between two
stationary stochastic processes Z and Z′ on the same finite
alphabet Z is [10, Ch. 10]
pZ ′n1 (z1n )
1 X
pZ ′n1 (z1n ) log
n→∞ n
pZ1n (z1n )
n
n

(3f)

z1 ∈Z

provided the limit exists. If the limit exists, it is finite if, for
all n and all z1n , pZ1n (z1n ) = 0 implies pZ ′ n1 (z1n ) = 0 (short:
pZ ′ n1 ≪ pZ1n ). In particular, if Z′ ∼ Mar(Z, P′ ) and Z ∼
Mar(Z, P), then [11]
D̄(Z′ ||Z) =

X

′
µz Pz→z
′ log

z,z ′ ∈Z

′
Pz→z
′
.
Pz→z′

(3g)

provided P ≪ P′ .
These information-theoretic quantities can be used to give
an equivalent definition of Markovity:
Lemma 1 ([12, Prop. 3]). Suppose the stochastic process Z
is stationary. Then, Z is Markov iff H̄(Z) = H(Z2 |Z1 ).
If Z is a stationary process on Z (not necessarily Markov),
then one can approximate this process by a Markov chain
Z̃ ∼ Mar(Z, P):
Lemma 2 ([10, Cor. 10.4]). Let Z be a stationary process
on Z, and let Z′ ∼ Mar(Z, P′ ) be any Markov chain on Z.
Then,
P = arg min D̄(Z||Z′ )

(2)

If the transition probability matrix P = [Pzn−1 →zn ] of a
time-homogeneous Markov chain is irreducible and aperiodic
(see [9] for terminology), then there exists a unique vector µ
such that µT = µT P, which represents the invariant distribution of Z. If the initial distribution pZ1 coincides with µ,
then Z is stationary and we denote this stationary, irreducible
and aperiodic first-order Markov chain by Z ∼ Mar(Z, P).
In this work we deal exclusively with first-order stationary,
irreducible and aperiodic Markov chains.
We use information-theoretic cost functions for the aggregation problem. The entropy of Z, the conditional entropy of

(3b)

z∈Z

D̄(Z′ ||Z) := lim

We denote vectors and matrices by bold lower case and
blackboard bold upper case letters, e.g., a and A. A diagonal
matrix with vector a on the main diagonal is denoted by
diag(a). The transpose of A is AT .
We denote random variables (RVs) by upper case letters,
e.g., Z, and their alphabet by calligraphic letters, e.g., Z.
In this work we will restrict ourselves to RVs with finite
alphabets, i.e., |Z| < ∞. Realizations are denoted by lower
case letters, e.g., z, where z ∈ Z. The probability mass
function (PMF) of Z is pZ , where pZ (z) := Pr(Z = z) for all
z ∈ Z. Joint and conditional PMFs are defined accordingly.
We denote a one-sided, discrete-time, stochastic process
with Z := (Z1 , Z2 , . . . ), where each Zk takes values from the
n
:= (Zm , Zm+1 , Zn ).
(finite) alphabet Z. We abbreviate Zm
We consider only stationary processes, i.e., PMFs are invariant
w.r.t. a time shift. In particular, the marginal distribution of Zk
is equal for all k and shall be denoted as pZ .
A first-order Markov chain is a process that satisfies, for all
n > 1 and all z1n ∈ Z n ,
pZn |Z n−1 (zn |z1n−1 ) = pZn |Zn−1 (zn |zn−1 ).

Z2 given Z1 , and the mutual information between Z1 and Z2
are defined by
X
H(Z) := −
pZ (z) log pZ (z)
(3a)

(4a)

P′

where
Pz→z′ = pZ2 |Z1 (z ′ |z).

(4b)

Moreover, for Z̃ ∼ Mar(Z, P),
D̄(Z||Z̃) = H(Z2 |Z1 ) − H̄(Z).

(4c)

By Lemma 1 we know that right-hand side of (4c) is 0 iff
Z is Markov. Hence, one can view the KLDR D̄(Z||Z̃) as a
measure of how close a process Z is to a Markov chain.

3

III. M ARKOV C HAIN AGGREGATION
Given a Markov chain X, Markov aggregation deals with
the problem of finding a Markov chain Ỹ on a given smaller
alphabet Y which is the optimal representation of X in the
¯
sense of minimizing a given cost function C(X,
Ỹ). This is
depicted in Fig. 1 by the diagonal arrow and is summarized
in the following definition:
Definition 1 (Markov Aggregation Problem). Let X ∼
¯ ·) be given.
Mar(X , P), Y, and an arbitrary cost function C(·,
The Markov aggregation problem concerns finding a minimizer of
¯
Ỹ)
(5)
min C(X,
Ỹ

where the optimization is over Markov chains on Y.
In this work we address the Markov aggregation problem
using the two-step approach depicted in Fig 1. The first step is
to use a (possibly stochastic) mapping from X to Y. Applying
this mapping to X leads to a stationary process Y which may
not be Markov (in fact, Y is a hidden Markov process). In
the second step we look for the optimal approximation Ỹ of
Y in the sense of Lemma 2.
This two-step approach is a popular method of Markov
aggregation and has been employed in various works including
[5], [7], [8]. In these references, the mapping in the first
step was restricted to be deterministic whereas in this work
we allow it to be stochastic. In other words, while these
references were looking for a partition of X induced by
a function g: X → Y, in this work we permit stochastic
mappings induced by a conditional distribution pY |X . We
represent pY |X as a row stochastic matrix W = [Wx→y ],
where Wx→y = pY |X (y|x).
With this notation, the following corollary to Lemma 2
solves the second of the two steps in our approach, i.e.,
it characterizes the optimal approximation Ỹ of the hidden
Markov process Y that we obtain by observing X through
W:
Corollary 1. Let X ∼ Mar(X , P) and let W denote a
conditional distribution from X to Y. Let Y be the hidden
Markov process obtained by observing X through W, and let
Ỹ ∼ Mar(Y, Q) be its best Markov approximation in the
sense of minimizing D̄(Y||Ỹ) (cf. Lemma 2). Then,
Q = UPW

(6)

where U := diag(ν)−1 WT diag(µ) with ν T := µT W being
the marginal distribution of Yk .
Note that this corollary extends [7, Lem. 3] from deterministic to stochastic mappings.
With the second step solved, the two-step approach to the
optimization problem stated in Definition 1 boils down to
optimization over the mapping W. We can thus restate the
Markov aggregation problem as follows:
Definition 2 (Markov Aggregation Problem Restated). Let
¯ ·) be
X ∼ Mar(X , P), Y, and an arbitrary cost function C(·,
given. Let
¯
C(X, W) = C(X,
Ỹ)
(7)

where Ỹ is the Markov approximation of the hidden Markov
process Y that is obtained by observing X through the
stochastic mapping W. The Markov aggregation problem
using the two-step approach concerns finding a minimizer of
min C(X, W)
W

(8)

where the optimization is over stochastic mappings from X to
Y. If the optimization is restricted over deterministic mappings
g, we abuse notation and write C(X, g) for the cost.
Note that Definition 1 and Definition 2 are not equivalent
in general, i.e., the optimal aggregated chain Ỹ obtained
by solving (5) is not the same as the optimal aggregated
chain Ỹ obtained by solving (8). The two formulations only
become equivalent when we restrict the optimization in (5) to
aggregated Markov chains which can be obtained as a result
of the aforementioned two-step approach.
A. Markov Aggregation via Lumpability
¯
The optimal mapping W depends on the cost function C.
One possible choice in the light of Lemma 2 is
¯
C(X,
Ỹ) = D̄(Y||Ỹ).

(9)

In other words, we wish to find a mapping W such that the
hidden Markov process Y is as close to a Markov chain
as possible in an information-theoretic sense. This may be
reasonable since it states that data obtained by simulating the
aggregated model Ỹ differs not too much from data obtained
by simulating the original model in conjunction with the
stochastic mapping, i.e., data obtained from Y.
There are two shortcomings of the cost (9). First, (9) focuses
only on getting Y close to Ỹ but not on preserving any form
of information in X. This gives rise to trivial optimal solutions:
If W is such that the conditional distribution does not depend
on the conditioning event, i.e., pY |X = pY , or W = 1αT
for some probability vector α, then Y is independent and
identically distributed (i.i.d.) and hence Markov. Indeed, in this
case H̄(Y) = H(Y2 |Y1 ) = H(Y ), from which D̄(Y||Ỹ) = 0
follows. The cost function is thus inappropriate for Markov
aggregation, unless it is regularized appropriately. The second
shortcoming is linked to the fact that D̄(Y||Ỹ) requires,
by (4c), the computation of the entropy rate H̄(Y) of a hidden
Markov process. This problem is inherently difficult [13],
and analytic expressions do not even exist for simple cases
(cf. [14]). In the following, we discuss two previously proposed relaxations of the Markov aggregation problem for
¯
C(X,
Ỹ) = D̄(Y||Ỹ).
The authors of [7] addressed the second shortcoming by
relaxing the cost via
CL (X, W) := H(Y2 |Y1 ) − H(Y2 |X1 ) ≥ D̄(Y||Ỹ).

(10)

This cost does not require computing H̄(Y) and is linked to
the phenomenon of lumpability, the fact that a function of
a Markov chain has the Markov property [12, Thm. 9]: If
CL (X, W) = 0, then Y is a Markov chain.
We now show that, at least for deterministic mappings, this
cost function also has a justification in approximate probabilistic bisimulations, or ε-bisimulations. More specifically, the

4

authors of [15] discussed bisimilarity of Markov processes and
showed that two Markov chains are bisimilar if one can be
described as a function of the other (see discussion after [15,
Def. 5.2]). In other words, if X ∼ Mar(X , P) is a Markov
chain, g: X → Y a surjective function, and Y ∼ Mar(Y, Q)
satisfies Yk = g(Xk ), then X and Y are bisimilar. Since
this is equivalent to lumpability, bisimilarity is implied by
CL (X, g) = 0.
Extending this line of reasoning, we give a justification of
the cost function CL (X, g) in terms of ε-bisimulation of pairs
of Markov chains, even in case X is not lumpable w.r.t. g. To
this end, we adapt [16, Def. 4 & 5] for our purposes:
Definition 3 (ε-Bisimulation). Consider two finite Markov
chains X ∼ Mar(X , P) and Ỹ ∼ Mar(Y, Q) and assume
w.l.o.g. that X and Y are disjoint. We say that X and Ỹ are
ε-bisimilar if there exists a relation Rε ⊆ (X ∪ Y) × (X ∪ Y)
such that for all x ∈ X and y ∈ Y for which (y, x) ∈ Rε ,
and all T ⊆ X ∪ Y we have
X
X
Px→x′ ≥
(11)
Qy→y′ − ε
x′ ∈Rε (T )∩X

y ′ ∈T ∩Y

where Rε (T ) := {s2 ∈ X ∪ Y: s1 ∈ T, (s1 , s2 ) ∈ Rε }.
The definitions of ε-bisimulations are typically given for
labeled [16, Def. 4 & 5] or controlled [2, Def. 4.4] Markov
processes with general alphabets and thus contain more restrictive conditions than our Definition 3. Our definition is
equivalent if the alphabets are finite and if the set of labels is
empty. We are now ready to state
Proposition 1. Let X ∼ Mar(X , P) and the surjective
function g: X → Y be given. Let Q be as in Corollary 1,
where Wx→y = 1 iff y = g(x). Let Ỹ ∼ Mar(Y, Q). Then,
X and Ỹ are ε-bisimilar with
s
ln(2)CL (X, g)
ε=
.
(12)
2 minx∈X µx
Proof: See Section VIII-A.
Despite this justification, the cost function CL (X, W) is
mainly of theoretical interest. The reason is that the shortcoming of leading to trivial solutions is inherited by CL (X, W),
since for W = 1αT one gets CL (X, W) = 0, regardless of α
and P. Even restricting W to be a deterministic partition, as
considered in [7], does not solve this problem: The combinatorial search over all partitions may have its global optimum at a
partition that makes Y close to an i.i.d. process. Indeed, if the
cardinality of Y is not constrained (or if g is not required to be
surjective), then the constant function g yields CL (X, g) = 0.
B. Markov Aggregation by Predictability
A different approach was taken by the authors of [5] who
proposed the following cost function (again with the focus on
deterministic partitions):
CP (X, W) := I(X1 ; X2 ) − I(Y1 ; Y2 )

(13)

The computation of CP (X, W) is simple as it does not require
computing H̄(Y). Furthermore, CP (X, W) reflects the wish

to preserve the temporal dependence structure of X, i.e., it is
connected to predicting future states of Y based on knowledge
of past states of Y. Since X is not i.i.d., observing Xk reveals
some information about Xk+1 . Minimizing CP (X, W) thus
tries to find a W such that Yk reveals as much information
about Yk+1 as possible, and hence does not lead to the
same trivial solutions as CL (X, W) and D̄(Y||Ỹ): A constant
function g or a soft partition W = 1αT render Y1 and Y2
independent, hence the cost is maximized at CP (X, W) =
I(X1 ; X2 ). Unfortunately, as it was shown in [7, Thm. 1],
we have CP (X, W) ≥ CL (X, W), i.e., (13) does not capture
Markovity of Y as well as the relaxation proposed by [7].
Note that although CP (X, W) does not lead to the same
trivial solutions as CL (X, W) and D̄(Y||Ỹ), it still tries to
preserve only part of the information about temporal dependence in X, i.e., information which is helpful in predicting the
next sample. Such a goal is justified in scenarios in which X
is quasi-static, i.e., runs on different time scales: The process
X moves quickly and randomly within a group of states, but
moves only slowly from one group of states to another.
Since all other information contained in X is not necessarily
preserved by minimizing CP (X, W), this cost function can
also lead to undesired solutions: For example, if X is i.i.d.
and hence does not contain any temporal dependence structure,
then I(X1 ; X2 ) = 0 and CP (X, W) = 0 for every mapping
W.
A third objective for optimization may be worth mentioning.
The information contained in X splits into a part describing its
temporal dependence structure (measured by its redundancy
rate I(X1 ; X2 )) and a part describing its new information
generated in each time step (measured by its entropy rate
H̄(X)). Indeed, we have
H(X) = H̄(X) + I(X1 ; X2 ).

(14)

While in this work we focus on preserving Markovity
via CL (X, W) and the temporal dependence structure via
CP (X, W), the authors of [12] investigated conditions such
that the newly generated information (measured by H̄(X)) is
preserved. Developing a Markov aggregation framework that
trades between three different goals – Markovity, temporal
dependence, generated information – is the object of future
work.
IV. R EGULARIZED M ARKOV C HAIN AGGREGATION
In this section we combine the approaches in [7] and [5]
to obtain a new cost function for Markov aggregation. As
discussed in Section III, the aim of Markov aggregation in
[7] is to get a process Y which is as Markov as possible.
This can be captured well by the cost function D̄(Y||Ỹ)
in the light of Lemma 2. The authors of [5] define the
Markov aggregation problem in terms of finding a mapping
W which preserves the temporal dependence structure in X.
This temporal dependence is captured well by the redundancy
rate of the process, which for a Markov chain X equals
R̄(X) = I(X1 ; X2 ). Preserving this temporal dependence
structure is thus well captured by maximizing the redundancy
rate of Y, i.e., by the following optimization problem:
min R̄(X) − R̄(Y)
W

(15)

5

Hence, to combine both the goal of Markovity and the goal of
preserving temporal information, one can define the following
Markov aggregation problem
min (1 − β)D̄(Y||Ỹ) + β(R̄(X) − R̄(Y))
W |
{z
}

(16)

=:δβ (X,W)

where 0 ≤
Definition 1.
ensures that
negative. We

β ≤ 1. Clearly, for β = 0 we are back at
For a general β, the data processing inequality
R̄(X) ≥ R̄(Y), hence the cost (16) is nonmoreover have

Lemma 3. δβ (X, W) is non-decreasing in β.
Proof: See Section VIII-B.
Although minimizing δβ (X, W) tries to preserve both
Markovity and the temporal information in X, the computation
of (16) requires computing H̄(Y). We thus take the approach
of [7] to relax (16). By [7, Thm. 1], we have
D̄(Y||Ỹ) ≤ CL (X, W).

(17)

By the data processing inequality we also have
R̄(X) − R̄(Y) ≤ CP (X, W).

were proposed to solve the respective Markov aggregation
problem.
1
• For β = 2 , optimizing (20) is equivalent to optimizing
CP (X, W). The authors of [5] proposed this cost function for deterministic aggregations, i.e., they proposed
optimizing CP (X, g). Note that this restriction to deterministic aggregations comes at the loss of optimality:
In [17, Example 3] a reversible, three-state Markov chain
was given for which the optimal aggregation to |Y| = 2
states is stochastic. For the bi-partition problem, i.e., for
|Y| = 2, the authors of [5] propose a relaxation to
a spectral, i.e., eigenvector-based optimization problem,
the solution of which has a computational complexity of
O(|X |3 ). In general, this relaxation leads to a further loss
of optimality, even among the search over all deterministic bi-partitions. For a general Y, they suggest to solve
the problem by repeated bi-partitioning, i.e., splitting sets
of states until the desired cardinality is achieved.
• For β = 1, the problem becomes equivalent to maximizing I(X1 ; Y2 ). This is exactly the information bottleneck
problem [18] for a Lagrangian parameter γ → ∞:

(18)

I(X2 ; Y2 ) − γI(X1 ; Y2 ).

Hence, combining the two we get an upper bound on
δβ (X, W) that does not require computing H̄(Y). Indeed, for
β ∈ [0, 1],
δβ (X, W) ≤ (1 − β)CL (X, W) + βCP (X, W).

(19)

•

Rather than the right-hand side of (19), we propose the
following cost function for Markov aggregation:
Cβ (X, W) := (1 − 2β)CL (X, W) + βCP (X, W).

where again β ∈ [0, 1]. One can justify going from (19) to
(20) by noticing that for every 0 ≤ β ≤ 1 for (19), one can
find a 0 ≤ β ≤ 0.5 for (20) such that the two optimization
problems are equivalent, i.e., they have the same optimizer W.
Furthermore, for β = 1, the cost function in (20) corresponds
to information bottleneck problem, a case that is not covered
by (19). Hence, not only is Cβ a strict generalization of δβ but
also has the information bottleneck problem as an interesting
corner case. In the following we summarize some of the
properties of Cβ .

Proof: See Section VIII-C.
V. R ELATED W ORK : S PECIAL C ASES

OF

Cβ (X, W)

We now show that specific settings of β lead to cost
functions that have been proposed previously in the literature.
We list these approached together with the algorithms that

Algorithmic approaches to solving this optimization problem are introduced in [19]. Note, that in this case the
optimal aggregation will be deterministic [17, Thm. 1].
For β = 0, the authors of [7] relaxed their cost function
C0 (X, g) = CL (X, g) as
C0 (X, g) = H(Y2 |Y1 ) − H(Y2 |X1 ) = I(Y2 ; X1 |Y1 )
≤ I(X2 ; X1 |Y1 ) = I(X2 ; X1 ) − I(X2 ; Y1 )
(22)

(20)

Lemma 4. For Cβ and 0 ≤ β ≤ 1 we have:
1) Cβ (X, W) ≥ 0
2) δ0.5 (X, W) = C0.5 (X, W) = 21 CP (X, W)
3) C1 (X, W) = CIB (X, W) := I(X1 ; X2 |Y2 )
4) For β ≤ 12 , βCP (X, W) ≤ δβ (X, W) ≤ Cβ (X, W)
5) For β ≥ 12 , Cβ (X, W) ≤ δβ (X, W) ≤ βCP (X, W)
6) If X is reversible, then Cβ (X, W) is non-decreasing in
β

(21)

•

and proposed using the agglomerative information bottleneck method [20] with the roles of X1 and X2 in (21) exchanged to solve this relaxed optimization problem. The
method has a computational complexity of O(|X |4 ) [19,
Sec. 3.4]. While the mapping minimizing CL (X, W)
may be stochastic, the mapping minimizing (22) will
be deterministic; hence, with this relaxation in mind,
the restriction to deterministic aggregations made in [7]
comes without an additional loss of optimality compared
to what is lost in the relaxation.
The authors of [6] proposed minimizing
I(X1 ; X2 ) − I(X2 ; Y1 ) − γH(Y2 |X1 ).

(23)

They suggested using a deterministic annealing approach,
reducing γ successively until γ = 0. In the limiting case,
the cost function then coincides with (22) and the optimal
aggregation is again deterministic. Note that, for reversible Markov chains, we have I(X1 ; Y2 ) = I(X2 ; Y1 ),
hence both (22) and (23) (for γ = 0) are equivalent to
C1 (X, W). Analyzing [6, Sec. III.B] shows that in each
annealing step the quantity
D(pX2 |X1 =x ||pX2 |Y1 =y )
X
pX |X (x′ |x)
:=
pX2 |X1 (x′ |x) log 2 1 ′
pX2 |Y1 (x |y)
′
x ∈X

(24)

6

has to be computed for every x and y. Hence, the computational complexity of this approach is O(|Y| · |X |2 )
in each annealing step.
VI. M ARKOV C HAIN AGGREGATION A LGORITHMS
We now propose an iterative method for optimizing (20)
over deterministic aggregations for general values of β. The
method consists of a sequential optimization algorithm (Algorithm 1) and an annealing procedure for β (Algorithm 2)
that prevents getting stuck in local optima. Since we focus
only on deterministic aggregations, in the remainder of this
section we can replace Cβ (X, W) by Cβ (X, g) for some
g: X → Y. Our algorithm has a computational complexity of
O(|Y| · |X |2 ) per iteration. Note, however, that the restriction
to deterministic aggregation functions comes, at least for some
values of β, with a loss of optimality, i.e., in general we have
minW Cβ (X, W) ≤ ming Cβ (X, g).
A. Sequential Algorithm
We briefly illustrate an iteration of Algorithm 1: Suppose
x ∈ X is mapped to the aggregate state y ∈ Y, i.e., g(x) = y.
We remove x from aggregate state y. We then assign x to
every aggregate state y ′ , y ′ ∈ Y, while keeping the rest of the
mapping g the same and evaluate the cost function. Finally, we
assign x to the aggregate state that minimized the cost function
(breaking ties, if necessary). This procedure is repeated for
every x ∈ X .
Algorithm 1 Sequential Generalized Information-Theoretic
Markov Aggregation.
1: function g = S GITMA(P, β, |Y|, #itermax , optional:
initial aggregation function ginit )
2:
if ginit is empty then
⊲ Initialization
3:
g ← Random Aggregation Function
4:
else
5:
g ← ginit
6:
end if
7:
#iter ← 0
8:
while #iter < #itermax do
⊲ Main Loop
9:
for all elements x ∈ X do
⊲ Optimizing g
10:
for all aggregate
( states y ∈ Y do
g(x′ )
x′ 6= x
11:
gy (x′ ) =
⊲ Assign x
y
x′ = x
to aggregate state y
12:
Cgy = Cβ (X, gy )
13:
end for
⊲ (break ties)
14:
g = arg min Cgy
gy

15:
16:
17:
18:

end for
#iter ← #iter + 1
end while
end function

It is easy to verify that the cost function is reduced in each
step of Algorithm 1, as a state is only assigned to a different
aggregate state if the cost function is reduced. Hence, the
algorithm modifies g in each iteration in order to reduce the

cost until it either reaches the maximum number of iterations
or until the cost converges.
Note that the algorithm is random in the sense that it
is started with a random aggregation function g. Depending
on the specific application, though, a tailored initialization
procedure may lead to performance improvements.
Finally, it is worth mentioning that for β = 1 our Algorithm 1 is equivalent to the sequential information bottleneck
algorithm proposed in [19, Sec. 3.4].
B. Annealing Procedure for β
Although Algorithm 1 is guaranteed to converge (with
proper tiebreaking), convergence to a global optimum is not
ensured. The algorithm may get stuck in poor local minima.
This happens particularly often for small values of β, as our
experiments in Section VII-B show. The reason is that, for
small β, Cβ (X, W) has many poor local minima and, randomly
initialized, the algorithm is more likely to get stuck in one of
them. In contrast, our results suggest that for larger values of
β the cost function has only few poor local minima and that
the algorithm converges to a good local or a global minimum
for a significant portion of random initializations.
A solution for small β would thus be to choose an initialization that is close to a “good” local optimum. A simple idea is
thus to re-use the function g obtained for a large value of β as
initial aggregation for smaller values of β. We thus propose the
following annealing algorithm: We initialize β = 1 to obtain
g. Then, in each iteration of the annealing procedure, β is
reduced and the aggregation function is updated, starting from
the result of the previous iteration. The procedure stops when β
reaches the desired value, βtarget . The β-annealing algorithm is
sketched as Algorithm 2. As is clear from the description, the
β-annealing algorithm closely follows graduated optimization
in spirit [21]. The results for synthetic datasets with and
without β-annealing are discussed in Section VII-B, which
show that without restarts one keeps getting stuck in bad local
optima for small β, while with β-annealing one is able to avoid
them. Furthermore in our experiments we have observed that
β-annealing achieves good results for random initializations,
hence tailoring initialization procedures is not necessary at
least for the scenarios we considered.
Algorithm 2 β-Annealing Information-Theoretic Markov Aggregation
1: function g = A NN ITMA(P, βtarget , |Y|, #itermax , ∆)
2:
β←1
3:
g = sGITMA(P, β, |Y|, #itermax )
⊲ Inizialization
4:
while β > βtarget do
5:
β ← max{β − ∆, βtarget }
6:
g = sGITMA(P, β, |Y|, #itermax , g)
7:
end while
8: end function
Note that the β-annealing algorithm admits producing results for a series of values of β at once: Keeping all intermediate aggregation functions, one obtains aggregations
for all values of β in the set {1, 1 − ∆, 1 − 2∆, . . . , 1 −

7

1−β

∆⌈ ∆target ⌉, βtarget }. The aggregations one obtains are exactly
those one would obtain from restarting A NN ITMA for each
value in this set, each time with the same random initial
partition. We used this fact in our experiments: If we were
interested in results for βtarget ranging between 0 and 1 in
steps of 0.05, rather than restarting A NN ITMA for each value
in this set, we started A NN ITMA for βtarget = 0 and ∆ = 0.05
once, keeping all intermediate results.
C. Computational Complexity of the Sequential Algorithm
Note that the asymptotic computational complexity of Algorithm 2 equals that of Algorithm 1, since the former simply
calls the latter ⌈(1 − βtarget )/∆⌉ + 1 times. We thus only
evaluate the complexity of Algorithm 1. To this end, we first
observe that the cost Cβ (X, W) can be expressed with only
three mutual information terms for any β and W:
Cβ (X, g) = βI(X1 ; X2 ) + (1 − 2β)I(X1 ; g(X2 ))
− (1 − β)I(g(X1 ); g(X2 )). (25)
The first term I(X1 ; X2 ) is constant regardless of the aggregation hence the computation of Cβ (X, g) depends upon the
computation of the other two terms.
In each iteration of the main loop, we evaluate Cβ (X, gy )
in line 12 for each x ∈ X and y ∈ Y. Note that gy differs
from the current g only for one element as defined in line
11. Thus, the joint PMF pX1 ,gy (X2 ) differs from pX1 ,g(X2 ) in
only two rows and hence can be computed from pX1 ,g(X2 )
in O(|X |) computations. Moreover, I(X1 ; gy (X2 )) can be
computed from I(X1 ; g(X2 )) in O(|X |) computations, cf. [20,
Prop. 1]. This is due to the fact that we can write [22,
eq. (2.28)]
I(X1 ; gy (X2 )) = I(X1 ; g(X2 ))


X
pX1 ,gy (X2 ) (x1 , y2 )
+
pX1 ,gy (X2 ) (x1 , y2 ) log
pX1 (x1 )pgy (X2 ) (y2 )
x1 ∈X
y2 ∈{y,g(x)}

−

X

pX1 ,g(X2 ) (x1 , y2 ) log

x1 ∈X
y2 ∈{y,g(x)}



pX1 ,g(X2 ) (x1 , y2 )
pX1 (x1 )pg(X2 ) (y2 )



.

(26)
The term I(gy (X1 ); gy (X2 )) can be computed from
I(g(X1 ); g(X2 )) in O(|Y|) computations, but requires the
updated joint PMF pgy (X1 ),gy (X2 ) . This PMF can be computed
from pg(X1 ),g(X2 ) in O(|X |) computations. Combining this
with the fact that line 12 is executed once for each aggregate
state in Y and once for each state in X in every iteration,
we get that optimizing g has a computational complexity of
O(|Y| · |X |2 ) per iteration.
VII. E XPERIMENTS

AND

R ESULTS

A. A Non-Reversible Markov Chain
The last property of Lemma 4 cannot be generalized to
non-reversible Markov chains. Specifically, as the proof of
Lemma 4 shows, Cβ is non-decreasing in β iff CP ≥ 2CL .
Since one can find also non-reversible Markov chains for

which this holds, reversibility is sufficient but not necessary
for Cβ to be non-decreasing in β. We next consider a nonreversible Markov chain X ∼ Mar({1, 2, 3}, P) with


0.4
0.3
0.3
P =  0.25 0.3
0.45 
(27)
0.15 0.425 0.425

and let g be such that g(1) = 1 and g(2) = g(3) = 2. Then,
CL = 0.0086 and CP = 0.0135, i.e., CP < 2CL . In this case,
Cβ is decreasing with increasing β.
B. Quasi-Lumpable and Nearly Completely Decomposable
Markov Chains
Suppose we have a partition {Xi }, i = 1, . . . , M , of X
with |Xi | = Ni . Then for any A′ and P′ij that are M × M
and Ni × Nj row stochastic matrices, respectively, define A =
[aij ] = (1 − α)A′ + αI, α ∈ [0, 1], and let


a11 P′11
a12 P′12 · · · a1M P′1M
 a21 P′21
a22 P′22 · · · a2M P′2M 


P′ = 
 . (28)
..
..
..
..


.
.
.
.
′
′
′
aM1 PM1 aM2 P22 · · · aMM PMM

Let further X ∼ Mar(X , P′ ). If g induces the partition
{Xi }, then it can be shown that Y is Markov with transition
probability matrix A, i.e., Ỹ ≡ Y, and C0 (X, Ỹ) = 0. The
Markov chain X is lumpable w.r.t. the partition g. The matrix
P′ is block stochastic and the parameter α specifies how
dominant the diagonal blocks are. Specifically, if α = 1, then
P′ is block diagonal and we call X completely decomposable.
Such a Markov chain is not irreducible. We hence look at
Markov chains X ∼ Mar(X , P) with
P = (1 − ε)P′ + εE

(29)

where ε ∈ [0, 1] and where E (which can be interpreted as
noise) is row stochastic and irreducible. For small values of
ε we call X nearly completely decomposable (NCD) if α is
close to one, otherwise we call it quasi-lumpable.
We now perform experiments with these types of Markov
chains. We set M = 3, N1 = N2 = 25, and N3 = 50,
and chose the parameters from α ∈ {0, 0.5, 0.95} and ε ∈
{0, 0.4, 0.8}. For each pair (α, ε), we generated 250 random
matrices A′ and P′ij . A selection of the corresponding matrices
P is shown in Fig. 2.
Note that in practice the states of even a completely decomposable Markov chain X are rarely ordered such that the
transition probability matrix is block diagonal. Rather, the state
labeling must be assumed to be random. In this case, P is
obtained by a random permutation of the rows and columns
of a block diagonal matrix (see Fig. 2(d)), which prevents
the optimal aggregation function being “read off” simply by
looking at P. That P has a block structure in our case does
not affect the performance of our algorithms, since they 1) are
unaware of this structure and 2) are initialized randomly.
We applied our aggregation algorithm both with and without
the annealing procedure for β ∈ {0, 0.1, . . . , 0.9, 1} with the
goal of retrieving the partition {Xi }. We measure the success,
i.e., the degree to which the function g obtained from the

8

10

20

20

20

40

40

40

60

60

60

80

80

80

100

100

100

20
30
40
50
60
70
80
90

20

40

60

80

100

20

(a) α = ε = 0

40

60

80

100

20

(b) α = 0.95, ε = 0.8

40

60

80

100

(c) α = 0.95, ε = 0.4

100
20

40

60

80

100

(d) α = 0.95, ε = 0.4; rows and
columns are permuted

0.2

ARI

Cβ

Cβ

1

0.2

0.15
0.1

0.1

0.5

5 · 10−2
0

0

0.2

0.4

0.6

0.8

0

1

0
0

0.2

0.4

β

0.8

0

1

α=0 (lumpable)
α=0.5
α=0.95 (NCD)
0.2

0.4

0.6

0.8

0.5

0
1

1

0.5

0
0

0.2

0.4

β
(h) ARI with β-annealing, ε = 0

0.8

1
ARI

ARI

0.5

0.6

(g) ARI without β-annealing, ε = 0.4

1

0

0.4
β

(f) Cost without β-annealing, ε = 0.4

1

0

0.2

β

(e) Cost with β-annealing, ε = 0.4

ARI

0.6

0.6

0.8

β
(i) ARI with β-annealing, ε = 0.4

1

0

0.2

0.4

0.6

0.8

1

β
(j) ARI with β-annealing, ε = 0.8

Fig. 2. (a)-(c): Colorplots of the transition probability matrices P for different values of α and ε. For large α the block diagonal structure becomes more
dominant. (d): A random permutation of the rows and columns hides the block structure. (e)-(j): Curves showing the cost function Cβ and the adjusted Rand
index (ARI) for different settings. Mean values (solid lines) are shown together with the standard deviation (shaded areas).

algorithm agrees with the partition {Xi }, using the adjusted
Rand index (ARI). An ARI of one indicates that the two
partitions are equivalent. Note that we always assume that the
number M of sets in the partition {Xi } is known.
The results are shown in Fig. 2. Specifically, Fig. 2(e) shows
that the cost for the aggregation found by our algorithm with
β-annealing decreases monotonically with decreasing β: We
obtain a partition for a given value of β. This partition has,
by assuming CP ≥ 2CL (cf. Section VII-A), an even lower
cost for a smaller value of β. Further optimization for this
smaller value of β further reduces the cost, leading to the
depicted phenomenon. In contrast, the sequential Algorithm 1
without the annealing procedure fails for values of β less than
0.5. This is apparent both in the cost in Fig. 2(f) (which has

a sharp jump around β = 0.5) and in the ARI in Fig. 2(g)
(which drops to zero). Apparently, the algorithm gets stuck in
a bad local optimum.
Figs. 2(h) to 2(j) show the ARI of the aggregations obtained
by our algorithm with β-annealing. First of all, it can be
seen that performance improves with increasing α, since
the dominant block structure makes discovering the correct
partition more easy. Moreover, it can be seen that for α = 0 the
optimum β lies at smaller values, typically smaller than 0.5.
The position of this optimum increases with increasing noise:
While in the noiseless case the correct partition is typically
obtained for β close to zero, in the highly noisy case of
ε = 0.8 we require β ≈ 0.4 to achieve good results. The
reason may be that the higher noise leads to more partitions

9

A GGREGATING A

β Value
Ref.
β=1
β = 0.8
β = 0.5
β=0

ARI
–
0.43
0.46
0.35
0.12

A GGREGATING A

β Value
Ref.
β=1
β = 0.8
β = 0.5
β=0

ARI
–
0.2
0.34
0.24
0.35
0.15
0.31
0.01
0.02

LETTER BI - GRAM MODEL .

TABLE I
T HE PARTITIONS ARE SHOWN TOGETHER WITH THE ARI ARI W. R . T. THE REFERENCE PARTITION ( FIRST
ROW ) FOR |Y| = 4

Partitions, shown for |Y| = 4
{ },{!"$’(),-.:;?[]},{aeiou},{0123456789},{AEIOU},{BCDFGHJKLMNPQRSTVWYZ},{bcdfghjklmnpqrstvwxyz}
{ !’),-.0:;?]},{aeioy},{"$(123456789ABCDEFGHIJKLMNOPQRSTUVWY[h},{Zbcdfgjklmnpqrstuvwxz}
{ !’),-.:;?Z]},{aeiouy},{"$(0123456789ABCDEFGHIJKLMNOPQRSTUVWY[h},{bcdfgjklmnpqrstvwxz}
{ !3?Z},{’2456789AOUaeiou},{"$(-01BCDEFGHIJKLMNPQRSTVWY[bhjqw},{),.:;]cdfgklmnprstvxyz}
{ -2CEFMPSTcfgopst},{"’456789AOUZaeiu},{!$1?BDGHJLNQRVW[bhjklmqrvwz},{(),.03:;IKY]dnxy}

LETTER BI - GRAM MODEL .

TABLE II
T HE PARTITIONS ARE SHOWN TOGETHER WITH THE ARI ARI W. R . T. THE REFERENCE PARTITION ( FIRST
ROW ) FOR |Y| ∈ {2, 7}

Partitions, shown for |Y| ∈ {2, 7}
{ },{!"$’(),-.:;?[]},{aeiou},{0123456789},{AEIOU},{BCDFGHJKLMNPQRSTVWYZ},{bcdfghjklmnpqrstvwxyz}
{ !"’),-.01235689:;?KU]aehioy},{$(47ABCDEFGHIJLMNOPQRSTVWYZ[bcdfgjklmnpqrstuvwxz}
{ !’),-.:;?]},{aeioy},{"$(0123456789ABCDEFGHIJLMNOPQRSTVWY[},{bcfjmpqstw},{dgx},{KUh},{Zklnruvz}
{ !"’),-.01235689:;?EU]aehiouy},{$(47ABCDFGHIJKLMNOPQRSTVWYZ[bcdfgjklmnpqrstvwxz}
{ !’),-.:;?]},{aeioy},{$"(0123456789ABCDEFGHIJLMNOPQRSTUVWYZ[},{Kh},{bcfjkmpqstw},{dg},{lnruvxz}
{ !’-12368?EOUZaeiou},{"$(),.04579:;ABCDFGHIJKLMNPQRSTVWY[]bcdfghjklmnpqrstvwxyz}
{ ’},{!),.:;?]dy},{aeiou},{"$(-0123589ACEIMOPRSTUWZ},{BDFGHJKLNQVYhj},{7[bcfgkmpqstw},{46lnrvxz}
{ !$(-0124578?ABCFHLMNOPRSTUVWaceglnostuwxz},{"’),.369:;DEGIJKQYZ[]bdfhijkmpqrvy}
{ 4689ao},{$’AKOiux},{!?HVZhjkmvz},{"(-25CEFLMNRUWY[egnprs},{37BPQbl},{1:;STctw},{),.0DGIJ]dfy}

being quasi-lumpable by leading to an i.i.d. Y, hence for small
values of β one may get drawn into these “false solutions”
more easily. In contrast, for NCD Markov chains (i.e., for
α = 0.95) sometimes noise helps in discovering the correct
partition. Comparing Figs. 2(h) and 2(i), one can see that a
noise of ε = 0.4 allows us to perfectly discover the partition.
We believe that a small amount of noise helps in escaping bad
local minima.
The fact that the β for which the highest ARI is achieved not
necessarily falls together with the values 0, 0.5, or 1 indicates
that our generalized aggregation framework has the potential to
strictly outperform aggregation cost functions and algorithms
that have been previously proposed (cf. Section V).
C. An Example from Natural Language Processing
We took the letter bi-gram model from [23], which was
obtained by analyzing the co-occurrence of letters in F. Scott
Fitzgerald’s book “The Great Gatsby”. The text was modified
by removing chapter headings, line breaks, underscores, and
by replacing é by e. With the remaining symbols, we obtained
a Markov chain with an alphabet size of N = 76 (upper and
lower case letters, numbers, punctuation, etc.).
We applied Algorithm 2 for |Y| ∈ {2, . . . , 7} and β ∈
{0, 0.1, . . . , 0.9, 1}. To get consistent results, we restarted the
algorithm 20 times for β = 1 and chose the aggregation g
that minimized C1 (X, g); we used this aggregation g as an
initialization for the β-annealing procedure.
Looking at the results for |Y| = 4 in Table I, one can
observe that the results for β = 0.8 appear to be most
meaningful when compared to other values of β such as
β = 1 (information bottleneck), β = 0.5 (as proposed
in [5]), and β = 0 (as proposed in [7]). Specifically, for
β = 0 it can be seen that not even the annealing procedure
was able to achieve meaningful results. This conclusion is
supported by calculating the ARI of these aggregations for
a plausible reference aggregation of the alphabet into upper

case vowels, upper case consonants, lower case vowels, lower
case consonants, numbers, punctuation, and the blank space as
shown in the first row of the Table I. The absolute ARI values
are not a good performance indicator in this case since we are
comparing to a reference partition with seven sets whereas
|Y| = 4.
In Table II the same experiment is repeated for |Y| ∈ {2, 7}.
We again observe that β = 0.8 leads to the most meaningful
results which is also supported by ARI values.
D. Clustering via Markov Aggregation
Data points are often described only by pairwise similarity
values, and these similarity values can be used to construct the
transition probability matrix of a Markov chain. Then, with
this probabilistic interpretation, our information-theoretic cost
functions for Markov aggregation can be used for clustering.
This approach has been taken by [3], [4].
We considered two different data sets: three linearly separable clusters and three concentric circles, as shown in Fig 3. The
three linearly separable clusters were obtained by placing 40,
20, and 40 points, drawn from circularly symmetric Gaussian
distributions with standard deviations 2.5, 0.5, and 1.5 at
horizontal coordinates -10, 0, and 10, respectively. The three
concentric circles were obtained by placing 40 points each at
uniformly random angles at radii {0.1, 7, 15}, and by adding
to each data point spherical Gaussian noise with a standard
deviation of 0.3. In both cases, we computed the transition
probability matrix P according to
Pi→j ∝ e

−

kxi −xj k2
2
σk

(30)

where xi and xj are the coordinates of the i-th and j-th data
point, k ·k22 is the squared Euclidean distance, and where σk is
a scale parameter. We set σk to the average squared Euclidean
distance between a data point and its k nearest neighbors (and

10

20

b
b

b

b

80

40

60

80

100

b b
b

bb

b
b

b

b
bb

bb
b

b

b

bb
b

b

20

30

30

40

40

50

50

60

60

70

70

80

80

b

b
b bb b

b

b b
b

bb

b

bb
b
bb

bb

40

60

80

100

(c) P, k = 120

b
b

b

b

b

b

b

b

b
b

b

b

b bb
b

b bb
b

b

b bb b
b

b

b

b

b bb b

b b

b
b

b

b

b

b

b

bb
b b b b bb bb b bb b b bb b b
b b

b bb

b

bb

b

bb
b
b

b
bb

b
b

bb

bb b
b
b

bbb
b

b

b

bb
b

b

b
b

b

b b
b

b

bb b

{0.2, 0.5, 0.8},

10

90

b

120

20

b

b

bb
b b b b bb bb b bb b b bb b b
b b

b

10

b bb
b

b bb
b

100

bbb

b

b bb

b

b

80

b bb b
b

b
60

b

20

(b) β ∈
k = 15

b

b

120

(a) P, k = 15

b
b

bb
b

b

b

b

40

b

b
b

b

b b

b

b
b bb b

b
20

b

b

b

bb
b b b b bb bb b b b bb bb b b
b b

b
bb b

120

b bb
b

b bb
b

100

b

b

20

b

b bb
b

b

b

b bb b
b

40

60

b
b

b

b b

bb

b

b

b

b

b

bbb
b

120

(d) β = 0.2, k = 120

b
b

b
b
b
b
b

b

b
b

bb

b
b

b
b

b

b
b
b

b

b
b

b
b
b

b
b

b
b

b
b

b

(e) β = 0.8, k = 120

b
b
b
b
b

b b b b b b bb
b bb b b b
b b
b

b

b
b

b b
b

b

b

b
bb

b

b

b
b
b b
b
b b b
b b bb bb
b
b
b bb
b

b

b b

b
b

90

100

100
20

40

60

80

(f) P, k = 15

100

20

40

60

80

100

(h) β ∈ {0.2, 0.5, 0.8}, k ∈ {15, 100}

(g) P, k = 100

Fig. 3. Clustering three circles (first row) and three linearly separable clusters (second row). For k = 15, the transition probability matrices (shown in (a)
and (f)) are nearly completely decomposable. The result for the three circles depends strongly on a careful setting of the parameters β and k ((b), (d), and
(e)), while the three linearly separable clusters were separated correctly for all parameter choices (h).

averaged this quantity over all data points). We set k either to
15 or to the total number of data points.
We applied our Algorithm 2 with the annealing procedure
for β. As in the previous experiment, we restarted the algorithm 50 times for β = 1 and chose the aggregation g
that minimized C1 (X, g); we used this aggregation g as an
initialization for the β-annealing procedure.
The results are shown in Fig. 3, together with a colorplot
of the respective transition probability matrices. It can be
seen that the three linearly separable clusters were detected
correctly for all chosen parameter values. This is not surprising
for k = 15, since in this case the resulting Markov chain
is nearly completely decomposable. Interestingly, though, the
same results were observed for k = 100 for which P is
structured, but not block diagonal. One may claim that these
results are due to Algorithm 2 getting stuck in a local optimum
for β = 1 which accidentally coincides with the correct
clustering, and that optimizing our cost function for values
of β larger than 0.5 but smaller than 1 may fail. Since we
reproduced these results by using Algorithm 1 (with 50 restarts
to escape bad local optima) for values of β greater than 0.5,
this claim can be refuted.
For the three concentric circles, things look different. We
correctly identified the clusters only for a nearly completely
decomposable P, i.e., for a careful setting of k (and we were
able to reproduce these results for β greater than 0.5 using
Algorithm 1). For k = 120, i.e., equal to the number of data
points, the three circles were not identified correctly.
Since we have reason to believe that the optimal k depends
strongly on the data set, we are hesitant to recommend this approach to cluster data points that are not linearly separable (in
which case a simpler method such as k-means would suffice).
Our preliminary analysis of [3] suggests that their approach

(in which X is a random walk on the k-nearest neighbor graph
of the data set and in which the authors chose β = 0.5) suffers
from similar problems. Finally, the authors of [4] suggest to
let X “relax” to some metastable point, i.e., take an r-th
power of P such that Pr is approximately a projection; their
approach is equivalent to ours for β = 1, with P replaced
by Pr . Nevertheless, also this approach requires setting r and
k for (30). Whether this relaxation to metastability can be
successfully combined with our generalized cost function for
Markov aggregation will be deferred to future investigations.
VIII. P ROOFS
A. Proof of Proposition 1
Consider the relation Rε = {(g(x), x): x ∈ X }. It can be
shown that
∀T ⊆ X ∪ Y: Rε (T ) = g −1 (T ∩ Y) ⊆ X .

(31)

We thus need to show that, for all x and all B ⊆ Y,
X
X
Px→x′ ≥
Qg(x)→y − ε.

(32)

Now let R = [Rx→y ] = PW, i.e., we have
X
Rx→y =
Px→x′ .

(33)

x′ ∈g−1 (B)

y∈B

x′ ∈g−1 (y)

One can show along the lines of [7, (65)–(68)] that
X
X
Rx→y
CL (X, W) =
µx
Rx→y log
Qg(x)→y
x∈X
y∈Y
{z
}
|

(34)

=:D(Rx→· ||Qg(x)→· )

from which we get that, for every x,
D(Rx→· ||Qg(x)→· ) ≤

CL (X, W)
.
minx µx

(35)

11

With Pinsker’s inequality [22, Lemma 12.6.1] and [22,
(12.137)] we thus get that, for every x and every B ⊆ Y,
s
X
ln(2)CL (X, W)
.
(36)
Rx→y − Qg(x)→y ≤
2 minx µx
y∈B

Combining this with (33) thus shows that (32) holds for
s
ln(2)CL (X, W)
.
(37)
ε=
2 minx µx

δβ (X, W) − βI(X1 ; X2 )
= (1 − β)H(Y2 |Y1 ) − (1 − 2β)H̄(Y) − βH(Y )
≤ (1 − β)H(Y2 |Y1 ) − (1 − 2β)H(Y2 |X1 ) − βH(Y )
= (1 − 2β)H(Y2 |Y1 ) − (1 − 2β)H(Y2 |X1 ) − βI(Y1 ; Y2 )
= (1 − 2β)CL (X, W) − βI(Y1 ; Y2 ).
The inequality is reversed for β ≥ 0.5.
For the fifth property, we repeat the last steps with
− (1 − 2β)H̄(Y) ≤ −(1 − 2β)H(Y2 |Y1 )

This completes the proof.
B. Proof of Lemma 3
We show that the derivative of δβ (X, W) w.r.t. β is positive.
Indeed,
d
δβ (X, W)
dβ
= R̄(X) − R̄(Y) − D̄(Y||Ỹ)
(38)
= I(X1 ; X2 ) − H(Y ) + H̄(Y) − H(Y2 |Y1 ) + H̄(Y). (39)
The entropy rate of the reversed process equals the entropy
rate of the original process, i.e.,
H̄(Y) = lim H(Yn |Y1n−1 ) = lim H(Y1 |Y2n ).
n→∞

The fourth property is obtained by observing that, if β ≤ 0.5

n→∞

(40)

We can now apply [22, Lem. 4.4.1] to both sides to get
H̄(Y) ≥ H(Y2 |X1 ) and H̄(Y) ≥ H(Y1 |X2 ). We use this
in the derivative to get
d
δβ (X, W)
dβ
≥ I(X1 ; X2 ) − H(Y ) + H(Y1 |X2 ) − H(Y2 |Y1 ) + H(Y2 |X1 )
(41)
= H(X|Y ) − H(X1 |Y1 , X2 ) − H(Y2 |Y1 ) + H(Y2 |X1 )
(42)
= I(X1 ; X2 |Y1 ) − I(X1 ; Y2 |Y1 ) ≥ 0

(43)

by data processing.
C. Proof of Lemma 4
The first property follows by recognizing that
Cβ (X, W)
= (1 − β)CL (X, W) + β(CP (X, W) − CL (X, W))

(44)

and that CP (X, W) ≥ CL (X, W).
The second property follows immediately from the definition of δβ (X, W) and CP (X, W).
For the third property, note that
C1 (X, W) = CP (X, W) − CL (X, W)
= I(X1 ; X2 ) − H(Y ) + H(Y2 |X1 )
= I(X1 ; X2 ) − I(X1 , Y2 ) = I(X1 ; X2 |Y2 ).

(45)

noticing that (1 − 2β) ≤ 0 if β ≥ 0.5. Again, the inequality
is reversed for β ≤ 0.5.
If X is reversible, then the PMFs do not change if the
order of the indices is reversed. As a consequence, we
have I(X1 ; X2 |Y2 ) = I(X2 ; X1 |Y1 ) = C1 (X, W). But
CL (X; W) = I(Y2 ; X1 |Y1 ) ≤ C1 (X, W) by data processing.
Thus, the sixth property follows by noting that, with (44),
Cβ (X, W) = (1 − β)CL (X, W) + βC1 (X, W).
ACKNOWLEDGMENTS
The authors thank Matthias Rungger and Majid Zamani,
both from Hybrid Control Systems Group, Technical University of Munich, for discussions suggesting the connection
between lumpability and bisimulation. The work of Rana Ali
Amjad was supported by the German Ministry of Education
and Research in the framework of an Alexander von Humboldt
Professorship. The work of Bernhard C. Geiger was funded
by the Erwin Schrödinger Fellowship J 3765 of the Austrian
Science Fund.
R EFERENCES
[1] M. A. Katsoulakis and J. Trashorras, “Information loss in coarsegraining of stochastic particle dynamics,” J. Stat. Phys., vol. 122, no. 1,
pp. 115–135, 2006.
[2] A. Abate, “Approximation metrics based on probabilistic bisimulations
for general state-space Markov processes: A survey,” Electronic Notes
in Theoretical Computer Science, vol. 297, pp. 3 – 25, 2013, Proc.
Workshop on Hybrid Autonomous Systems.
[3] A. Alush, A. Friedman, and J. Goldberger, “Pairwise clustering based on
the mutual-information criterion,” Neurocomputing, vol. 182, pp. 284–
293, 2016.
[4] N. Tishby and N. Slonim, “Data clustering by Markovian relaxation and
the information bottleneck method,” in Advances in Neural Information
Processing Systems (NIPS), Denver, CO, Nov. 2000. [Online]. Available:
http://citeseer.ist.psu.edu/viewdoc/summary?doi=10.1.1.24.3488
[5] K. Deng, P. G. Mehta, and S. P. Meyn, “Optimal Kullback-Leibler
aggregation via spectral theory of Markov chains,” IEEE Trans. Autom.
Control, vol. 56, no. 12, pp. 2793–2808, Dec. 2011.
[6] Y. Xu, S. M. Salapaka, and C. L. Beck, “Aggregation of graph models
and Markov chains by deterministic annealing,” IEEE Trans. Autom.
Control, vol. 59, no. 10, pp. 2807–2812, Oct. 2014.
[7] B. C. Geiger, T. Petrov, G. Kubin, and H. Koeppl, “Optimal KullbackLeibler aggregation via information bottleneck,” IEEE Trans. Autom.
Control, vol. 60, no. 4, pp. 1010–1022, Apr. 2015, open-access:
arXiv:1304.6603 [].
[8] M. Vidyasagar, “Reduced-order modeling of Markov and hidden Markov
processes via aggregation,” in Proc. IEEE Conf. on Decision and Control
(CDC), Atlanta, GA, Dec. 2010, pp. 1810–1815.
[9] J. G. Kemeny and J. L. Snell, Finite Markov Chains, 2nd ed. Springer,
1976.
[10] R. M. Gray, Entropy and Information Theory. New York, NY: Springer,
1990.

12

[11] Z. Rached, F. Alajaji, and L. L. Campbell, “The Kullback-Leibler
divergence rate between Markov sources,” IEEE Trans. Inf. Theory,
vol. 50, no. 5, pp. 917–921, May 2004.
[12] B. C. Geiger and C. Temmel, “Lumpings of Markov chains, entropy rate preservation, and higher-order lumpability,” J. Appl. Probab.,
vol. 51, no. 4, pp. 1114–1132, Dec. 2014, extended version available:
arXiv:1212.4375 [].
[13] D. Blackwell, “The entropy of functions of finite-state {M}arkov
chains,” in Trans. first {P}rague Conf. Inf. theory, {S}tatistical Decis. Funct. random Process. held {L}iblice near {P}rague from
{N}ovember 28 to 30, 1956.
Prague: Publishing House of the
Czechoslovak Academy of Sciences, 1957, pp. 13–20.
[14] O. Ordentlich, “Novel lower bounds on the entropy rate of binary hidden
Markov processes,” in Proc. IEEE Int. Sym. on Information Theory
(ISIT), Jul. 2016, pp. 690–694.
[15] J. Desharnais, A. Edalat, and P. Panangaden, “Bisimulation for labelled
Markov processes,” Information and Computation, vol. 179, no. 2, pp.
163 – 193, 2002.
[16] G. Bian and A. Abate, “On the relationship between bisimulation and
trace equivalence in an approximate probabilistic context,” in Proc. Int.
Conf. on Foundations of Software Science and Computation Structure
(FOSSACS), J. Esparza and A. S. Murawski, Eds. Uppsala: Springer
Berlin Heidelberg, Apr. 2017, pp. 321–337.
[17] B. C. Geiger and R. A. Amjad, “Mutual information-based clustering:
Hard or soft?” in Proc. Int. ITG Conf. on Systems, Communications and Coding (SCC), Hamburg, Feb. 2017, pp. 1–6, open-access:
arXiv:1608.04872 [].
[18] N. Tishby, F. C. Pereira, and W. Bialek, “The information bottleneck
method,” in Proc. Allerton Conf. on Communication, Control, and
Computing, Monticello, IL, Sep. 1999, pp. 368–377.
[19] N. Slonim, “The information bottleneck: Theory and applications,” Ph.D.
dissertation, Hebrew University of Jerusalem, 2002.
[20] N. Slonim and N. Tishby, “Agglomerative information bottleneck,” in
Advances in Neural Information Processing Systems (NIPS), Denver,
CO, Nov. 1999, pp. 617–623.
[21] A. Blake and A. Zisserman, Visual Reconstruction. Cambridge, MA,
USA: MIT Press, 1987.
[22] T. M. Cover and J. A. Thomas, Elements of Information Theory, 1st ed.
Wiley Interscience, 1991.
[23] B. C. Geiger and Y. Wu, “Higher-order optimal Kullback-Leibler
aggregation of Markov chains,” in Proc. Int. ITG Conf. on Systems,
Communications and Coding (SCC), Hamburg, Feb. 2017, pp. 1–6,
open-access: arXiv:1608.04637 [].

