# Markov Chain Aggregation with Error Bounds on Transient Distributions

Fabian Michel  Markus Siegle

(08/02/2024)

X 2pt 2pt

Markov Chain Aggregation with Error Bounds on Transient Distributions
Fabian Michel, Markus Siegle

## Abstract

We extend the existing theory of formal error bounds for the transient
distribution of an aggregated (or lumped) Markov chain when compared to the
transient distribution of the original chain, for both discrete- and
continuous-time Markov chains. In the discrete-time setting, we bound the
stepwise increment of the error, and in the continuous-time setting, we bound
the rate at which the error grows. We then compare these error bounds with
relevant concepts in the literature such as exact and ordinary lumpability as
well as deflatability and aggregatability. These concepts define stricter than
necessary conditions to identify settings in which the aggregation error is
zero. We also consider possible algorithms for finding suitable aggregations
for which the formal error bounds are low, and we analyse first experiments
with these algorithms on different models.

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries
    1. 2.1 Aggregation of Markov chains
    2. 2.2 Exact aggregation
    3. 2.3 The choice of aggregation
    4. 2.4 Aggregated transition and generator matrices
    5. 2.5 The choice of the conditional distributions
  3. 3 Bounding the approximation error
    1. 3.1 Error bounds for DTMCs
      1. 3.1.1 Same bounds, alternative notation
    2. 3.2 Error bounds for CTMCs
    3. 3.3 When is the error bound 0?
    4. 3.4 Tightness of the error bound
    5. 3.5 The median-based scheme and the error bounds
  4. 4 Lumpability and aggregatability
  5. 5 Choosing the aggregates
    1. 5.1 No aggregation
    2. 5.2 Almost aggregatability
      1. 5.2.1 Implementation of SVD dir
    3. 5.3 $\varepsilon$-almost exact lumpability
    4. 5.4 Experiments
  6. 6 Conclusion & outlook

##  Introduction

State aggregation in dynamic systems has been studied extensively since the
1960s (see [14]). Due to the curse of dimensionality, models with large state
spaces are often computationally intractable without state space reduction,
and one basic reduction technique is to aggregate multiple states into a
single state in the aggregated model. Conditions under which an aggregated
Markov chain is again a Markov chain are well known (see strong and weak
lumpability in [11]), and various cases where exact transient or stationary
probabilities of the original model can be derived from an aggregated model
have been analysed (see, e.g. [3]).

However, formal error bounds for the approximation error when exact
aggregation is not possible have only been studied rarely. [3] already gave
upper and lower bounds for the transient distribution of a Markov chain which
are derived from an aggregated model. Much later, [1] has presented improved
bounds for the transient distributions of discrete-time Markov chains, which
can also be applied to continuous-time Markov chains via uniformisation. We
will extend the theory developed in [1] to support a more general way of
disaggregation in Section 3.1 and to the continuous-time domain without
falling back on uniformisation in Section 3.2.

Subsequently, we will analyse the cases where the error bounds are zero in
Section 3.3 and compare them with lumpability concepts from [3, 9, 2, 11] in
Section 3.3 and Section 4. Many of these types of lumpability imply, but are
not equivalent to the error bound being zero. Furthermore, we show that the
error bounds from [1] are tight in general in Section 3.4. In Section 5, we
present two different algorithms, one from [2] and a newly developed one based
on [4], with the goal to identify an aggregation resulting in low error
bounds. We then apply these algorithms to sample models and analyse the
results.

##  Preliminaries

###  Aggregation of Markov chains

We consider discrete-time and continuous-time Markov chains (DTMCs and CTMCs)
on the discrete, finite state space $S=\\{1,\ldots,n\\}$. The dynamics are
given by a stochastic transition matrix $P\in\mathbb{R}^{n\times n}$ for
DTMCs, where
$P(r,s)=\mathbb{P}\left[X_{k+1}=s\mathrel{}\middle|\mathrel{}X_{k}=r\right]$
if $X_{k}$ denotes the state of the DTMC at time $k$. For CTMCs, the dynamics
are defined via the generator matrix $Q\in\mathbb{R}^{n\times n}$, where
$Q(r,s)$ is the transition rate from $r$ to $s$, and $Q(r,r)$ is the negative
of the exit rate of state $r$, i.e. $Q(r,r)=-\sum_{s\neq r}Q(r,s)$. Given an
initial distribution $p_{0}\in\mathbb{R}^{n}$, the transient distribution of a
DTMC (respectively a CTMC) is given by $p_{k}=p_{0}P^{k}$ (respectively
$p_{t}=p_{0}e^{Qt}$) if we interpret $p_{k}$ as a row vector.

An aggregation of the state space consists of a set
$\Omega=\\{\Omega_{1},\ldots,\Omega_{m}\\}$ of $m$ aggregates, where $\Omega$
is a partition of the state space $S$ and thus an aggregate state
$\rho\in\Omega$ is a subset of $S$ which represents all states belonging to
one aggregate. We define the aggregation function $\omega:S\to\Omega$ which
maps a state $s$ to the aggregate to which $s$ belongs, i.e. $s\in\omega(s)$.
In order to approximate the dynamics of the original Markov chain, we can
define a transition or generator matrix on the aggregated state space, i.e. a
transition matrix $\Pi\in\mathbb{R}^{m\times m}$ for the aggregated state
space of a DTMC and a generator matrix $\Theta\in\mathbb{R}^{m\times m}$ for
CTMCs. $\Pi(\rho,\sigma)$ for $\rho,\sigma\in\Omega$ should be an
approximation of the probability to transition from one aggregate state into
another, that is, an approximation of
$\mathbb{P}\left[X_{k+1}\in\sigma\mathrel{}\middle|\mathrel{}X_{k}\in\rho\right]$.
Note that this probability may now depend on $k$, in contrast to the
probability
$\mathbb{P}\left[X_{k+1}=s\mathrel{}\middle|\mathrel{}X_{k}=r\right]$ for
$r,s\in S$, which is time-independent by the Markov property. However, we
consider only time-independent approximations of
$\mathbb{P}\left[X_{k+1}\in\sigma\mathrel{}\middle|\mathrel{}X_{k}\in\rho\right]$.
Simlarly, for CTMCs, we should have

$\displaystyle\Theta(\rho,\sigma)\approx\lim_{u\to
0}\frac{\mathbb{P}\left[X_{t+u}\in\sigma\mathrel{}\middle|\mathrel{}X_{t}\in\rho\right]}{u}$

where the right hand side may depend on $t$, but where we only consider time-
independent $\Theta$. Actually, $\Pi$ and $\Theta$ need not be stochastic or
generator matrices in general, but we only consider stochastic matrices $\Pi$
and generator matrices $\Theta$ because the aggregated Markov chain can then
be seen again as a Markov chain on a reduced state space. Such an
approximation of the dynamics provides the computational benefit of a lower
number of states for which the dynamics have to be evaluated, possibly at the
cost of lower precision.

The aggregation can also be applied to initial and transient distributions. We
call $\pi_{0}\in\mathbb{R}^{m}$, defined via
$\pi_{0}(\sigma)=\sum_{s\in\sigma}p_{0}(s)$ for $\sigma\in\Omega$, the
aggregated initial distribution and define aggregated transient distributions
via $\pi_{k}=\pi_{0}\Pi^{k}$ (discrete time) and $\pi_{t}=\pi_{0}e^{\Theta t}$
(continuous time). In order to obtain an approximation of the transient
probability for a given state in the original chain, we consider distributions
$\alpha_{\sigma}\in\mathbb{R}^{n}$ with support on $\sigma\in\Omega$. As a
shorthand, we write $\alpha(s):=\alpha_{\omega(s)}(s)$. The value $\alpha(s)$
should approximate the conditional probability of being in state $s$ when we
know that we are in aggregate $\omega(s)$, i.e. the probability
$\mathbb{P}\left[X_{k}=s\mathrel{}\middle|\mathrel{}X_{k}\in\omega(s)\right]$.
Again, this probability is in general dependent on time, but we only consider
time-independent approximations $\alpha$. We then arrive at an approximated
transient distribution $\widetilde{p}_{k}\in\mathbb{R}^{n}$ by the following
formula: $\widetilde{p}_{k}(s)=\alpha(s)\pi_{k}(\omega(s))$, or, for
continuous time, $\widetilde{p}_{t}(s)=\alpha(s)\pi_{t}(\omega(s))$.

|  DTMCs |  CTMCs
---|---|---
state space | $S=\\{1,\ldots,n\\}$
states | $r,s\in S$
transition / generator matrix | $P\in\mathbb{R}^{n\times n}$ | $Q\in\mathbb{R}^{n\times n}$
initial distribution | $p_{0}\in\mathbb{R}^{n}$
transient distribution | $p_{k}=p_{0}P^{k}$ | $p_{t}=p_{0}e^{Qt}$
aggregates |  $\Omega$ with $\left|\Omega\right|=m$ and $\Omega$ a partition of the state space
aggregation function |  $\omega:S\to\Omega$ such that $s\in\omega(s)$
aggregated states | $\rho,\sigma\in\Omega$
aggregated transition / generator matrix | $\Pi\in\mathbb{R}^{m\times m}$ | $\Theta\in\mathbb{R}^{m\times m}$
aggregated initial distribution | $\pi_{0}(\sigma)=\sum_{s\in\sigma}p_{0}(s)$
aggregated transient distribution | $\pi_{k}=\pi_{0}\Pi^{k}$ | $\pi_{t}=\pi_{0}e^{\Theta t}$
approximated conditional distributions |  for every $\sigma\in\Omega$: $\alpha_{\sigma}\in\mathbb{R}^{n}$ is a probability distribution with support on $\sigma$. Shorthand: $\alpha(s)=\alpha_{\omega(s)}(s)$
disaggregation | $\widetilde{p}_{k}(s)=\alpha(s)\pi_{k}(\omega(s))$ | $\widetilde{p}_{t}(s)=\alpha(s)\pi_{t}(\omega(s))$
Table 1: Notation overview

We further define the aggregation matrix $\Lambda$ and the disaggregation
matrix $A$ as follows:

$\displaystyle\Lambda=\begin{pmatrix}|&&|\\\
\mathbbm{1}_{\Omega_{1}}&\ldots&\mathbbm{1}_{\Omega_{m}}\\\
|&&|\end{pmatrix}\in\mathbb{R}^{n\times m},\qquad
A=\begin{pmatrix}\rotatebox[origin={c}]{90.0}{$|$}\,\alpha_{\Omega_{1}}\,\rotatebox[origin={c}]{90.0}{$|$}\\\
\vdots\\\
\rotatebox[origin={c}]{90.0}{$|$}\,\alpha_{\Omega_{m}}\,\rotatebox[origin={c}]{90.0}{$|$}\end{pmatrix}\in\mathbb{R}^{m\times
n}\qquad\textrm{(note: }A\Lambda=I\textrm{)}$

We set $\widetilde{P}=\Lambda\Pi A$ and $\widetilde{Q}=\Lambda\Theta A$. This
implies in particular that

$\displaystyle\widetilde{p}_{k}\widetilde{P}^{l}$
$\displaystyle=\underbrace{\widetilde{p}_{k}\Lambda}_{\pi_{k}}\Pi^{l}A=\pi_{k+l}A=\widetilde{p}_{k+l}$
$\displaystyle\widetilde{p}_{t}e^{\widetilde{Q}u}$
$\displaystyle=\widetilde{p}_{t}\sum_{k=0}^{\infty}\frac{u^{k}}{k!}\widetilde{Q}^{k}=\widetilde{p}_{t}\sum_{k=0}^{\infty}\frac{u^{k}}{k!}\Lambda\Theta^{k}A=\underbrace{\widetilde{p}_{t}\Lambda}_{\pi_{t}}\left(\sum_{k=0}^{\infty}\frac{u^{k}}{k!}\Theta^{k}\right)A=\pi_{t}e^{\Theta
u}A=\pi_{t+u}A=\widetilde{p}_{t+u}$

Hence, $\widetilde{P}$ (which is stochastic) describes the approximated
dynamics of the aggregated chain if we blow it up again to the original state
space $S$. It holds that
$\widetilde{P}(r,s)=\alpha(s)\Pi(\omega(r),\omega(s))$, so we essentially
approximate the transition probability $P(r,s)$ by the transition probability
from aggregate $\omega(r)$ to aggregate $\omega(s)$, weighted in addition with
the conditional probability $\alpha(s)$ of being in state $s$ within aggregate
$\omega(s)$.

With $\widetilde{Q}$, the case is a bit more complicated. We still have
$\widetilde{Q}(r,s)=\alpha(s)\Theta(\omega(r),\omega(s))$ and $\widetilde{Q}$
also describes the approximated dynamics in a sense, via the equation
$\widetilde{p}_{t}e^{\widetilde{Q}u}=\widetilde{p}_{t+u}$. However,
$\widetilde{Q}$ is no longer a generator matrix. While the row sums of
$\widetilde{Q}$ are $0$, the negative entries are no longer confined to the
diagonal. Instead, if two states $r,s$ both belong to aggregate $\sigma$, then
$\widetilde{Q}(r,s)=\alpha(s)\Theta(\sigma,\sigma)<0$ as long as $\Theta$ is
the generator matrix of an irreducible CTMC with more than one state. In this
case, we cannot see $\widetilde{Q}(r,s)$ as a rate approximating the
transition rate $Q(r,s)\geq 0$. $\widetilde{Q}(r,s)$ is rather (the
approximation of) the negative of the rate at which mass is exiting $\sigma$,
coming from $s$.

###  Exact aggregation

###### Definition 1

Given a partition $\Omega$ of the state space of a DTMC or CTMC, we call a
probability distribution $p$ on the state space $S$ compatible with
distributions $\alpha_{\sigma}$ with support on $\sigma\in\Omega$ if $p\Lambda
A=p$.

Compatibility of $p$ and the distributions $\alpha$ means that

$\displaystyle\alpha(s)=\frac{p(s)}{\sum_{s^{\prime}\in\omega(s)}p(s^{\prime})}\textrm{
for all }s\in S\textrm{ s.t. }\sum_{s^{\prime}\in\omega(s)}p(s^{\prime})>0$

###### Definition 2

We call an aggregation $\Omega=\\{\Omega_{1},\ldots,\Omega_{m}\\}$ of the
state space of a DTMC (respectively CTMC) with distributions
$\alpha_{\Omega_{1}},\ldots,\alpha_{\Omega_{m}}$ and aggregated transition
matrix $\Pi$ (respectively $\Theta$) dynamic-exact if $\Pi A=AP$ (respectively
$\Theta A=AQ$).

If the initial distribution $p_{0}$ is further compatible with the $\alpha$
distributions, i.e. if $\widetilde{p}_{0}=p_{0}\Lambda A=p_{0}$, then we call
the aggregation exact.

If $\Omega=\\{\Omega_{1},\ldots,\Omega_{m}\\}$,
$\alpha_{\Omega_{1}},\ldots,\alpha_{\Omega_{m}}$ and $\Pi$ are an exact
aggregation, then $\widetilde{p}_{k}=p_{k}$ for all $k$. This can be proven
easily via induction and will also be shown later in 7. Indeed, $\Pi A=AP$
implies $AP^{k}=APP^{k-1}=\Pi APP^{k-2}=\ldots=\Pi^{k}A$. Hence, for any
$p_{0}$ compatible with the $\alpha$ distributions, we have

$\displaystyle p_{k}$ $\displaystyle=p_{0}P^{k}=p_{0}\Lambda
AP^{k}=p_{0}\Lambda\Pi^{k}A=\pi_{0}\Pi^{k}A=\pi_{k}A=\widetilde{p}_{k}$ (2.1)

Note that $\widetilde{p}_{0}=p_{0}\Lambda A$ is not necessarily equal to
$p_{0}$, even if $P=\widetilde{P}$. It is true, however, that
$\widetilde{p}_{k}=\widetilde{p}_{0}P^{k}$ if the aggregation is dynamic-exact
(this follows from 7). Note the difference to
$\widetilde{p}_{k}=\widetilde{p}_{0}\widetilde{P}^{k}$, which is always true.
This is the reason for the term “dynamic-exact”: the step dynamics are
correctly represented by $\widetilde{P}$, but the initial distribution might
be wrongly approximated. The corresponding statements also hold for CTMCs.

###### Remark

Note that the condition $\widetilde{P}=\Lambda\Pi A=P$ is stronger than (i.e.
implies) the condition $\Pi A=AP$ given in 2. Essentially, $\widetilde{P}=P$
requires that any probability distribution on the states of the DTMC is pushed
forward correctly by one step by applying $\widetilde{P}$. $\Pi A=AP$, on the
other hand, only requires that $\widetilde{P}=\Lambda AP$ which means that
only those probability distributions which are compatible with the conditional
distributions $\alpha$ are pushed forward correctly. For CTMCs, the condition
$\widetilde{Q}=Q$ can actually never hold if the partition $\Omega$ is not the
partition into single states (i.e. if it is true that $m<n$), and if the CTMC
is irreducible. It might be argued, therefore, that [2, Definition 2.1] is too
restrictive since it requires $\Lambda\Pi A=P$ in order to call $P$
aggregatable.

The condition $\Pi A=AP$ has appeared in the literature before. Indeed,
equation (4) on page 135 of [11] states that, if $\Pi$ is set as in (2.2)
below, and if the $\alpha$ distributions are compatible with the initial
distribution $p_{0}$, then $\Pi A=AP$ implies weak lumpability of the DTMC. A
DTMC is called weakly lumpable for a given partition $\Omega$ if there exists
an initial distribution $p_{0}$ such that the process $Y_{k}$, defined by
$Y_{k}=\sigma\in\Omega\iff X_{k}\in\sigma$, is a Markov chain. For such an
initial distribution, the probabilities $\pi_{k}(\sigma)$ are then exactly
equal to
$\mathbb{P}\left[Y_{k}=\sigma\right]=\mathbb{P}\left[X_{k}\in\sigma\right]$.
However, the concept of weak lumpability makes no statement about whether the
probability $\mathbb{P}\left[X_{k}=s\right]$ for $s\in\sigma$ can be
accurately derived from the knowledge of
$\mathbb{P}\left[X_{k}\in\sigma\right]$.

Again under the condition that $\Pi$ is set as in (2.2), [12, Definition 2.2]
defined the matrix $P$ to be $A$-lumpable if $\Pi A=AP$. In the subsequent
remarks, [12] then noted that, given $\pi_{k}$, an exact recovery of the
probabilities $p_{k}(s)$ is possible provided that the initial distribution
$p_{0}$ is compatible with the $\alpha$ distributions. Actually, [12, equation
(2.4)] corresponds almost exactly to (2.1).

Furthermore, an exact aggregation is called backward bisimulation of type 2 in
[5, Definition 4.3], yet again if $\Pi$ is set as in (2.2). [5] already
remarks that the coarsest exact aggregation cannot be found by the usual
partition refinement approach, which we will see later in Section 5.3.
Efficiently finding partitions $\Omega$ (and distributions $\alpha$) which are
an exact aggregation for a given DTMC or CTMC therefore currently still is an
open research problem to the best of our knowledge.

###  The choice of aggregation

We consider different possibilities for choosing the aggregates in detail
later in this paper. The quality of the approximated dynamics given by the
aggregated chain is strongly dependent on the choice of aggregation. For now,
however, we assume that we are given an aggregation in the form of
$\Omega=\\{\Omega_{1},\ldots,\Omega_{m}\\}$, and we want to choose $\Pi$
(respectively $\Theta$) and $\alpha$ in a way which results in a good
approximation of the dynamics of the original chain.

###  Aggregated transition and generator matrices

We consider the following aggregated transition probabilities and rates for
$\rho,\sigma\in\Omega$ where it may hold that $\rho=\sigma$:

$\displaystyle\Pi(\rho,\sigma)=\sum_{r\in\rho}\alpha(r)\sum_{s\in\sigma}P(r,s)\qquad\qquad\Theta(\rho,\sigma)=\sum_{r\in\rho}\alpha(r)\sum_{s\in\sigma}Q(r,s)$
(2.2)

The probability
$\mathbb{P}\left[X_{k+1}\in\sigma\mathrel{}\middle|\mathrel{}X_{k}\in\rho\right]$
is thus approximated via the weighted average (with weights $\alpha(r)$) of
the probability to transition from a single state $r\in\rho$ into any of the
states in $\sigma$. The reasoning behind this is as follows: we assume that if
we are in aggregate $\rho$, the probability to be in state $r\in\rho$ is given
(approximately) by $\alpha(r)$. For the continuous-time case, we take the
weighted average of the transition rates. In matrix notation, setting $\Pi$
and $\Theta$ as in (2.2) corresponds to $\Pi=AP\Lambda$ and
$\Theta=AQ\Lambda$. Note that $\Pi$ is again stochastic, and $\Theta$ is again
a generator matrix:

$\displaystyle\sum_{\sigma\in\Omega}\Pi(\rho,\sigma)=\sum_{r\in\rho}\alpha(r)\underbrace{\sum_{s\in
S}P(r,s)}_{1}=1\qquad\qquad\sum_{\sigma\in\Omega}\Theta(\rho,\sigma)=\sum_{r\in\rho}\alpha(r)\underbrace{\sum_{s\in
S}Q(r,s)}_{0}=0$

It is easy to see that we also have $\Pi(\rho,\sigma)\geq 0$ for all
$\rho,\sigma$ as well as $\Theta(\rho,\sigma)\geq 0$ for $\rho\neq\sigma$ and
$\Theta(\sigma,\sigma)\leq 0$. Different choices for $\Pi$ are discussed in
[1], and the above choice yields good approximations in terms of the transient
distribution in the experiments done in [1]. Only a so-called “median-based
scheme” (see [1, equation (21)]) performs better in some settings. As (2.2)
provides the additional benefit that the aggregated chain with the aggregated
dynamics is again a Markov chain (as $\Pi$ is stochastic), we will mostly
consider this choice of $\Pi$ and $\Theta$ in the present paper, which is also
a very intuitive choice.

The median-based scheme might perform better and can also be used with the
distributions $\alpha_{\sigma}$, and not just with the implicitly assumed
uniform distributions on the aggregates considered in [1]. It works as
follows: assume that the partition $\Omega$ and the distributions
$\alpha_{\sigma}$ are fixed. We then set $\Pi(\rho,\sigma)$ to the weighted
median of the (multi-)set

$\displaystyle a_{s}$
$\displaystyle=\frac{1}{\alpha(s)}\sum_{r\in\rho}\alpha(r)P(r,s)\qquad(s\in\sigma\textrm{
with }\alpha(s)>0)$ (2.3)

weighted with the weights $\alpha(s)$ for every state $s\in\sigma$. The
weighted median of the (multi-)set $\\{a_{s}\;|\;s\in\sigma\textrm{ with
}\alpha(s)>0\\}$ is the element $a_{s_{\textrm{med}}}$ with the property that

$\displaystyle\sum_{\begin{subarray}{c}s\in\sigma:\alpha(s)>0\\\
a_{s}<a_{s_{\textrm{med}}}\end{subarray}}\alpha(s)\leq\frac{1}{2}\quad$
$\displaystyle\textrm{and}\quad\sum_{\begin{subarray}{c}s\in\sigma:\alpha(s)>0\\\
a_{s}>a_{s_{\textrm{med}}}\end{subarray}}\alpha(s)\leq\frac{1}{2}$ (2.4)

It can be that there are multiple elements in the (multi-)set
$\\{a_{s}\;|\;s\in\sigma\textrm{ with }\alpha(s)>0\\}$ which satisfy (2.4). In
such a case, the weighted median is non-unique, and any of the elements
satisfying (2.4) may be chosen. For a fixed partition $\Omega$ and fixed
$\alpha_{\sigma}$ distributions, setting $\Pi(\rho,\sigma)$ to this weighted
median minimizes the error bounds which are introduced later in Section 3.1
and Section 3.2 (see 10). We usually use the more intuitive (2.2) to determine
$\Pi$ and $\Theta$, but some of the experiments at the end compare (2.2) to
the median-based scheme.

###### Remark

Using the median-based scheme can result in matrices $\Pi$ and $\Theta$ which
are no longer stochastic or a generator. Therefore, the aggregated transient
distributions $\pi_{k}$ and $\pi_{t}$ will, in fact, no longer be probability
distributions in this case. The error bounds given in Section 3.1 and Section
3.2 are still valid, however.

###  The choice of the conditional distributions

Another question is how to choose the conditional distributions
$\alpha_{\sigma}$. For DTMCs, the following definitions provided good results,
and are compatible with the aggregation techniques considered later. We
discuss this compatibility later in detail. We consider two possibilities for
setting $\alpha(s)$ for DTMCs:

  * •

The first possibility, called proportional $\alpha$, is given by

$\displaystyle\alpha(s)$ $\displaystyle=\frac{\sum_{r\in S}P(r,s)}{\sum_{r\in
S}\sum_{s^{\prime}\in\omega(s)}P(r,s^{\prime})}$ (2.5)

We set $\alpha(s)$ to the fraction of the incoming probabilities from all
states to $s$ divided by the total incoming probabilities to any state in the
aggregate of $s$. This is the same as the probability of being in state $s$,
conditioned on being in the aggregate of $s$, after the Markov chain took a
single step, starting with a uniform distribution. Intuitively, the
distributions $\alpha_{\sigma}$ should be approximations of this type of
conditional probabilities, with the exception that we do not necessarily start
with a uniform distribution. Note that $\alpha(s)$ is well-defined (the
denominator is never $0$) if we assume irreducibility, which also implies that
$\alpha(s)>0$ for all $s\in S$.

If we do not assume irreducibility, the denominator might be $0$. This would
imply, however, that none of the states in $\omega(s)$ can be reached from any
other state in the chain. In such a case, we could set
$\alpha(s)=\frac{p_{0}(s)}{\sum_{s^{\prime}\in\omega(s)}p_{0}(s^{\prime})}$ if
this is well defined and otherwise just choose an arbitrary $\alpha(s)$. The
probability to transition into aggregate $\omega(s)$ will be $0$ if the
denominator in (2.5) is $0$, and hence $\pi_{k}(\omega(s))>0$ can only happen
for $k=0$ (for any sensible choice of $\Pi$).

  * •

The second possibility, called uniform $\alpha$, is given by

$\displaystyle\alpha(s)$ $\displaystyle=\frac{1}{\left|\omega(s)\right|}$
(2.6)

This results in distributions $\alpha_{\sigma}$ which are uniform
distributions on the respective aggregate $\sigma$. This choice of $\alpha$ is
discussed later in more detail.

An alternative to these two possibilities would be to consider the aggregate
of $s$ as a Markov chain on its own and compute its stationary distribution,
which is then used for the conditional distribution $\alpha_{\sigma}$. This
makes sense in particular if the transition probabilities between aggregates
are low and the probability to stay for a long time within an aggregate is
high (which is the case if $P$ is nearly completely decomposable, see e.g.
[14]). We will, however, use the definitions given in (2.5) and (2.6) for
their simplicity and computational efficiency. We also compare the different
possibilities in experiments.

For CTMCs, the formula (2.5) cannot be used. The problem here is that the sum
$\sum_{r\in S}Q(r,s)$ could be negative since the diagonal entries of the
matrix $Q$ are negative, and while the row sums of $Q$ are all $0$, the column
sums can take arbitrary values. To adapt (2.5) for CTMCs, one possibility
would be to set $\alpha$ as follows:

$\displaystyle\alpha(s)$
$\displaystyle=\frac{\displaystyle\sum_{\begin{subarray}{c}r\in S\\\
r\notin\omega(s)\end{subarray}}Q(r,s)}{\displaystyle\sum_{\begin{subarray}{c}r\in
S\\\
r\notin\omega(s)\end{subarray}}\sum_{s^{\prime}\in\omega(s)}Q(r,s^{\prime})}$
(2.7)

Note that (2.7) only makes sense if there is more than one aggregate. See the
discrete-time case for what to do when the denominator in (2.7) is $0$.

###### Remark

As will become apparent later in the proof of 14, only excluding the diagonal
elements in (2.7), i.e. setting

$\displaystyle\alpha(s)$
$\displaystyle=\frac{\displaystyle\sum_{\begin{subarray}{c}r\in S\\\ r\neq
s\end{subarray}}Q(r,s)}{\displaystyle\sum_{r\in
S}\sum_{\begin{subarray}{c}s^{\prime}\in\omega(s)\\\ s^{\prime}\neq
r\end{subarray}}Q(r,s^{\prime})}$

results in some undesirable properties for the aggregation schemes discussed
later.

##  Bounding the approximation error

###  Error bounds for DTMCs

In this section, we follow [1] to derive formal error bounds for the
difference between the actual transient distribution $p_{k}$ of the Markov
chain and the approximation $\widetilde{p}_{k}$ of this distribution, obtained
after state space aggregation. [3, Theorem 11] gave similar bounds, but the
bounds in [1] are slightly better. For now, we consider the aggregation
$\Omega$, the aggregated transition probabilities $\Pi$, and the distributions
$\alpha_{\sigma}$ as given. They can be arbitrary as long as $\Omega$
partitions the state space into $m$ aggregates, $\Pi$ is a matrix in
$\mathbb{R}^{m\times m}$, and the $\alpha_{\sigma}$ are probability
distributions with support on the aggregate $\sigma\in\Omega$. The
calculations below do not use the particular forms of $\Pi$ and $\alpha$ given
in (2.2) and (2.5), but they hold for all choices of $\Pi$ and $\alpha$.

The following is very similar to what was already done in [1], with the
difference that [1] did not use distributions $\alpha_{\sigma}$, but simply
set $\alpha(s)=\frac{1}{\left|\omega(s)\right|}$ implicitly. We include the
calculations here to demonstrate that arbitrary choices of $\alpha$ do not
significantly change the error bounds derived in [1], and we will later reuse
them for the continuous-time case.

We call $e_{k}=\widetilde{p}_{k}-p_{k}$ the error after step $k$, and we want
to bound $\left\lVert
e_{k}\right\rVert_{1}=\sum_{i=1}^{n}\left|e_{k}(i)\right|$ where $e_{k}(i)$ is
the $i$-th entry of $e_{k}\in\mathbb{R}^{n}$. To this end, note that

$\displaystyle\begin{split}e_{k}&=\widetilde{p}_{k}-p_{k}=\widetilde{p}_{k-1}\cdot\widetilde{P}-p_{k-1}\cdot
P=\widetilde{p}_{k-1}\cdot\left(\widetilde{P}-P+P\right)-p_{k-1}\cdot P\\\
&=\widetilde{p}_{k-1}\cdot\left(\widetilde{P}-P\right)+\underbrace{\left(\widetilde{p}_{k-1}-p_{k-1}\right)}_{e_{k-1}}\cdot
P\end{split}$ $\displaystyle\implies\left\lVert e_{k}\right\rVert_{1}$
$\displaystyle\leq\left\lVert\widetilde{p}_{k-1}\cdot\left(\widetilde{P}-P\right)\right\rVert_{1}+\left\lVert
e_{k-1}\cdot P\right\rVert_{1}$ (3.1)

To further simplify the term above, we will use the following lemma.

###### Lemma 3

Let $e\in\mathbb{R}^{k}$ be an arbitrary (row) vector and
$P\in\mathbb{R}^{k\times k}$ be an arbitrary stochastic matrix. Then

$\displaystyle\left\lVert e\cdot P\right\rVert_{1}$
$\displaystyle\leq\left\lVert e\right\rVert_{1}$

###### Proof

It holds that

$\displaystyle\left\lVert e\cdot P\right\rVert_{1}$
$\displaystyle=\sum_{s=1}^{k}\left|\sum_{r=1}^{k}e(r)\cdot
P(r,s)\right|\leq\sum_{s=1}^{k}\sum_{r=1}^{k}\left|e(r)\right|\cdot P(r,s)$
$\displaystyle=\sum_{r=1}^{k}\left|e(r)\right|\cdot\sum_{s=1}^{k}P(r,s)=\sum_{r=1}^{k}\left|e(r)\right|=\left\lVert
e\right\rVert_{1}$

As a consequence of 3 and (3.1), we have

$\displaystyle\left\lVert e_{k}\right\rVert_{1}$
$\displaystyle\leq\underbrace{\left\lVert
e_{k-1}\right\rVert_{1}}_{\textrm{previous
error}}+\underbrace{\left\lVert\widetilde{p}_{k-1}\cdot\left(\widetilde{P}-P\right)\right\rVert_{1}}_{\textrm{error
from using approximated transition probabilities}}$

We bound the second term as follows (see [1, page 15-17]):

$\displaystyle\begin{split}\left\lVert\widetilde{p}_{k-1}\cdot\left(\widetilde{P}-P\right)\right\rVert_{1}&=\sum_{\sigma\in\Omega}\sum_{s\in\sigma}\left|\sum_{\rho\in\Omega}\sum_{r\in\rho}\widetilde{p}_{k-1}(r)\cdot\left(\widetilde{P}(r,s)-P(r,s)\right)\right|\\\
&=\sum_{\sigma\in\Omega}\sum_{s\in\sigma}\left|\sum_{\rho\in\Omega}\sum_{r\in\rho}\alpha(r)\pi_{k-1}(\rho)\cdot\left(\alpha(s)\Pi(\rho,\sigma)-P(r,s)\vphantom{\widetilde{P}}\right)\right|\\\
&=\sum_{\sigma\in\Omega}\sum_{s\in\sigma}\left|\sum_{\rho\in\Omega}\pi_{k-1}(\rho)\cdot\left(\alpha(s)\Pi(\rho,\sigma)\sum_{r\in\rho}\alpha(r)-\sum_{r\in\rho}\alpha(r)P(r,s)\right)\right|\\\
&=\sum_{\sigma\in\Omega}\sum_{s\in\sigma}\left|\sum_{\rho\in\Omega}\pi_{k-1}(\rho)\cdot\left(\alpha(s)\Pi(\rho,\sigma)-\sum_{r\in\rho}\alpha(r)P(r,s)\right)\right|\\\
&\leq\sum_{\sigma\in\Omega}\sum_{s\in\sigma}\sum_{\rho\in\Omega}\pi_{k-1}(\rho)\cdot\left|\alpha(s)\Pi(\rho,\sigma)-\sum_{r\in\rho}\alpha(r)P(r,s)\right|\\\
&=\sum_{\rho\in\Omega}\pi_{k-1}(\rho)\cdot\underbrace{\sum_{\sigma\in\Omega}\sum_{s\in\sigma}\left|\alpha(s)\Pi(\rho,\sigma)-\sum_{r\in\rho}\alpha(r)P(r,s)\right|}_{\displaystyle=:\tau(\rho)}\end{split}$
(3.2)

In particular, it follows that

$\displaystyle\begin{split}\left\lVert e_{k}\right\rVert_{1}&\leq\left\lVert
e_{k-1}\right\rVert_{1}+\sum_{\rho\in\Omega}\pi_{k-1}(\rho)\cdot\tau(\rho)\leq\ldots\\\
&\leq\left\lVert
e_{0}\right\rVert_{1}+\sum_{i=0}^{k-1}\underbrace{\sum_{\rho\in\Omega}\pi_{i}(\rho)\cdot\tau(\rho)}_{\leq\max_{\rho\in\Omega}\tau(\rho)}\leq\left\lVert
e_{0}\right\rVert_{1}+k\cdot\max_{\rho\in\Omega}\tau(\rho)\end{split}$ (3.3)

Calculating $\tau(\rho)$ for every aggregate $\rho\in\Omega$ once thus allows
us to bound the error $\left\lVert e_{k}\right\rVert_{1}$ of the approximated
transient probabilities at time $k$ using just the error at time $0$ and the
aggregated transition probabilities $\pi_{i}$ at times $i=0,\ldots,k-1$. As a
shorthand notation, we will use

$\displaystyle\left<\pi_{k},\tau\right>:=\sum_{\rho\in\Omega}\pi_{k}(\rho)\cdot\tau(\rho)$

to denote the error made at the time step $k\to k+1$. For an even simpler but
less precise bound, $\max_{\rho\in\Omega}\tau(\rho)$ gives an upper bound for
the error growth in every step which gives an instant linear error bound on
the error after $k$ steps.

Note that we can also write (3.3) in another way using the matrices $\Pi A$
and $AP$ which are both matrices in $\mathbb{R}^{m\times n}$. We have that,
with $\omega(s)=\sigma$:

$\displaystyle\alpha(s)\Pi(\rho,\sigma)=\underbrace{(\Pi
A)}_{\in\mathbb{R}^{m\times n}}(\rho,s)\qquad\textrm{ and
}\qquad\sum_{r\in\rho}\alpha(r)P(r,s)=(AP)(\rho,s)$

Hence,

$\displaystyle\left\lVert\Pi A-AP\right\rVert_{\infty}$
$\displaystyle=\max_{\rho\in\Omega}\sum_{s\in S}\left|(\Pi
A)(\rho,s)-(AP)(\rho,s)\right|$
$\displaystyle=\max_{\rho\in\Omega}\underbrace{\sum_{\sigma\in\Omega}\sum_{s\in\sigma}\left|\alpha(s)\Pi(\rho,\sigma)-\sum_{r\in\rho}\alpha(r)P(r,s)\right|}_{\tau(\rho)}$

We can plug this into (3.3) and arrive at

$\displaystyle\left\lVert e_{k}\right\rVert_{1}$ $\displaystyle\leq\left\lVert
e_{0}\right\rVert_{1}+k\cdot\left\lVert\Pi
A-AP\right\rVert_{\infty}=\left\lVert\pi_{0}A-p_{0}\right\rVert_{1}+k\cdot\left\lVert\Pi
A-AP\right\rVert_{\infty}$ (3.4)

In fact, we should aim to minimize the right hand term in (3.4) over
$\pi_{0}$, $\Pi$ and $A$ for an optimal result. This is the subject of ongoing
research. An approach similar to the one in [6] seems promising. Actually, we
can rewrite the whole process of bounding $\left\lVert e_{k}\right\rVert_{1}$,
as shown below.

####  Same bounds, alternative notation

Only in this subsection, we consider the following setting:

  * •

We are given a DTMC via its transition matrix $P\in\mathbb{R}^{n\times n}$ and
initial distribution $p_{0}\in\mathbb{R}^{n}$ (a row vector). We do not assume
that the DTMC possesses any special properties, we only require that $P$ is a
stochastic matrix and that $p_{0}$ is a probability distribution.

  * •

We want to approximate the behavior of this DTMC in a lower-dimensional
“state” space with $m$ “states”. This is done as follows:

  * •

We choose an initial (row) vector $\pi_{0}\in\mathbb{R}^{m}$. We do not
restrict the choice of $\pi_{0}$ in any way – it might have negative entries,
and it need not be a distribution.

  * •

We choose a “transition” matrix $\Pi\in\mathbb{R}^{m\times m}$ again without
any restrictions. We then set $\pi_{k}:=\pi_{0}\Pi^{k}$.

  * •

We choose a disaggregation matrix $A\in\mathbb{R}^{m\times n}$ again without
any restrictions and approximate $p_{k}=p_{0}P^{k}$ by
$\widetilde{p}_{k}:=\pi_{k}A=\pi_{0}\Pi^{k}A$.

In applications, it can be sensible to choose a probability distribution
$\pi_{0}$, a stochastic matrix $\Pi$, and a matrix $A$ with non-negative
entries and with each row summing to $1$. However, we explicitly do not
require these assumptions here because the following calculations hold in the
general case.

We now want to give a bound on $\left\lVert e_{k}\right\rVert_{1}$, where we
set $e_{k}=\widetilde{p}_{k}-p_{k}$, as before. We have

$\displaystyle\left\lVert e_{k}\right\rVert_{1}$
$\displaystyle=\left\lVert\pi_{0}\Pi^{k}A-p_{0}P^{k}\right\rVert_{1}=\left\lVert\vphantom{\left(\Pi^{k-1}A-p_{0}P^{k-1}\right)}\right.\underbrace{\left(\pi_{0}\Pi^{k-1}A-p_{0}P^{k-1}\right)}_{e_{k-1}}\cdot
P+\underbrace{\pi_{0}\Pi^{k-1}}_{\pi_{k-1}}\left(\Pi
A-AP\right)\left.\vphantom{\left(\Pi^{k-1}A-p_{0}P^{k-1}\right)}\right\rVert_{1}$
$\displaystyle\leq\left\lVert e_{k-1}\cdot
P\right\rVert_{1}+\left\lVert\pi_{k-1}\cdot\left(\Pi
A-AP\right)\right\rVert_{1}\overset{\textrm{\autoref{lem:stoch_1norm_mult}}}{\leq}\left\lVert
e_{k-1}\right\rVert_{1}+\left\lVert\pi_{k-1}\cdot\left(\Pi
A-AP\right)\right\rVert_{1}$

For general row vectors $v\in\mathbb{R}^{m}$ and matrices
$A\in\mathbb{R}^{m\times n}$, it holds that

$\displaystyle\left\lVert vA\right\rVert_{1}$
$\displaystyle=\sum_{s=1}^{n}\left|\sum_{\rho=1}^{m}v(\rho)\cdot
A(\rho,s)\right|\leq\sum_{s=1}^{n}\sum_{\rho=1}^{m}\left|v(\rho)\right|\cdot\left|A(\rho,s)\right|$
$\displaystyle=\sum_{\rho=1}^{m}\left|v(\rho)\right|\sum_{s=1}^{n}\left|A(\rho,s)\right|\leq\left(\sum_{\rho=1}^{m}\left|v(\rho)\right|\right)\left(\max_{\rho=1,\ldots,m}\sum_{s=1}^{n}\left|A(\rho,s)\right|\right)=\left\lVert
v\right\rVert_{1}\cdot\left\lVert A\right\rVert_{\infty}$

Hence:

$\displaystyle\left\lVert\pi_{k-1}\cdot\left(\Pi A-AP\right)\right\rVert_{1}$
$\displaystyle\leq\left\lVert\pi_{k-1}\right\rVert_{1}\cdot\left\lVert\Pi
A-AP\right\rVert_{\infty}=\left\lVert\pi_{0}\cdot\Pi^{k-1}\right\rVert_{1}\cdot\left\lVert\Pi
A-AP\right\rVert_{\infty}$
$\displaystyle\leq\left\lVert\pi_{0}\right\rVert_{1}\cdot\left\lVert\Pi\right\rVert^{k-1}_{\infty}\cdot\left\lVert\Pi
A-AP\right\rVert_{\infty}$

In total, we thus get

$\displaystyle\left\lVert e_{k}\right\rVert_{1}$ $\displaystyle\leq\left\lVert
e_{0}\right\rVert_{1}+\sum_{j=1}^{k}\left\lVert\pi_{0}\right\rVert_{1}\cdot\left\lVert\Pi\right\rVert^{j-1}_{\infty}\cdot\left\lVert\Pi
A-AP\right\rVert_{\infty}$
$\displaystyle=\left\lVert\pi_{0}A-p_{0}\right\rVert_{1}+\left\lVert\pi_{0}\right\rVert_{1}\cdot\left\lVert\Pi
A-AP\right\rVert_{\infty}\cdot\sum_{j=0}^{k-1}\left\lVert\Pi\right\rVert^{j}_{\infty}$
$\displaystyle=\left\lVert\pi_{0}A-p_{0}\right\rVert_{1}+\left\lVert\pi_{0}\right\rVert_{1}\cdot\left\lVert\Pi
A-AP\right\rVert_{\infty}\cdot\begin{cases}\frac{1-\left\lVert\Pi\right\rVert^{k}_{\infty}}{1-\left\lVert\Pi\right\rVert_{\infty}}&\textrm{
if }\left\lVert\Pi\right\rVert_{\infty}\neq 1\\\ k&\textrm{
otherwise}\end{cases}$

The above equation simplifies to (3.4) if $\pi_{0}$ is a probability
distribution and if $\Pi$ is stochastic.

###  Error bounds for CTMCs

We will now extend the above setting to continuous-time Markov chains in the
following way: we set

$\displaystyle\tau(\rho)$
$\displaystyle:=\sum_{\sigma\in\Omega}\sum_{s\in\sigma}\left|\alpha(s)\Theta(\rho,\sigma)-\sum_{r\in\rho}\alpha(r)Q(r,s)\right|$
(3.5)

which corresponds exactly to the definition in the case of discrete time, only
with $P$ replaced by $Q$, the generator matrix, and with $\Pi$ replaced by
$\Theta$. In [1], the following was already shown: define $\Theta$ as in (2.2)
and calculate $\tau$ as above in (3.5). Consider a uniformisation of the
original CTMC with uniformisation rate $q$, apply the same aggregation to the
resulting DTMC and set the aggregated transition probabilities of the DTMC as
indicated in (2.2). Then, the error growth in step $k\to k+1$ in the DTMC can
be bounded by $\frac{1}{q}\cdot\left<\pi_{k},\tau\right>$ where $\tau$ is
derived for the aggregated CTMC as in (3.5).

We now want to drop the detour via the uniformisation and claim that we can
directly bound the approximation error of the transient distribution in
continuous time via the factors $\tau$ which can now be interpreted as a rate
of error growth (instead of error growth per step).

###### Theorem 4

Let $e_{t}$ be the vector of component-wise error of the approximated
transient distribution at time $t$, i.e.

$\displaystyle e_{t}=\widetilde{p}_{t}-p_{t}=\pi_{0}\cdot e^{\Theta
t}A-p_{0}\cdot e^{Qt}$

Then, the following hold:

  1. (i)

We have the following bound for the error at time $t$:

$\displaystyle\left\lVert e_{t}\right\rVert_{1}$ $\displaystyle\leq\left\lVert
e_{0}\right\rVert_{1}+\int_{0}^{t}\left<\pi_{s},\tau\right>\;\mathrm{d}s\qquad\qquad\textrm{where
}\left<\pi_{s},\tau\right>=\sum_{\rho\in\Omega}\pi_{s}(\rho)\cdot\tau(\rho)$
$\displaystyle\textrm{and }\quad\left\lVert e_{t}\right\rVert_{1}$
$\displaystyle\leq\left\lVert
e_{0}\right\rVert_{1}+t\cdot\max_{\rho\in\Omega}\tau(\rho)\quad\textrm{ as
well as }\quad\left\lVert e_{t}\right\rVert_{1}\leq\left\lVert
e_{0}\right\rVert_{1}+t\cdot\left\lVert\Theta A-AQ\right\rVert_{\infty}$

  2. (ii)

$\left\lVert e_{t}\right\rVert_{1}$ is absolutely continuous in $t$ and thus
almost everywhere (a.e.) differentiable, and

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\left\lVert e_{t}\right\rVert_{1}$
$\displaystyle\leq\left<\pi_{t},\tau\right>\textrm{ a.e.}\quad\textrm{ and
}\quad\limsup_{u\to 0}\frac{\left\lVert e_{t+u}\right\rVert_{1}-\left\lVert
e_{t}\right\rVert_{1}}{u}\leq\left<\pi_{t},\tau\right>\textrm{ everywhere}$

Before being able to prove the above theorem, we need another lemma.

###### Lemma 5

Assume that $f:\mathbb{R}\to\mathbb{R}$ is differentiable in $0$, and that
$f(0)=0$. Then

$\displaystyle\limsup_{u\to
0}\frac{\left|f(u)\right|}{u}=\lim_{\begin{subarray}{c}u\to 0\\\
u>0\end{subarray}}\frac{\left|f(u)\right|}{u}$
$\displaystyle=\left|f^{\prime}(0)\right|$

###### Proof

We have:

$\displaystyle\left|f^{\prime}(0)\right|$ $\displaystyle=\left|\lim_{u\to
0}\frac{f(u)}{u}\right|=\lim_{u\to
0}\frac{\left|f(u)\right|}{\left|u\right|}=\lim_{\begin{subarray}{c}u\to 0\\\
u>0\end{subarray}}\frac{\left|f(u)\right|}{u}=\limsup_{u\to
0}\frac{\left|f(u)\right|}{u}$

######  Proof of (Theorem 4)

As ii immediately implies i, it suffices to prove ii. Note that, as shown for
the discrete-time case in (3.4), we also have
$\max_{\rho\in\Omega}\tau(\rho)=\left\lVert\Theta A-AQ\right\rVert_{\infty}$
by the exact same calculation in the continuous-time case, which explains the
last statement in i.

First, note the following: every component of $e_{t}$ is continuously
differentiable in $t$, as both $p_{t}$ and $\widetilde{p}_{t}$ are
continuously differentiable with respect to $t$. Indeed, calculating the
derivative of all components of $e_{t}$ simultaneously, we get (here again, we
interpret $\widetilde{p}_{t}$ and $p_{t}$ as row vectors)

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}(\widetilde{p}_{t}-p_{t})$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}t}\left(\widetilde{p}_{0}e^{\widetilde{Q}t}-p_{0}e^{Qt}\right)=\widetilde{p}_{0}e^{\widetilde{Q}t}\widetilde{Q}-p_{0}e^{Qt}Q=\widetilde{p}_{t}\widetilde{Q}-p_{t}Q$
(3.6)
$\displaystyle\left\lVert\frac{\mathrm{d}}{\mathrm{d}t}(\widetilde{p}_{t}-p_{t})\right\rVert_{1}$
$\displaystyle=\left\lVert\widetilde{p}_{t}\widetilde{Q}-p_{t}Q\right\rVert_{1}\leq\left\lVert\widetilde{p}_{t}\widetilde{Q}\right\rVert_{1}+\left\lVert
p_{t}Q\right\rVert_{1}\leq\left|S\right|\cdot\left(\max_{r,s\in
S}\left|\widetilde{Q}(r,s)\right|+\max_{r,s\in S}\left|Q(r,s)\right|\right)$

As every component of $e_{t}$ is continuously differentiable with bounded
derivative, $\left\lVert e_{t}\right\rVert_{1}$ is absolutely continuous and
differentiable almost everywhere (see, e.g., [13, Section 5.4 on page 108]).

We have

$\displaystyle\left\lVert e_{t+u}\right\rVert_{1}$
$\displaystyle=\left\lVert\widetilde{p}_{t+u}-p_{t+u}\right\rVert_{1}=\left\lVert\widetilde{p}_{t}\left(e^{\widetilde{Q}u}-e^{Qu}\right)+\left(\widetilde{p}_{t}-p_{t}\right)e^{Qu}\right\rVert_{1}$
$\displaystyle\leq\left\lVert\widetilde{p}_{t}\left(e^{\widetilde{Q}u}-e^{Qu}\right)\right\rVert_{1}+\left\lVert\left(\widetilde{p}_{t}-p_{t}\right)e^{Qu}\right\rVert_{1}$
$\displaystyle\overset{\circledcirc}{\leq}\left\lVert\widetilde{p}_{t}\left(e^{\widetilde{Q}u}-e^{Qu}\right)\right\rVert_{1}+\left\lVert\widetilde{p}_{t}-p_{t}\right\rVert_{1}=\left\lVert\widetilde{p}_{t}\left(e^{\widetilde{Q}u}-e^{Qu}\right)\right\rVert_{1}+\left\lVert
e_{t}\right\rVert_{1}$ $\displaystyle\implies\left\lVert
e_{t+u}\right\rVert_{1}-\left\lVert e_{t}\right\rVert_{1}$
$\displaystyle\leq\left\lVert\widetilde{p}_{t}\left(e^{\widetilde{Q}u}-e^{Qu}\right)\right\rVert_{1}=\sum_{s\in
S}\left|\left(\widetilde{p}_{t}\left(e^{\widetilde{Q}u}-e^{Qu}\right)\right)(s)\right|$
(3.7)

where $\circledcirc$ follows from 3 since $e^{Qu}$ is a stochastic matrix. We
now take a closer look at the right hand side. In particular, we are
interested in

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}u}\widetilde{p}_{t}\left(e^{\widetilde{Q}u}-e^{Qu}\right)$
$\displaystyle=\widetilde{p}_{t}\left(e^{\widetilde{Q}u}\widetilde{Q}-e^{Qu}Q\right)$
$\displaystyle=\widetilde{p}_{t}\left(\widetilde{Q}-Q\right)\textrm{ for }u=0$
$\displaystyle\implies\left.\frac{\mathrm{d}}{\mathrm{d}u}\right|_{u=0}\left(\widetilde{p}_{t}\left(e^{\widetilde{Q}u}-e^{Qu}\right)\right)(s)$
$\displaystyle=\left(\widetilde{p}_{t}\left(\widetilde{Q}-Q\right)\right)(s)=\sum_{r\in
S}\widetilde{p}_{t}(r)\left(\widetilde{Q}(r,s)-Q(r,s)\right)$

Hence

$\displaystyle\limsup_{u\to 0}\frac{\left\lVert
e_{t+u}\right\rVert_{1}-\left\lVert e_{t}\right\rVert_{1}}{u}$
$\displaystyle\overset{\textrm{\eqref{eq:errdiff_bound}}}{\leq}\limsup_{u\to
0}\sum_{s\in
S}\frac{\left|\left(\widetilde{p}_{t}\left(e^{\widetilde{Q}u}-e^{Qu}\right)\right)(s)\right|}{u}$
$\displaystyle\leq\sum_{s\in S}\limsup_{u\to
0}\frac{\left|\left(\widetilde{p}_{t}\left(e^{\widetilde{Q}u}-e^{Qu}\right)\right)(s)\right|}{u}$
$\displaystyle\hskip-11.38092pt\overset{\textrm{\autoref{lem:absfctlim}}}{=}\sum_{s\in
S}\left|\sum_{r\in
S}\widetilde{p}_{t}(r)\left(\widetilde{Q}(r,s)-Q(r,s)\right)\right|\overset{\circledast}{\leq}\sum_{\rho\in\Omega}\pi_{t}(\rho)\cdot\tau(\rho)=\left<\pi_{t},\tau\right>$

where $\circledast$ follows from the exact same calculation as in the case of
discrete time, see (3.2).

Looking at the equation

$\displaystyle\left\lVert e_{t}\right\rVert_{1}$ $\displaystyle\leq\left\lVert
e_{0}\right\rVert_{1}+\int_{0}^{t}\left<\pi_{s},\tau\right>\;\mathrm{d}s$

one can see that, as was the case for discrete time, the $\tau$ factors
provide a bound for the “dynamic” error, resulting from using approximated
transition rates between states, which accumulates over time, and which is
separate from the “static” error occurring when approximating $p_{0}$.
Therefore, low $\tau$ factors only provide an indication that the error
accumulated over time is low, but it can still be that the error in the
initial distribution $\widetilde{p}_{0}$ is high.

###  When is the error bound 0?

Later, it is useful to further subdivide the error factors $\tau$ into error
factors for single transitions. We will therefore use the notation

$\displaystyle\tau(\rho,\sigma)$
$\displaystyle:=\sum_{s\in\sigma}\left|\alpha(s)\Pi(\rho,\sigma)-\sum_{r\in\rho}\alpha(r)P(r,s)\right|$
$\displaystyle\textrm{or, for continuous time: }\quad\tau(\rho,\sigma)$
$\displaystyle:=\sum_{s\in\sigma}\left|\alpha(s)\Theta(\rho,\sigma)-\sum_{r\in\rho}\alpha(r)Q(r,s)\right|$

In order to better understand the error bound which arises from the $\tau$
factors, this section analyses in which cases the error bound is equal to $0$.

###### Lemma 6

Given a DTMC, a partition $\Omega$ of its state space, arbitrary distributions
$\alpha_{\sigma}$ with support on $\sigma\in\Omega$, and arbitrary $\Pi$, it
holds that

$\displaystyle\forall\rho\in\Omega:\tau(\rho)=0\iff\Pi
A=AP\;\mathrel{\mathop{\iff}\limits^{\vbox
to5.69054pt{\kern-2.0pt\hbox{$\scriptstyle\textrm{\autoref{def:exactagg}}$}\vss}}}\;\textrm{the
aggregation is dynamic-exact}$

The same statement holds for CTMCs with $\Pi$ replaced by $\Theta$ and $P$
replaced by $Q$.

###### Proof

Note that $\tau(\rho)=0$ if, and only if, $\tau(\rho,\sigma)=0$ for all
$\sigma$. Hence, we have that

$\displaystyle\begin{split}\forall\rho\in\Omega:\tau(\rho)=0&\iff\forall\rho,\sigma\in\Omega:\tau(\rho,\sigma)=0\\\
&\iff\forall\rho,\sigma\in\Omega:\sum_{s\in\sigma}\left|\alpha(s)\Pi(\rho,\sigma)-\sum_{r\in\rho}\alpha(r)P(r,s)\right|=0\\\
&\iff\forall\rho,\sigma\in\Omega:\forall
s\in\sigma:\;\;\alpha(s)\Pi(\rho,\sigma)=\sum_{r\in\rho}\alpha(r)P(r,s)\end{split}$
(3.8)

This already proves 6 since we have, with $\omega(s)=\sigma$:

$\displaystyle\alpha(s)\Pi(\rho,\sigma)=\underbrace{(\Pi
A)}_{\in\mathbb{R}^{m\times n}}(\rho,s)\qquad\textrm{ and
}\qquad\sum_{r\in\rho}\alpha(r)P(r,s)=(AP)(\rho,s)$

The same calculation is true with $\Pi$ replaced by $\Theta$ and $P$ replaced
by $Q$ (also compare with (3.4), where the same was basically already proven).

###### Corollary 7

Given a DTMC or CTMC, a partition $\Omega$ of its state space, arbitrary
distributions $\alpha_{\sigma}$ with support on $\sigma\in\Omega$, and
arbitrary $\Pi$ (respectively $\Theta$) such that $\Pi A=AP$ (respectively
$\Theta A=AQ$, i.e. a dynamic-exact aggregation), it holds that

$\displaystyle\left\lVert\widetilde{p}_{k}-p_{k}\right\rVert_{1}\leq\left\lVert\widetilde{p}_{0}-p_{0}\right\rVert_{1}\textrm{
or, for continuous time,
}\left\lVert\widetilde{p}_{t}-p_{t}\right\rVert_{1}\leq\left\lVert\widetilde{p}_{0}-p_{0}\right\rVert_{1}$

In particular, if an aggregation is exact, then
$\left\lVert\widetilde{p}_{k}-p_{k}\right\rVert_{1}=0$ for all $k$
(respectively $\left\lVert\widetilde{p}_{t}-p_{t}\right\rVert_{1}=0$ for all
$t$).

###### Proof

This follows from 6 and (3.3) (respectively Theorem 4 for CTMCs).

We next compare the notion of exactness and dynamic-exactness as defined in 2
with a similar notion from [9]. In particular, [9, Theorem 2] and [9, Theorem
9] state that the distributions $\widetilde{p}_{k}$ (respectively
$\widetilde{p}_{t}$) are equal to $p_{k}$ (respectively $p_{t}$) under the
following conditions:

  1. (i)

$\displaystyle\forall s,s^{\prime}\in S\textrm{ s.t.
}\omega(s)=\omega(s^{\prime}):\forall\rho\in\Omega:\frac{\sum_{r\in\rho}\alpha(r)P(r,s)}{\alpha(s)}=\frac{\sum_{r\in\rho}\alpha(r)P(r,s^{\prime})}{\alpha(s^{\prime})}$

  2. (ii)

the initial distribution $p_{0}$ is compatible with the $\alpha$ distributions
(this is called $p_{0}$ respects the $\alpha$ distributions in [9])

The next proposition, 8, shows that i is equivalent to
$\forall\rho\in\Omega:\tau(\rho)=0$ and thus by 6 equivalent to dynamic-
exactness. It follows that conditions i and ii together are equivalent to
exactness, and hence imply (see 7) that the transient distributions
$\widetilde{p}_{k}$ (respectively $\widetilde{p}_{t}$) agree with $p_{k}$
(respectively $p_{t}$), which is exactly the statement of [9, Theorem 2] and
[9, Theorem 9], but obtained via a different proof.

###### Proposition 8

Given an irreducible DTMC, a partition $\Omega$ of its state space, and
arbitrary distributions $\alpha_{\sigma}$ with support on $\sigma\in\Omega$,
let $\Pi$ as in (2.2). Then:

$\displaystyle\forall\rho\in\Omega:\tau(\rho)=0$ $\displaystyle\iff$
$\displaystyle\forall s\in S:\alpha(s)>0\textrm{ and }$ $\displaystyle\forall
s,s^{\prime}\in S\textrm{ s.t.
}\omega(s)=\omega(s^{\prime}):\forall\rho\in\Omega:\frac{\sum_{r\in\rho}\alpha(r)P(r,s)}{\alpha(s)}=\frac{\sum_{r\in\rho}\alpha(r)P(r,s^{\prime})}{\alpha(s^{\prime})}$

The same statement holds for irreducible CTMCs with $\Pi$ replaced by $\Theta$
and $P$ replaced by $Q$.

###### Proof

We have that (compare with (3.8))

$\displaystyle\forall\rho\in\Omega:\tau(\rho)=0$
$\displaystyle\iff\forall\rho,\sigma\in\Omega:\forall
s\in\sigma:\;\;\alpha(s)\Pi(\rho,\sigma)=\sum_{r\in\rho}\alpha(r)P(r,s)$ (3.9)

We first show: $\forall\rho\in\Omega:\tau(\rho)=0$ implies $\forall s\in
S:\alpha(s)>0$. Call $S_{0}:=\\{s_{0}\in S:\alpha(s_{0})=0\\}\neq S$ (we have
$S_{0}\neq S$ because the $\alpha_{\sigma}$ are probability distributions, so
there must be at least one state $s_{\neg 0}$ with $\alpha(s_{\neg 0})>0$). If
$S_{0}\neq\varnothing$, then, by our assumption of irreducibility, there must
be some state $s\in S_{0}$ and some state $r\in S\setminus S_{0}$ (hence
$\alpha(r)>0$) with $P(r,s)>0$. In particular, we have that
$\sum_{r^{\prime}\in\omega(r)}\alpha(r^{\prime})P(r^{\prime},s)>0$ while
$\alpha(s)\Pi(\omega(r),\omega(s))=0$. Thus, since the last statement of (3.9)
does not hold in this case (it is violated for $\rho=\omega(r)$,
$\sigma=\omega(s)$ and $s$), it follows that there is some $\rho\in\Omega$
such that $\tau(\rho)\neq 0$. $S_{0}\neq\varnothing$ is thus impossible when
all $\tau$ factors are zero, and hence, writing $\circledast$ for $\forall
s\in S:\alpha(s)>0$:

$\displaystyle\forall\rho\in\Omega:\tau(\rho)=0$
$\displaystyle\iff\circledast\textrm{ and }\forall\rho,\sigma\in\Omega:\forall
s\in\sigma:\;\;\alpha(s)\Pi(\rho,\sigma)=\sum_{r\in\rho}\alpha(r)P(r,s)$
$\displaystyle\iff\circledast\textrm{ and }\forall\rho,\sigma\in\Omega:\forall
s\in\sigma:\;\;\Pi(\rho,\sigma)=\frac{\sum_{r\in\rho}\alpha(r)P(r,s)}{\alpha(s)}$
$\displaystyle\overset{\circledcirc}{\iff}\circledast\textrm{ and
}\forall\rho,\sigma\in\Omega:\forall
s,s^{\prime}\in\sigma:\;\;\frac{\sum_{r\in\rho}\alpha(r)P(r,s)}{\alpha(s)}=\frac{\sum_{r\in\rho}\alpha(r)P(r,s^{\prime})}{\alpha(s^{\prime})}$

We show $\circledcirc$:

  * •

$\implies$: this is immediately clear since $\Pi(\rho,\sigma)$ does not depend
on $s$.

  * •

$\impliedby$: if we have

$\displaystyle\frac{\sum_{r\in\rho}\alpha(r)P(r,s)}{\alpha(s)}=\frac{\sum_{r\in\rho}\alpha(r)P(r,s^{\prime})}{\alpha(s^{\prime})}\;\;\forall
s,s^{\prime}\in\sigma$

then, it holds that, for any $s\in\sigma$,

$\displaystyle\frac{\sum_{r\in\rho}\alpha(r)P(r,s)}{\alpha(s)}=\sum_{s^{\prime}\in\sigma}\alpha(s^{\prime})\cdot\frac{\sum_{r\in\rho}\alpha(r)P(r,s^{\prime})}{\alpha(s^{\prime})}=\sum_{r\in\rho}\alpha(r)\sum_{s^{\prime}\in\sigma}P(r,s^{\prime})=\Pi(\rho,\sigma)$

where the first equality holds since the sum on the right hand side consists
of a weighted sum of terms which are equal, and since the weights sum up to
$1$.

Note that the above proof also works for the continuous-time case.

###### Remark

If we want to state 8 in words, we could say: if we set $\Pi$ as in (2.2),
then the overall error bound is $0$ for DTMCs if, and only if,

$\displaystyle\frac{\alpha(s)}{\alpha(s^{\prime})}=\frac{\sum_{r\in\rho}\alpha(r)P(r,s)}{\sum_{r\in\rho}\alpha(r)P(r,s^{\prime})}$

for any two states $s$ and $s^{\prime}$ in the same aggregate $\sigma$ and any
aggregate $\rho$ such that the probability to go from $\rho$ to $s^{\prime}$
is positive.

Note that for CTMCs, $\sum_{r\in\rho}\alpha(r)Q(r,s)$ can be negative if
$s\in\rho$. However, if, in addition, $s^{\prime}\in\rho$ and if
$\tau(\sigma)=0$ for all $\sigma$, then
$\sum_{r\in\rho}\alpha(r)Q(r,s^{\prime})$ must be negative as well by 8,
because $\alpha(s)>0$ and $\alpha(s^{\prime})>0$. Usually,
$\sum_{r\in\rho}\alpha(r)Q(r,s)$ can be interpreted as the rate at which
probability mass is flowing from $\rho$ to $s$, under the assumption that
$s\notin\rho$ and that the mass within $\rho$ is distributed according to
$\alpha$. Now, if $s\in\rho$, this is no longer the case. Instead,

$\displaystyle\sum_{r\in\rho}\alpha(r)Q(r,s)=\underbrace{\alpha(s)Q(s,s)}_{\textrm{negative
of the rate at which mass is exiting
}s}+\underbrace{\sum_{\begin{subarray}{c}r\in\rho\\\ r\neq
s\end{subarray}}\alpha(r)Q(r,s)}_{\textrm{rate at which mass is flowing from
}\rho\setminus\\{s\\}\textrm{ to }s}$

again under the assumption that the mass within $\rho$ is distributed
according to $\alpha$.

With the help of 8, we can now apply other results from [9] (we assume that
$\Pi$, respectively $\Theta$, is set as in (2.2)):

  * •

By [9, Theorem 7], if $\tau(\rho)=0$ for all $\rho$ and if the DTMC is
aperiodic and irreducible, then we have
$\left\lVert\pi_{k}-p_{k}\Lambda\right\rVert_{1}\to 0$ for $k\to\infty$ (i.e.
the approximate aggregate probabilities $\pi_{k}(\sigma)$ converge to the
exact aggregate probability $\mathbb{P}\left[X_{k}\in\sigma\right]$) and
$\frac{p_{k}(s)}{\pi_{k}(\omega(s))}\to\alpha(s)$. For periodic chains, we
have to consider $\frac{1}{k}\sum_{i=1}^{k}p_{i}(s)$ instead of $p_{k}(s)$,
see [9, Theorem 6]. In addition, the stationary distribution $\mu$ of the
irreducible DTMC satisfies $\mu=\widehat{\mu}A$ where $\widehat{\mu}$ is the
stationary distribution of the aggregated chain in both periodic and aperiodic
cases.

  * •

By [9, Theorem 11], if $\tau(\rho)=0$ for all $\rho$ and if the CTMC is
irreducible, then we have $\left\lVert\pi_{t}-p_{t}\Lambda\right\rVert_{1}\to
0$ for $t\to\infty$ and $\frac{p_{t}(s)}{\pi_{t}(\omega(s))}\to\alpha(s)$. In
addition, the stationary distribution $\mu$ of the CTMC satisfies
$\mu=\widehat{\mu}A$ where $\widehat{\mu}$ is the stationary distribution of
the aggregated chain.

We have seen that $\tau(\rho)=0$ for all $\rho$ entails a number of desirable
properties of the aggregation, and that other sources (see, e.g., [9] or [2,
Definition 2.1], where the notion of aggregatability is defined, a concept
which implies dynamic-exactness) have identified similar, sometimes
equivalent, conditions under which aggregation without loss of precision is
possible.

###  Tightness of the error bound

We now want to show that the bounds given by the $\tau$ factors are also
tight, in the sense that no better bounds are possible in general, meaning
that $\tau(\rho)>0$ for some $\rho$ implies that the $\widetilde{p}_{k}$ are
not always equal to $p_{k}$.

###### Theorem 9

Given a DTMC or CTMC, a partition $\Omega$, arbitrary distributions
$\alpha_{\sigma}$ with support on $\sigma\in\Omega$, and arbitrary $\Pi$ or
$\Theta$, assume that $\tau(\rho)>0$ for some $\rho$. Then, there exists an
initial distribution $p_{0}$ which is compatible with the $\alpha$
distributions such that
$\left\lVert\widetilde{p}_{1}-p_{1}\right\rVert_{1}=\left<\pi_{0},\tau\right>$
or, for CTMCs, $\lim_{t\to
0,t>0}\frac{1}{t}\left\lVert\widetilde{p}_{t}-p_{t}\right\rVert_{1}=\left<\pi_{0},\tau\right>$.

###### Proof

There must be some aggregate $\rho\in\Omega$ with $\tau(\rho)>0$. We choose
$p_{0}=\alpha_{\rho}$, which is clearly compatible with the $\alpha$
distributions. In the discrete-time case, we have $p_{1}=\alpha_{\rho}P$ and
$\widetilde{p}_{1}=\alpha_{\rho}\Lambda\Pi A$. Hence

$\displaystyle\left\lVert\widetilde{p}_{1}-p_{1}\right\rVert_{1}$
$\displaystyle=\sum_{s\in
S}\left|\widetilde{p}_{1}(s)-p_{1}(s)\right|=\sum_{s\in
S}\left|\alpha(s)\Pi(\rho,\omega(s))-\sum_{r\in\rho}\alpha(r)P(r,s)\right|$
$\displaystyle=\sum_{\sigma\in\Omega}\sum_{s\in\sigma}\left|\alpha(s)\Pi(\rho,\sigma)-\sum_{r\in\rho}\alpha(r)P(r,s)\right|=\tau(\rho)=\left<\pi_{0},\tau\right>$

since $\pi_{0}$ is the Dirac measure on $\rho$. We proceed simlarly for the
continuous-time case, again setting $p_{0}=\alpha_{\rho}$. Note that, as
already established in (3.6), we have

$\displaystyle\left.\frac{\textrm{d}}{\textrm{d}t}\right|_{t=0}(\widetilde{p}_{t}-p_{t})(s)$
$\displaystyle=\left(\widetilde{p}_{0}\widetilde{Q}-p_{0}Q\right)(s)\;\;\overset{\widetilde{p}_{0}=p_{0}=\alpha_{\rho}}{=}\;\;\left(\alpha_{\rho}\widetilde{Q}-\alpha_{\rho}Q\right)(s)$
$\displaystyle=\left(\alpha_{\rho}\Lambda\Theta A-\alpha_{\rho}Q\right)(s)$

Hence, noting that $\left\lVert\widetilde{p}_{0}-p_{0}\right\rVert_{1}=0$, we
obtain

$\displaystyle\lim_{\begin{subarray}{c}t\to 0\\\
t>0\end{subarray}}\frac{1}{t}\left\lVert\widetilde{p}_{t}-p_{t}\right\rVert_{1}$
$\displaystyle=\sum_{s\in S}\lim_{\begin{subarray}{c}t\to 0\\\
t>0\end{subarray}}\frac{\left|(\widetilde{p}_{t}-p_{t})(s)\right|}{t}\;\;\overset{\textrm{\autoref{lem:absfctlim}}}{=}\;\;\sum_{s\in
S}\left|\left(\alpha_{\rho}\Lambda\Theta A-\alpha_{\rho}Q\right)(s)\right|$
$\displaystyle=\ldots\textrm{ (as in the discrete-time case)
}=\left<\pi_{0},\tau\right>$

###### Remark

If $\tau(\rho)>0$ for some $\rho$, it can still be true that the stationary
distribution $\mu$ of the Markov chain satisfies $\mu=\widehat{\mu}A$ where
$\widehat{\mu}$ is the stationary distribution of the aggregated chain (see
next paragraph). The $\tau$ factors can therefore not be an optimal bound for
$\left\lVert\mu-\widehat{\mu}A\right\rVert_{1}$, i.e. a bound for how far off
the approximated stationary distribution, obtained via the aggregated chain,
is from the actual stationary distribution of the chain.

Indeed, if we take an arbitrary DTMC with stationary distribution $\mu$ and an
arbitrary partition $\Omega$, then we can choose the distributions $\alpha$
such that they are compatible with $\mu$ by simply setting
$\alpha(s)=\mu(s)(\sum_{r\in\omega(s)}\mu(r))^{-1}$. But then, setting
$\widehat{\mu}:=\mu\Lambda$ (which implies $\widehat{\mu}A=\mu$ by the
compatibility), we obtain

$\displaystyle\left(\widehat{\mu}\Pi\right)(\sigma)$
$\displaystyle=\sum_{\rho\in\Omega}\widehat{\mu}(\rho)\Pi(\rho,\sigma)=\sum_{\rho\in\Omega}\widehat{\mu}(\rho)\sum_{r\in\rho}\alpha(r)\sum_{s\in\sigma}P(r,s)$
$\displaystyle=\sum_{\rho\in\Omega}\sum_{r\in\rho}\underbrace{\alpha(r)\widehat{\mu}(\rho)}_{=\mu(r)}\sum_{s\in\sigma}P(r,s)=\sum_{s\in\sigma}\underbrace{\sum_{\rho\in\Omega}\sum_{r\in\rho}\mu(r)P(r,s)}_{\mu(s)}=\sum_{s\in\sigma}\mu(s)=\widehat{\mu}(\sigma)$

i.e. $\widehat{\mu}$ is indeed a stationary measure for the aggregated chain.
But since our choice of DTMC and partition was arbitrary, it is of course not
necessarily the case that $\tau(\rho)=0$ for all $\rho$ (take, for example,
the partition $\Omega=\\{S\\}$).

###  The median-based scheme and the error bounds

In this section, we give a proof of the statement that the median-based scheme
minimizes the error bounds for fixed partition $\Omega$ and fixed $\alpha$
distributions.

###### Proposition 10

Assume that the partition $\Omega$ is fixed, and that the distributions
$\alpha$ are fixed as well (but arbitrary). Then, setting $\Pi$ (or $\Theta$)
according to (2.3) and (2.4) minimizes $\tau(\rho)$ for every $\rho\in\Omega$
among all matrices $\Pi\in\mathbb{R}^{m\times m}$ (or
$\Theta\in\mathbb{R}^{m\times m}$).

###### Proof

Recall the definition of $\tau(\rho)$:

$\displaystyle\tau(\rho)=\sum_{\sigma\in\Omega}\underbrace{\sum_{s\in\sigma}\left|\alpha(s)\Pi(\rho,\sigma)-\sum_{r\in\rho}\alpha(r)P(r,s)\right|}_{=\tau(\rho,\sigma)}=\sum_{\sigma\in\Omega}\sum_{s\in\sigma}\alpha(s)\cdot\left|\Pi(\rho,\sigma)-\frac{1}{\alpha(s)}\sum_{r\in\rho}\alpha(r)P(r,s)\right|$

Minimizing $\tau(\rho)$ (with $\alpha$ and $\Omega$ fixed) amounts to
minimizing $\tau(\rho,\sigma)$ separately for each $\sigma\in\Omega$, since we
can choose a different value for $\Pi(\rho,\sigma)$ for every $\sigma$. We
thus want to set $\Pi(\rho,\sigma)$ by solving the following optimization
problem:

$\displaystyle\Pi(\rho,\sigma):=\operatorname*{arg\,min}_{x\in\mathbb{R}}\sum_{s\in\sigma}\alpha(s)\cdot\left|x-\vphantom{\frac{1}{\alpha(s)}\sum_{r\in\rho}}\right.\underbrace{\frac{1}{\alpha(s)}\sum_{r\in\rho}\alpha(r)P(r,s)}_{a_{s}}\left.\vphantom{\frac{1}{\alpha(s)}\sum_{r\in\rho}}\right|$

Note: it is easy to see that we can assume $\alpha(s)>0$ for all $s\in S$
w.l.o.g. in this proof. We show that the weighted median of the (multi-)set
$\\{a_{s}\;|\;s\in\sigma\\}$ with weights $\alpha(s),\;s\in\sigma$ does indeed
solve the given minimization problem. For ease of notation and w.l.o.g., we
write $\\{a_{s}\;|\;s\in\sigma\\}=\\{a_{1},a_{2},\ldots,a_{k}\\}$, we assume
$a_{1}<a_{2}<\ldots<a_{k}$ (we can merge elements $a_{i}$ and $a_{i+1}$ if
$a_{i}=a_{i+1}$ into a new element whose weight is the sum of the two weights)
and we denote the corresponding weights by $\alpha(1),\ldots,\alpha(k)$.
Consider the function
$f:\mathbb{R}\to\mathbb{R},f(x)=\sum_{i=1}^{k}\alpha(i)|x-a_{i}|$. Note that
this function is piecewise linear: for $x\in(a_{i},a_{i+1})$, $f$ is a linear
combination of linear functions and thus linear. On the other hand, when
$x=a_{i}$ for some $i$, the slope of $f$ might change.

By this piecewise linearity and since $\lim_{x\to\pm\infty}f(x)=\infty$, we
know that $f$ must take a global minimum for some $x$ in the set
$\\{a_{1},\ldots,a_{k}\\}$. Denote the value of $x$ for which the minimum is
taken by $a_{i^{\ast}}$. Assume for a contradiction that $a_{i^{\ast}}$ is not
the weighted median. W.l.o.g. we consider the case where
$\sum_{i>i^{\ast}}\alpha(i)>\frac{1}{2}$, i.e. the case where the total weight
of the elements which are bigger than $a_{i^{\ast}}$ is bigger than
$\frac{1}{2}$. Note that we also assumed above (w.l.o.g.) that no two elements
in the set $\\{a_{s}\;|\;s\in\sigma\\}$ are identical, i.e. we have
$a_{i^{\ast}+1}>a_{i^{\ast}}$.

We proceed to show that $f(a_{i^{\ast}+1})<f(a_{i^{\ast}})$, and hence the
minimum is not taken at $a_{i^{\ast}}$ and our assumption that $a_{i^{\ast}}$
is not the weighted median must have been wrong. To see that this is true,
consider the following:

$\displaystyle f(a_{i^{\ast}+1})$
$\displaystyle=\sum_{i=1}^{k}\alpha(i)|a_{i^{\ast}+1}-a_{i}|$
$\displaystyle=\sum_{i=1}^{i^{\ast}}\alpha(i)|a_{i^{\ast}+1}-a_{i}|+\underbrace{\alpha(i^{\ast}+1)|a_{i^{\ast}+1}-a_{i^{\ast}+1}|}_{0}+\sum_{i=i^{\ast}+2}^{k}\alpha(i)|a_{i^{\ast}+1}-a_{i}|$
$\displaystyle=\sum_{i=1}^{i^{\ast}}\alpha(i)\left(|a_{i^{\ast}}-a_{i}|+|a_{i^{\ast}+1}-a_{i^{\ast}}|\right)+\sum_{i=i^{\ast}+2}^{k}\alpha(i)\left(|a_{i^{\ast}}-a_{i}|-|a_{i^{\ast}+1}-a_{i^{\ast}}|\right)$
$\displaystyle=\underbrace{\sum_{i=1}^{k}\alpha(i)|a_{i^{\ast}}-a_{i}|}_{f(a_{i^{\ast}})}+|a_{i^{\ast}+1}-a_{i^{\ast}}|\underbrace{\left(\vphantom{\sum_{i=1}^{i^{\ast}}}\right.\underbrace{\sum_{i=1}^{i^{\ast}}\alpha(i)}_{<\frac{1}{2}}-\underbrace{\sum_{i=i^{\ast}+1}^{k}\alpha(i)}_{>\frac{1}{2}}\left.\vphantom{\sum_{i=1}^{i^{\ast}}}\right)}_{<0}$

Note that the additional term (if compared to the previous line) in the
leftmost sum on the last line is compensated by the additional term in the
rightmost sum.

The proof above also works for the continuous-time case.

Summarizing this section, we have introduced the $\tau$ factors which give a
formal error bound on $\left\lVert\widetilde{p}_{k}-p_{k}\right\rVert_{1}$ or
$\left\lVert\widetilde{p}_{t}-p_{t}\right\rVert_{1}$. In general, these $\tau$
factors provide the best possible error bounds for the transient
distributions. $\tau(\rho)=0$ for all $\rho$ is equivalent to dynamic-
exactness which implies that the error
$\left\lVert\widetilde{p}_{k}-p_{k}\right\rVert_{1}$ or
$\left\lVert\widetilde{p}_{t}-p_{t}\right\rVert_{1}$ is bounded by the initial
error as well as that the exact stationary distribution $\mu$ of the chain can
be derived from the stationary distribution $\widehat{\mu}$ of the aggregated
chain via $\mu=\widehat{\mu}A$.

##  Lumpability and aggregatability

In this section, we will consider how various definitions of lumpability and
aggregatability from the literature affect the error bounds $\tau$. This
analysis also provides insights into why the proportional and uniform choices
for the distributions $\alpha$ are sensible.

The following definition was given in [3, Definition 1]:

###### Definition 11

A partition $\Omega=\\{\Omega_{1},\ldots,\Omega_{m}\\}$ of the state space of
a DTMC is called ordinarily lumpable if

$\displaystyle\forall r,r^{\prime}\in S\textrm{ s.t.\leavevmode\nobreak\
}\omega(r)=\omega(r^{\prime}):\forall\sigma\in\Omega:\qquad\sum_{s\in\sigma}P(r,s)=\sum_{s\in\sigma}P(r^{\prime},s)$
(4.1)

That is, for any two states in the same aggregate, the outgoing probabilities
to any other aggregate must be identical. For CTMCs, a partition is called
ordinarily lumpable if (4.1) holds with $P$ replaced by $Q$.

For an ordinarily lumpable partition, it is easy to show that
$\pi_{k}(\sigma)=\sum_{s\in\sigma}\widetilde{p}_{k}(s)=\sum_{s\in\sigma}p_{k}(s)$
if $\Pi$ is set as in (2.2), for any $\sigma\in\Omega$, all $k$, and any
initial distribution $p_{0}$ (and independently of the choice of $\alpha$),
i.e. the probability $\mathbb{P}\left[X_{k}\in\sigma\right]$ is exactly equal
to $\pi_{k}(\sigma)$. A simple induction suffices to prove this, see [3,
Theorem 5]. The same holds in the continuous-time case.

Note that ordinary lumpability is called strong lumpability in [11], which
considers only DTMCs. As before, define the process $Y_{k}$ by
$Y_{k}=\sigma\in\Omega\iff X_{k}\in\Omega$. By [11, Theorem 6.3.2], ordinary
lumpability is equivalent to the following: for every initial distribution
$p_{0}$, $Y_{k}$ is a Markov chain (whose transition probabilities do not
depend on the choice of $p_{0}$).

[3, Definition 1] also defines exact lumpability:

###### Definition 12

A partition $\Omega=\\{\Omega_{1},\ldots,\Omega_{m}\\}$ of the state space of
a DTMC is called exactly lumpable if

$\displaystyle\forall s,s^{\prime}\in S\textrm{ s.t.\leavevmode\nobreak\
}\omega(s)=\omega(s^{\prime}):\forall\rho\in\Omega:\qquad\sum_{r\in\rho}P(r,s)=\sum_{r\in\rho}P(r,s^{\prime})$
(4.2)

That is, for any two states in the same aggregate, the incoming probabilities
from any other aggregate must be identical. For CTMCs, a partition is called
exactly lumpable if (4.2) holds with $P$ replaced by $Q$.

Furthermore, a partition $\Omega$ is called strictly lumpable if it is both
ordinarily and exactly lumpable.

If a partition is ordinarily (exactly) lumpable for a CTMC, then the partition
is also ordinarily (exactly) lumpable for any uniformisation of the CTMC,
regardless of the uniformisation rate, see the remarks after [3, Definition
1].

[2, Definition 2.1] also defines lumpability. Note that this definition of
lumpability agrees with the definition of ordinary lumpability given above.
[2, Definition 2.1] further defines deflatability and aggregatability:

###### Definition 13

A partition $\Omega=\\{\Omega_{1},\ldots,\Omega_{m}\\}$ of the state space of
a DTMC, together with distributions $\alpha_{\sigma}\in\mathbb{R}^{n}$ with
support on $\sigma\in\Omega$, is called deflatable if

$\displaystyle\forall r\in S:\forall s\in S:\qquad
P(r,s)=\alpha(s)\cdot\sum_{s^{\prime}\in\omega(s)}P(r,s^{\prime})$ (4.3)

In words: the probability to go from $r$ to $s$ only depends on $r$ and the
aggregated state $\omega(s)$ as well as a factor which depends on $s$, but not
on $r$. Another description: after a jump into a partition element $\sigma$,
the particular target state can be chosen according to the probability measure
$\alpha_{\sigma}$, independently of where the jump started.

The partition $\Omega$, together with distributions $\alpha$, is further
called aggregatable if it is deflatable and if $\Omega$ is ordinarily
lumpable.

Note that aggregatability implies that $\widetilde{P}=P$ by [2, Proposition
2.6]. 13 cannot be extended to CTMCs easily. It therefore might make more
sense to consider dynamic-exact aggregations instead of aggregatable
partitions, i.e. a partition together with distributions $\alpha_{\sigma}$
such that $\Pi A=AP$ (or $\Theta A=AQ$), motivated by 2 and the following
remarks. Dynamic-exactness can be defined for both DTMCs and CTMCs.

###### Proposition 14

Assume an irreducible DTMC or CTMC is exactly lumpable w.r.t. the partition
$\Omega$. When setting $\Pi$ (respectively $\Theta$) as in (2.2) and $\alpha$
as in (2.5) (respectively (2.7) (proportional $\alpha$)) _or_ as in (2.6)
(uniform $\alpha$), then it holds that

$\displaystyle\tau(\rho)$
$\displaystyle=0\;\;\forall\rho\in\Omega\qquad\textrm{and}\qquad\alpha(s)=\frac{1}{\left|\omega(s)\right|}\;\;\forall
s\in S$

###### Proof

We first show that $\alpha(s)=\frac{1}{\left|\omega(s)\right|}$ if the Markov
chain is exactly lumpable. Let $s,s^{\prime}\in S$ such that
$\omega(s)=\omega(s^{\prime})=\sigma$ and assume exact lumpability. First, we
look at DTMCs and assume that $\alpha$ is set as in (2.5). Then:

$\displaystyle\alpha(s^{\prime})\cdot\sum_{r\in
S}\sum_{s^{\prime\prime}\in\sigma}P(r,s^{\prime\prime})$
$\displaystyle=\sum_{r\in
S}P(r,s^{\prime})=\sum_{\rho\in\Omega}\sum_{r\in\rho}P(r,s^{\prime})\;\overset{\textrm{ex.\leavevmode\nobreak\
lump.}}{=}\;\sum_{\rho\in\Omega}\sum_{r\in\rho}P(r,s)$
$\displaystyle=\sum_{r\in S}P(r,s)=\alpha(s)\cdot\underbrace{\sum_{r\in
S}\sum_{s^{\prime\prime}\in\sigma}P(r,s^{\prime\prime})}_{>0\textrm{ by
irreducibility}}$

This implies $\alpha(s)=\alpha(s^{\prime})$ and hence, $\alpha_{\sigma}$ must
be the uniform distribution on $\sigma$. For CTMCs and $\alpha$ set as in
(2.7), an almost identical calculation but summing over $r\in
S\setminus\sigma$ instead of $r\in S$ yields the desired result (which would
not be the case if we only exclude the diagonal elements of $Q$ in (2.7)
instead of the “diagonal blocks”).

We now want to show $\tau(\rho)=0$ for all $\rho$. By 8, this is equivalent to
(in the following, we use that $\alpha(s)=\frac{1}{\left|\omega(s)\right|}$)

$\displaystyle\begin{split}&\forall s,s^{\prime}\in S\textrm{
s.t.\leavevmode\nobreak\
}\omega(s)=\omega(s^{\prime}):\forall\rho\in\Omega:\alpha(s^{\prime})\sum_{r\in\rho}\alpha(r)P(r,s)=\alpha(s)\sum_{r\in\rho}\alpha(r)P(r,s^{\prime})\\\
\hskip-5.0pt\iff\;&\forall s,s^{\prime}\in S\textrm{ s.t.\leavevmode\nobreak\
}\omega(s)=\omega(s^{\prime}):\forall\rho\in\Omega:\frac{1}{\left|\omega(s^{\prime})\right|}\sum_{r\in\rho}\frac{1}{\left|\rho\right|}P(r,s)=\frac{1}{\left|\omega(s)\right|}\sum_{r\in\rho}\frac{1}{\left|\rho\right|}P(r,s^{\prime})\\\
\hskip-5.0pt\iff\;&\forall s,s^{\prime}\in S\textrm{ s.t.\leavevmode\nobreak\
}\omega(s)=\omega(s^{\prime}):\forall\rho\in\Omega:\sum_{r\in\rho}P(r,s)=\sum_{r\in\rho}P(r,s^{\prime})\\\
\hskip-5.0pt\iff\;&\Omega\textrm{ is exactly lumpable}\end{split}$ (4.4)

The same calculation holds for CTMCs.

###### Proposition 15

Assume a partition $\Omega$ is given for an irreducible DTMC or CTMC. When
setting $\Pi$ (respectively $\Theta$) as in (2.2) and $\alpha$ as in (2.6)
(uniform $\alpha$), then it holds that

$\displaystyle\forall\rho\in\Omega:\tau(\rho)$
$\displaystyle=0\;\;\iff\;\;\Omega\textrm{ is exactly lumpable}$

###### Proof

See (4.4) and surrounding remarks.

###### Proposition 16

Given an irreducible DTMC and a partition $\Omega$ of its state space, assume
that

$\displaystyle\forall r\in S:\forall s\in S:\qquad
P(r,s)=c(s)\cdot\sum_{s^{\prime}\in\omega(s)}P(r,s^{\prime})$

for constants $c(s)\in[0,1]$ which only depend on $s\in S$ (this is equivalent
to the existence of distributions $\alpha$ s.t. $\Omega$ and $\alpha$ are
deflatable).

When setting $\Pi$ as in (2.2) and $\alpha$ as in (2.5) (proportional
$\alpha$), then it holds that

$\displaystyle\tau(\rho)$
$\displaystyle=0\;\;\forall\rho\in\Omega\qquad\textrm{and}\qquad\alpha(s)=c(s)\;\;\forall
s\in S$

###### Proof

First, we show $\alpha(s)=c(s)$. Indeed, by (2.5),

$\displaystyle\alpha(s)$ $\displaystyle=\frac{\sum_{r\in S}P(r,s)}{\sum_{r\in
S}\sum_{s^{\prime}\in\omega(s)}P(r,s^{\prime})}=\frac{\sum_{r\in
S}c(s)\sum_{s^{\prime}\in\omega(s)}P(r,s^{\prime})}{\sum_{r\in
S}\sum_{s^{\prime}\in\omega(s)}P(r,s^{\prime})}=\frac{c(s)\sum_{r\in
S}\sum_{s^{\prime}\in\omega(s)}P(r,s^{\prime})}{\sum_{r\in
S}\sum_{s^{\prime}\in\omega(s)}P(r,s^{\prime})}=c(s)$

Again by 8, $\tau(\rho)=0$ for all $\rho$ is equivalent to

$\displaystyle\overbrace{\forall s,s^{\prime}\in S\textrm{
s.t.\leavevmode\nobreak\
}\omega(s)=\omega(s^{\prime}):\forall\rho\in\Omega:}^{\forall\ldots}\;\alpha(s^{\prime})\sum_{r\in\rho}\alpha(r)P(r,s)=\alpha(s)\sum_{r\in\rho}\alpha(r)P(r,s^{\prime})$
$\displaystyle\iff\;$
$\displaystyle\forall\ldots:\;\;c(s^{\prime})\sum_{r\in\rho}c(r)P(r,s)=c(s)\sum_{r\in\rho}c(r)P(r,s^{\prime})$
$\displaystyle\iff\;$
$\displaystyle\forall\ldots:\;\;c(s^{\prime})\sum_{r\in\rho}c(r)c(s)\sum_{\widehat{s}\in\omega(s)}P(r,\widehat{s})=c(s)\sum_{r\in\rho}c(r)c(s^{\prime})\sum_{\widehat{s}\in\omega(s^{\prime})}P(r,\widehat{s})$
$\displaystyle\iff\;$
$\displaystyle\forall\ldots:\;\;c(s^{\prime})c(s)\sum_{r\in\rho}c(r)\sum_{\widehat{s}\in\omega(s)}P(r,\widehat{s})=c(s)c(s^{\prime})\sum_{r\in\rho}c(r)\sum_{\widehat{s}\in\omega(s^{\prime})}P(r,\widehat{s})$

The statement in the last line ist true since $\omega(s)=\omega(s^{\prime})$
by the assumption on $s$ and $s^{\prime}$.

We next show that none of the lumpability concepts above are necessary
conditions for $\tau(\rho)=0$ for all $\rho$. Hence, except for [9], a large
part of the literature has treated stricter than necessary conditions in order
for dynamic-exact aggregation to be possible. None of the definitions of
ordinary and exact lumpability as well as deflatability take into account the
initial distribution, so none of these conditions are sufficient for an exact
aggregation. Exact lumpability and deflatability only imply dynamic-exactness.
At the same time, all three concepts (ordinary & exact lumpability,
deflatability) are still useful since they are easier to check
computationally, and can thus be relevant for practical applications. This
will be discussed in more detail in the next section.

###### Proposition 17

There are partitions $\Omega$ of the state space of a DTMC and probability
distributions $\alpha_{\sigma}$ with support on $\sigma\in\Omega$ which are
dynamic-exact (when $\Pi$ is set as in (2.2)), but where $\Omega$ is neither
ordinary lumpable, nor exactly lumpable, nor are $\Omega$ and the
distributions $\alpha$ deflatable.

###### Proof

We consider the state space $S=\\{1,2,3\\}$, the aggregation
$\Omega=\\{\\{1\\},\\{2,3\\}\\}$ and
$\alpha(1)=1,\alpha(2)=\frac{1}{4},\alpha(3)=\frac{3}{4}$ as well as the DTMC
given by the following transition matrix:

$\displaystyle P=\begin{pmatrix}0&\frac{1}{4}&\frac{3}{4}\\\
0&\frac{1}{2}&\frac{1}{2}\\\
\frac{4}{9}&\frac{1}{18}&\frac{1}{2}\end{pmatrix}$

Calculating $\Pi$ as in (2.2), we get

$\displaystyle\Pi=\begin{pmatrix}0&1\\\
\frac{1}{3}&\frac{2}{3}\end{pmatrix},\qquad A=\begin{pmatrix}1&0&0\\\
0&\frac{1}{4}&\frac{3}{4}\end{pmatrix},\qquad\Lambda=\begin{pmatrix}1&0\\\
0&1\\\ 0&1\end{pmatrix}$

It is easy to check that $\Pi A=AP$ holds, i.e. this aggregation is dynamic-
exact and $\tau(\\{1\\})=\tau(\\{2,3\\})=0$. However, we have

$\displaystyle\widetilde{P}=\Lambda\Pi
A=\begin{pmatrix}0&\frac{1}{4}&\frac{3}{4}\\\
\frac{1}{3}&\frac{1}{6}&\frac{1}{2}\\\
\frac{1}{3}&\frac{1}{6}&\frac{1}{2}\end{pmatrix}\neq P$

This already shows that $\Omega$ and $\alpha$ are not aggregatable. We now
show that none of the stated properties hold for $\Omega$ and $\alpha$:

  * •

ordinary lumpability: since $\omega(2)=\omega(3)$, by (4.1) in 11, we would
need

$\displaystyle\underbrace{P(2,2)+P(2,3)}_{=1}=\sum_{s\in\\{2,3\\}}P(2,s)=\sum_{s\in\\{2,3\\}}P(3,s)=\underbrace{P(3,2)+P(3,3)}_{=\frac{5}{9}}$

which is clearly not true.

  * •

exact lumpability: since $\omega(2)=\omega(3)$, by (4.2) in 12, we would need

$\displaystyle\underbrace{P(2,2)+P(3,2)}_{=\frac{5}{9}}=\sum_{r\in\\{2,3\\}}P(r,2)=\sum_{r\in\\{2,3\\}}P(r,3)=\underbrace{P(2,3)+P(3,3)}_{=1}$

which is clearly not true.

  * •

deflatability: since $\omega(2)=\omega(3)$, by (4.3) in 13, we would need

$\displaystyle\frac{1}{2}=P(2,2)=\alpha(2)\sum_{s\in\\{2,3\\}}P(2,s)=\alpha(2)\cdot
1\implies\alpha(2)=\frac{1}{2}$
$\displaystyle\frac{1}{18}=P(3,2)=\alpha(2)\sum_{s\in\\{2,3\\}}P(3,s)=\alpha(2)\cdot\frac{5}{9}\implies\alpha(2)=\frac{1}{10}$

so the given $\Omega$ and $\alpha$ are not deflatable, and there is even no
other choice of $\alpha$ such that $\Omega$ and $\alpha$ would be deflatable.

In fact, there is no aggregation $\Omega$ with $\Omega\neq\\{S\\}$ and
$\Omega\neq\\{\\{1\\},\\{2\\},\\{3\\}\\}$ (the trivial partitions) which is
ordinary or exactly lumpable, or for which deflatable $\alpha$ distributions
exist.

###### Remark

Note the following: for the example given in the proof above, neither setting
$\alpha$ as in (2.5) (proportional $\alpha$) nor setting $\alpha$ as in (2.6)
(uniform $\alpha$) will deliver the choice of $\alpha$ which results in a
dynamic-exact aggregation. For uniform $\alpha$, this is immediately clear,
and for proportional $\alpha$, we get

$\displaystyle\alpha(1)=1,\qquad\alpha(2)=\frac{29}{92}\approx
0.315,\qquad\alpha(3)=\frac{63}{92}\approx 0.685$

This shows that the proposed ways of calculating $\alpha$ can only be seen as
an approximation of the optimal choice for $\alpha$. By 14 and 16, the
proportional $\alpha$ calculation will, however, at least identify the optimal
choice if the aggregation is exactly lumpable or if a deflatable distribution
exists, which is why this choice of calculating $\alpha$ seems to be a good
compromise between optimality and computational efficiency.

Summarizing this section, we have seen that ordinary lumpability implies that
the aggregate probabilities $\pi_{k}$ (and $\pi_{t}$) are correct, but they
cannot necessarily be split into the correct probabilities for single states,
and we have seen that exact lumpability and deflatability both imply dynamic-
exactness. On the other hand, dynamic-exactness implies none of these three
concepts.

##  Choosing the aggregates

We want to choose a partition $\Omega$ of the state space $S$ of a DTMC or
CTMC in such a way that the error bounds $\tau(\rho)$ for $\rho\in\Omega$ are
small, as this will result in a good approximation $\widetilde{p}_{k}$ of the
transient distributions $p_{k}$. At the same time, we would like
$\left|\Omega\right|=m\ll n=\left|S\right|$ in order to reduce the
computational effort required to calculate $\widetilde{p}_{k}$. An ideal
algorithm would receive a parameter $\varepsilon$ as input and determine the
partition $\Omega$ with the fewest aggregates satisfying
$\sum_{\rho\in\Omega}\tau(\rho)<\varepsilon$ or
$\max_{\rho\in\Omega}\tau(\rho)<\varepsilon$ (the latter motivated by (3.2),
(3.3) and Theorem 4 i). This would guarantee a stepwise error (or error growth
rate for continuous time chains) of at most $\varepsilon$.

Solving the above problem exactly will in general result in a runtime
exceeding the time needed to simply compute $p_{k}$ exactly for the original
chain. We will therefore consider different ways of choosing an $\Omega$ which
is somehow close to the optimal solution. The presented algorithms can be used
for general irreducible Markov chains, without any assumption on the structure
of the chain. However, they can only be expected to perform well if a certain
structure _is_ present: Section 3.3 and Section 4 characterize settings in
which the error bounds are low, and which can be identified with comparatively
little computational effort, in contrast to finding a partition which
satisfies 2. The algorithms presented in this section will thus try to
identify aggregates which are close to fulfilling conditions such as exact
lumpability or aggregatability. We will start with some remarks on the time
needed to calculate the exact transient distribution of a Markov chain, in
order to be able to assess the speed-up resulting from aggregation.

###  No aggregation

If we perform no aggregation and compute exact transient distributions, then
we get the following runtimes:

  * •

For DTMCs, computing $p_{k}$ amounts to $k$ vector-matrix multiplications
(each of the form $p_{i+1}=p_{i}P$) of a vector of length $n$ with a matrix of
size $n\times n$. Each such multiplication has a runtime of order
$\mathcal{O}(n^{2})$, and since we need $k$ of those, the total runtime
amounts to $\mathcal{O}(kn^{2})$.

  * •

For CTMCs, computing the exact transient distribution $p_{t}$ at time $t$ is
not feasible in general. Instead, using uniformisation and truncation (see [1,
Section 2.5] and [7]), the transient distribution can be computed up to a pre-
defined error $\varepsilon$ (we will neglect the dependence of the runtime on
$\varepsilon$, which we assume to be fixed). The runtime of this computation
is $\mathcal{O}(qtn^{2})$ where $q:=\max_{s\in S}\left|Q(s,s)\right|$.

###  Almost aggregatability

In 13, aggregatability for DTMCs was defined, and we have seen in 16 that
deflatability actually already implies that the error bound is zero. The first
algorithm, taken from [2], will thus try to identify a partition which is
close to an aggregatable partition (or “almost aggregatable”, see [2,
Definition 2.8]). As the definition of deflatability only applies to DTMCs,
this partitioning algorithm only works for DTMCs a priori.

We quickly summarize the algorithm from [2] here for better understanding of
some of the choices available when implementing this algorithm. Note that in
[2], the transpose of the transition matrix $P$ is considered, and all
transient distributions are column vectors instead of row vectors as in this
paper. In the following summary, we will stick to our notation, and hence
present a transposed version of the results in [2].

We first consider the case of an aggregatable matrix $P$, which means that
$\widetilde{P}\;\overset{\textrm{def.}}{=}\;\Lambda\Pi A=P$ by [2, Proposition
2.6]. The idea is to find a connection between the singular value
decomposition of $\Pi$ and $P$ which allows for identification of the
aggregates by analysing the singular value decomposition of $P$, without
knowledge of $\Pi$ or $\Omega$. For an almost aggregatable $P$, the changes in
the singular value decomposition compared to the closest aggregatable matrix
will be small, and the algorithm can thus also be applied if a matrix is only
almost aggregatable. For a detailed justification, we refer to [2].

In order to understand the connection between the singular value decomposition
and $\Omega$, we assume for a moment that $P$ is aggregatable and that
$\Omega$ and the $\alpha$ distributions are given. Consider matrices

$\displaystyle
D_{A}=\begin{pmatrix}\left\lVert\alpha_{\Omega_{1}}\right\rVert_{2}&&\\\
&\ddots&\\\
&&\left\lVert\alpha_{\Omega_{m}}\right\rVert_{2}\end{pmatrix}\qquad\qquad
D_{\Lambda}$ $\displaystyle=\begin{pmatrix}\sqrt{\left|\Omega_{1}\right|}&&\\\
&\ddots&\\\ &&\sqrt{\left|\Omega_{m}\right|}\end{pmatrix}$

Let $\widehat{U}\widehat{\Sigma}\widehat{V}$ be the singular value
decomposition of $D_{\Lambda}\Pi D_{A}$ (note: we do _not_ transpose
$\widehat{V}$), i.e. $\widehat{U}\widehat{\Sigma}\widehat{V}=D_{\Lambda}\Pi
D_{A}$, $\widehat{U}$ and $\widehat{V}$ are orthogonal (their rows and columns
are orthonormal vectors), and $\widehat{\Sigma}$ contains the singular values
on its diagonal, ordered from largest to smallest value. Then, the singular
value decomposition $U\Sigma V$ of $P$ can be written as

$\displaystyle U=\left(\begin{array}[]{c|c}\Lambda
D_{\Lambda}^{-1}\widehat{U}&U^{(2)}\end{array}\right)\in\mathbb{R}^{n\times
n}\qquad\Sigma=\left(\begin{array}[]{c|c}\widehat{\Sigma}&0\\\ \hline\cr
0&0\end{array}\right)\in\mathbb{R}^{n\times n}\qquad
V=\left(\begin{array}[]{c}\widehat{V}D_{A}^{-1}A\\\ \hline\cr
V^{(2)}\end{array}\right)\in\mathbb{R}^{n\times n}$

where $U^{(2)}$ and $V^{(2)}$ are such that $U$ and $V$ are orthogonal. This
is verified by computing:

$\displaystyle U\Sigma
V=\Lambda\overbrace{D_{\Lambda}^{-1}\underbrace{\widehat{U}\widehat{\Sigma}\widehat{V}}_{D_{\Lambda}\Pi
D_{A}}D_{A}^{-1}}^{\Pi}A=\widetilde{P}=P$

Details can be found in [2], where orthogonality of $U$ and $V$ is shown as
well. Note that we used the matrices $D_{\Lambda}$ and $D_{A}$ only in order
to scale the singular values of $\Pi$ such that they agree with the singular
values of $P$, they have no other relevance. For simplicity, we will assume
that the singular value decomposition of $P$ is unique if the singular values
in $\Sigma$ are ordered by size and if we ignore $U^{(2)}$ and $V^{(2)}$ which
are anyway not used anymore. This is true, for example, when $\Pi$ has full
rank (i.e. rank $m$) and all non-zero singular values of $P$ have multiplicity
one.

We now consider the first $m$ rows of $V$ (which are the first $m$ right-
singular vectors of $P$), i.e. $\widehat{V}D_{A}^{-1}A\in\mathbb{R}^{m\times
n}$. We call the $i$-th column of this submatrix
$v(i)=\left(V(1,i),\ldots,V(m,i)\right)\in\mathbb{R}^{m}$. Note that $i$
corresponds to state $i$ of the Markov chain, and hence $\omega(i)$ is the
corresponding aggregate. Using that the upper $m$ rows of $V$ are equal to
$\widehat{V}D_{A}^{-1}A$, we can write

$\displaystyle
v(i)=\left\lVert\alpha_{\omega(i)}\right\rVert_{2}^{-1}\alpha(i)\cdot\underbrace{\left(\widehat{V}(1,\omega(i)),\ldots,\widehat{V}(m,\omega(i))\right)}_{\omega(i)\textrm{-th
column of }\widehat{V}}\in\mathbb{R}^{m}$

This implies that for two states $r,s$ in the same aggregate, $v(r)$ and
$v(s)$ will point in the same direction – in fact,

$\displaystyle\forall r,s\in S\textrm{ s.t.\leavevmode\nobreak\
}\omega(r)=\omega(s):\qquad\mathbb{R}^{m}\ni v(r)$
$\displaystyle=\underbrace{\frac{\alpha(r)}{\alpha(s)}}_{\in\mathbb{R}}\cdot
v(s)\in\mathbb{R}^{m}$ (5.1)

Furthermore, by orthogonality of the columns of $\widehat{V}$, $v(r)$ and
$v(s)$ will be orthogonal if $\omega(r)\neq\omega(s)$. This structure of the
submatrix can be exploited to recover $\Omega$ (and even $\alpha$). For almost
aggregatable $P$, we can only expect (5.1) to hold approximately. This leads
to the following three possible algorithms to compute $\Omega$ if $m$ is
already fixed:

  * •

SVD sgn: Proposed as a very simple algorithm in [2] with only limited
practical applicability due to its numerical instability. The aggregates are
recovered by putting two states $r$ and $s$ into the same aggregate if the
sign structure of the vectors $v(r)$ and $v(s)$ is identical. By (5.1) and by
orthogonality of the vectors $v$ for states in different aggregates, this
yields the correct partition $\Omega$ if $P$ is aggregatable. However, for
almost aggregatable $P$, perturbed values in the vectors $v$ can lead to the
sign of an entry changing, resulting in instability.

  * •

SVD seba: Proposed as a more stable algorithm in [2] via a combination with
[8]. The sparse eigenbasis approximation algorithm proposed in [8] is applied
to the first $m$ rows of $V$. This results in an approximate sparse basis of
the space spanned by the first $m$ rows of $V$. In the case of aggregatability
of $P$, the space spanned by the first $m$ rows of $V$ is spanned by the
vectors $\alpha_{\Omega_{1}},\ldots,\alpha_{\Omega_{m}}$ as a consequence of
(5.1). Therefore, the sparse basis obtained from applying the algorithm of [8]
should approximately correspond to the $\alpha$ distributions for almost
aggregatable $P$.

  * •

SVD dir: A new proposal presented in this paper, with the intention to fully
exploit (5.1). States are clustered such that the distance between the
corresponding $v$ vectors within a cluster should be low, where the distance
is measured as follows: the shorter vector is projected onto the longer
vector, and we then measure the euclidean distance between the shorter
original and the projected vector (see Algorithm 1 below). The intention
behind this is to measure whether two vectors $v(r)$ and $v(s)$ point in
approximately the same direction, as should be the case for $r$ and $s$ in the
same aggregate by (5.1). The distance between vectors is not measured as the
angle between vectors because this approach would suffer from the same
numerical instability as SVD sgn: if the entries of a $v$ vector are close to
$0$, small perturbations can lead to huge changes in the angle of the vector.
An additional ordering of the $v$ vectors by length is applied to increase
stability. We give the implementation details in Section 5.2.1.

With given $m$, the algorithms above can recover $\Omega$ from the singular
value decomposition of an almost aggregatable matrix $P$. However, $m$ is in
general not known in advance. In order to choose $m$, it makes sense to
analyse the spectrum of singular values of $P$. One way would be to identify a
gap in this spectrum and set $m$ to the number of singular values above the
gap. We chose a different approach which delivered better results in
experiments. Call the singular values $\gamma_{1}\geq\ldots\geq\gamma_{n}\geq
0$. Given a threshold parameter $\varepsilon$, $m$ is set to the smallest
value such that

$\displaystyle\sum_{i=1}^{m}\gamma_{i}$
$\displaystyle\geq(1-\varepsilon)\sum_{i=1}^{n}\gamma_{i}$ (5.2)

i.e. the first $m$ singular values sum to at least $1-\varepsilon$ times the
sum of all singular values.

The runtime of the SVD algorithm depends on the variant chosen. For the simple
SVD sgn, the computational cost is dominated by the singular value
decomposition, which needs $\mathcal{O}(n^{3})$ time (see [10, Figure 8.6.1]).
We saw that computing $p_{k}$ without aggregation results in a complexity of
$\mathcal{O}(kn^{2})$. Asymptotically (for large $k$ and $n$), applying the
SVD algorithm therefore only makes sense if $k\gg n$. As noted in [2], we can
reduce the complexity of the SVD approach by random sampling of the entries of
$P$ under some assumptions, which makes the algorithm more attractive. See [2]
for details.

As mentioned before, the three variants of the SVD algorithm can only be
applied to DTMCs. They are expected to perform well for almost aggregatable
$P$, and this is confirmed in the experiments. However, in different settings,
identifying almost exactly lumpable partitions may provide better results,
with the additional benefit that this is also possible for CTMCs. This is
discussed in Section 5.3.

####  Implementation of SVD dir

We briefly summarize the implementation details of SVD dir for completeness.
First, we define the function which computes if two vectors point in
approximately the same direction, which should also work if vectors with
entries close to zero are perturbed randomly. In particular, just using the
angle between vectors is not a good distance measure for such perturbed
vectors, since the angle of a vector with small entries can change
dramatically when the vector entries are perturbed. Instead, we use the
distance measure vdist, calculated as in the following Algorithm 1.

Input: $\textrm{vectors }v_{1},v_{2}\in\mathbb{R}^{m}\textrm{
s.t.\leavevmode\nobreak\ }\left\lVert v_{1}\right\rVert_{2}\geq\left\lVert
v_{2}\right\rVert_{2}$
Output: $\textrm{a real number }\geq 0\textrm{, which roughly measures a
distance between vector angles}$

1:$v_{\textrm{proj}}\leftarrow\textrm{projection of }v_{2}\textrm{ onto
}v_{1}$

2:return $\left\lVert v_{\textrm{proj}}-v_{2}\right\rVert_{2}$

Algorithm 1 vdist (measuring whether vectors point in a similar direction)

We can now look at the implementation details of SVD dir, given in Algorithm
2. After calculating the singular value decomposition of $P$ and after
determining the number of rows of $V$ to consider with the aid of (5.2), the
algorithm iterates over all states $s\in S$ in descending order of the length
of $v(s)$ (in line 13 of Algorithm 2). This ordering is used to increase
stability: the larger the entries of a vector, the less its angle varies under
small perturbations of the vector components. Hence, the first assignments of
states to aggregates are based on longer vectors whose angles are most
reliable.

The assignment to aggregates proceeds as follows: when processing state $s$,
the minimum distance $d_{\textrm{min}}$ of $v(s)$ (according to vdist, see
Algorithm 1) to the vectors $v(s_{1}),\ldots,v(s_{i})$ is calculated, where
$s_{1},\ldots,s_{i}$ are all states which have already been processed. If this
distance is smaller than the parameter $\delta$ (i.e. if
$d_{\textrm{min}}<\delta$), then $s$ is assigned to the aggregate of
$s_{i^{\ast}}$ where the distance $\textrm{vdist}(v(s),v(s_{i^{\ast}}))$ is
minimal. In another effort to increase numerical stability, the distance
$d_{\textrm{min}}$ is first only calculated for the states in
$S_{\textrm{reliable}}=\\{s^{\prime}\in\\{s_{1},\ldots,s_{i}\\}:\left\lVert
v(s^{\prime})\right\rVert_{2}>2\delta\\}$, and an assignment to an aggregate
of one of the states in $S_{\textrm{reliable}}$ is performed if
$d_{\textrm{min}}<\delta$ on this set, and only if this is not possible, we
calculate distances to all previously processed states. This additional step
follows the same reasoning as before: the angle of short perturbed vectors is
uncertain, so the aggregate assignment is more certain if we assign $s$ to an
aggregate which contains a state with a longer vector pointing in a similar
direction as $v(s)$. If no previously processed state $s^{\prime}$ satisfies
$\textrm{vdist}(v(s),v(s^{\prime}))<\delta$, then $s$ is assigned to a new
aggregate.

Input: $\textrm{a Markov chain, defined via its transition matrix }P\textrm{
on state space }S=\\{1,\ldots,n\\}\textrm{,}$
Input: $\textrm{the parameter }\varepsilon\textrm{, and the parameter }\delta$
Output: $\textrm{an aggregation function }\omega$
Output: whose corresponding partition is close to an aggregatable partition

1:$\omega\leftarrow\left((s\in S)\mapsto 1\right)$ $\triangleright$
aggregation function

2:$U,\Sigma,V\leftarrow\textrm{singular value decomposition of }P$
$\triangleright$ s.t. $P=U\Sigma V$, _not_ $P=U\Sigma V^{\mathsf{T}}$, with
ordered singular values in $\Sigma$

3:$\gamma_{1},\ldots,\gamma_{n}\leftarrow\textrm{values on diagonal of
}\Sigma$ $\triangleright$ $\gamma_{1}\geq\ldots\geq\gamma_{n}\geq 0$ are the
singular values of $P$

4:$m\leftarrow 0$ $\triangleright$ number of considered rows of $V$

5:while $\sum_{i=1}^{m}\gamma_{i}<(1-\varepsilon)\sum_{i=1}^{n}\gamma_{i}$ do
$\triangleright$ determine $m$ according to (5.2)

6: $m\leftarrow m+1$

7:end while

8:for all $s\in\\{1,\ldots,n\\}$ do

9: $v(s)\leftarrow(V(1,s),\ldots,V(m,s))\in\mathbb{R}^{m}$ $\triangleright$
extract cropped columns from $V$

10:end for

11:$m_{\textrm{current}}\leftarrow 0$ $\triangleright$ current number of
aggregates

12:$(s_{1},\ldots,s_{n})\;\leftarrow\;$ permutation of $(1,\ldots,n)$ s.t.
$\left\lVert v(s_{1})\right\rVert_{2}\geq\ldots\geq\left\lVert
v(s_{n})\right\rVert_{2}$

13:for all $i\in\\{1,\ldots,n\\}\textrm{ (in order)}$ do $\triangleright$ go
through $v(s_{i})$ by descending length

14:
$S_{\textrm{reliable}}\leftarrow\\{s\in\\{s_{1},\ldots,s_{i-1}\\}:\left\lVert
v(s)\right\rVert_{2}>2\delta\\}$

15: $S_{\lnot\textrm{reliable}}\leftarrow\\{s_{1},\ldots,s_{i-1}\\}\setminus
S_{\textrm{reliable}}$

16: $d_{\textrm{min}}\leftarrow\min_{s\in
S_{\textrm{reliable}}}\textrm{vdist}(v(s_{i}),v(s))$ $\triangleright$ see
Algo. 1, $d_{\textrm{min}}=\infty$ for $S_{\textrm{reliable}}=\varnothing$

17: if $d_{\textrm{min}}<\delta$ then

18: $\omega(s_{i})\leftarrow\omega\left(\operatorname*{arg\,min}_{s\in
S_{\textrm{reliable}}}\textrm{vdist}(v(s_{i}),v(s))\right)$ $\triangleright$
$s_{i}$ added to aggregate which contains closest vector

19: else

20: $d_{\textrm{min}}\leftarrow\min_{s\in
S_{\lnot\textrm{reliable}}}\textrm{vdist}(v(s_{i}),v(s))$ $\triangleright$ try
again with less reliable vectors

21: if $d_{\textrm{min}}<\delta$ then

22: $\omega(s_{i})\leftarrow\omega\left(\operatorname*{arg\,min}_{s\in
S_{\lnot\textrm{reliable}}}\textrm{vdist}(v(s_{i}),v(s))\right)$

23: else

24: $m_{\textrm{current}}\leftarrow m_{\textrm{current}}+1$

25: $\omega(s_{i})\leftarrow m_{\textrm{current}}$ $\triangleright$ $s_{i}$
added to new aggregate

26: end if

27: end if

28:end for

29:return $\omega$

Algorithm 2 Calculating almost aggregatable partitions with SVD dir

As input for Algorithm 2, $\varepsilon$ should be chosen according to how
close the resulting aggregation should be to an aggregatable partition. For
$\varepsilon=0$ (and $\delta=0$), the algorithm finds aggregatable partitions.
For $\varepsilon=1$, all states are assigned to a single aggregate. We used
$\delta=0.05$ in all our experiments. Using a constant $\delta$ regardless of
the size of the state space makes sense insofar as that the columns of $V$ are
vectors of unit length. However, since we crop the columns to dimension $m$
(which results in shorter vectors in general), a $\delta$ depending on $m$
might also be a good choice. Further fine tuning is necessary here. The same
holds for the choice to put a state $s$ into $S_{\textrm{reliable}}$ if
$\left\lVert v(s)\right\rVert_{2}>2\delta$ in line 14 of Algorithm 2. Other
choices should be investigated. An in-depth analysis of the numerical
stability of Algorithm 2 is also missing, but we observed at least a
convincing performance in our experiments.

We conclude this section with a short runtime analysis of Algorithm 2. Line 2
(the singular value decomposition) takes $\mathcal{O}(n^{3})$ time, as stated
before (see [10, Figure 8.6.1]). Lines 1 and 3 to 11 are negligible in
comparison. Line 12 (sorting the vectors by descending length) takes time
$\mathcal{O}(n\log(n))$ after having computed all vector lengths in time
$\mathcal{O}(nm)$; both runtimes are also smaller than $\mathcal{O}(n^{3})$.
In line 13, we loop over all $n$ states. Determining $S_{\textrm{reliable}}$
(lines 14 and 15) can be done in time $\mathcal{O}(1)$ by successively adding
each processed state with suitable length to $S_{\textrm{reliable}}$. Line 16
takes $\mathcal{O}(nm)$ time at most (we need to call vdist at most $n$ times,
and vdist is applied to vectors of dimension $m$). The same holds for line 20.
The remaining part of the loop body is negligible. We arrive at
$\mathcal{O}(nm)$ for the loop body, giving a total runtime of the loop in
lines 13 to 28 of $\mathcal{O}(n^{2}m)$ which is also smaller than
$\mathcal{O}(n^{3})$. The overall runtime is thus still dominated by the
singular value decomposition with $\mathcal{O}(n^{3})$.

###  $\varepsilon$-almost exact lumpability

By 14, if $\Omega$ is exactly lumpable, $\Pi$ is set as in (2.2) and $\alpha$
as in (2.5) or (2.6) (proportional or uniform $\alpha$), then the error bound
is zero. However, we cannot expect an exactly lumpable partition to exist for
a general Markov chain. It might be, though, that a partition exists which is
close to being exactly lumpable. This motivates the following definition:

###### Definition 18

We call a partition $\Omega$ $\varepsilon$-almost exactly lumpable if:

$\displaystyle\forall s,s^{\prime}\in S\textrm{ s.t.\leavevmode\nobreak\
}\omega(s)$
$\displaystyle=\omega(s^{\prime}):\qquad\sum_{\rho\in\Omega}\left|\sum_{r\in\rho}P(r,s)-\sum_{r\in\rho}P(r,s^{\prime})\right|\leq\varepsilon$

The incoming probabilities to two states in the same aggregate from another
aggregate are not required to be identical anymore in 18, but they are close
to being identical. For the SVD partitioning algorithm, which also takes a
(different) parameter $\varepsilon$ used for cutting off the smallest singular
values (see (5.2)), it seems to be difficult to derive a bound on the $\tau$
factors depending on the input parameter $\varepsilon$. For
$\varepsilon$-almost exactly lumpable partitions, we have at least the
following result:

###### Proposition 19

Given a partition $\Omega$ of the state space of a DTMC or CTMC which is
$\varepsilon$-almost exactly lumpable, let $\alpha$ as in (2.6) (uniform
$\alpha$) and $\Pi$ or $\Theta$ as in (2.2). Then

$\displaystyle\sum_{\rho\in\Omega}\tau(\rho)$
$\displaystyle\leq\left|\Omega\right|\cdot\frac{\max_{\rho\in\Omega}\left|\rho\right|}{\min_{\rho\in\Omega}\left|\rho\right|}\cdot\varepsilon$

Hence, given an algorithm which takes $\varepsilon$ as input and outputs an
$\varepsilon$-almost exactly lumpable partition, we can choose some
$\varepsilon$ which guarantees a desired error bound $\tau$ in advance without
actually running the algorithm and calculating the $\tau$ factors for the
resulting partition. In practice, the bound in 19 always seemed to be much
larger than the actual value of $\sum_{\rho\in\Omega}\tau(\rho)$, though.

######  Proof of (19)

We have

$\displaystyle\begin{split}\sum_{\rho\in\Omega}\tau(\rho)&=\sum_{\rho\in\Omega}\sum_{\sigma\in\Omega}\sum_{s\in\sigma}\left|\alpha(s)\Pi(\rho,\sigma)-\sum_{r\in\rho}\alpha(r)P(r,s)\right|\\\
&=\sum_{\rho\in\Omega}\sum_{\sigma\in\Omega}\sum_{s\in\sigma}\left|\frac{1}{\left|\sigma\right|}\sum_{r\in\rho}\frac{1}{\left|\rho\right|}\sum_{s^{\prime}\in\sigma}P(r,s^{\prime})-\sum_{r\in\rho}\frac{1}{\left|\rho\right|}P(r,s)\right|\\\
&=\sum_{\rho\in\Omega}\sum_{\sigma\in\Omega}\sum_{s\in\sigma}\left|\frac{1}{\left|\sigma\right|}\frac{1}{\left|\rho\right|}\sum_{s^{\prime}\in\sigma}\sum_{r\in\rho}P(r,s^{\prime})-\frac{1}{\left|\sigma\right|}\frac{1}{\left|\rho\right|}\sum_{s^{\prime}\in\sigma}\sum_{r\in\rho}P(r,s)\right|\\\
&=\sum_{\rho\in\Omega}\sum_{\sigma\in\Omega}\frac{1}{\left|\sigma\right|}\frac{1}{\left|\rho\right|}\sum_{s\in\sigma}\left|\sum_{s^{\prime}\in\sigma}\sum_{r\in\rho}\left(P(r,s^{\prime})-P(r,s)\right)\right|\\\
&\leq\sum_{\sigma\in\Omega}\frac{1}{\left|\sigma\right|}\sum_{s\in\sigma}\sum_{s^{\prime}\in\sigma}\sum_{\rho\in\Omega}\frac{1}{\left|\rho\right|}\left|\sum_{r\in\rho}P(r,s^{\prime})-\sum_{r\in\rho}P(r,s)\right|\\\
&\leq\sum_{\sigma\in\Omega}\frac{1}{\left|\sigma\right|}\sum_{s\in\sigma}\sum_{s^{\prime}\in\sigma}\frac{1}{\min_{\rho\in\Omega}\left|\rho\right|}\underbrace{\sum_{\rho\in\Omega}\left|\sum_{r\in\rho}P(r,s^{\prime})-\sum_{r\in\rho}P(r,s)\right|}_{\leq\varepsilon}\\\
&\leq\sum_{\sigma\in\Omega}\frac{\left|\sigma\right|}{\min_{\rho\in\Omega}\left|\rho\right|}\cdot\varepsilon\leq\left|\Omega\right|\cdot\frac{\max_{\rho\in\Omega}\left|\rho\right|}{\min_{\rho\in\Omega}\left|\rho\right|}\cdot\varepsilon\end{split}$
(5.3)

The same calculation holds for continuous time with $P$ replaced by $Q$ and
$\Pi$ replaced by $\Theta$. Also note that the bound can be improved to

$\displaystyle\sum_{\rho\in\Omega}\tau(\rho)$
$\displaystyle\leq\left|\Omega\right|\cdot\frac{\max_{\rho\in\Omega}\left|\rho\right|-1}{\min_{\rho\in\Omega}\left|\rho\right|}\cdot\varepsilon$
(5.4)

with the same calculation by noting that the double sum
$\sum_{s\in\sigma}\sum_{s^{\prime}\in\sigma}$ in (5.3) actually sums over
elements which are zero for $s=s^{\prime}$.

###### Example 20

The bound given in 19 cannot be significantly improved. We can give an example
which shows that $\sum_{\rho\in\Omega}\tau(\rho)$ can grow linearly in
$\left|\Omega\right|$ for constant $\varepsilon$. Consider the state space
$S=\\{1,\ldots,n\\}$ with $n=2m$, the partition
$\Omega=\\{\\{1,2\\},\\{3,4\\},\ldots,\\{2m-1,2m\\}\\}$ and $\alpha$ set as in
(2.6) (uniform $\alpha$). Define $P(r,2s)=\frac{2+\varepsilon}{2n}$ and
$P(r,2s-1)=\frac{2-\varepsilon}{2n}$ for $r\in S$ and $s\in\\{1,\ldots,m\\}$.
This partition is $\varepsilon$-almost exactly lumpable. Indeed, consider two
states $s^{\prime},s^{\prime\prime}$ in the same aggregate. Then, w.l.o.g., we
have $s^{\prime}=2s$ and $s^{\prime\prime}=2s-1$ for some
$s\in\\{1,\ldots,m\\}$. Hence

$\displaystyle\sum_{\rho\in\Omega}\left|\sum_{r\in\rho}P(r,s^{\prime})-\sum_{r\in\rho}P(r,s^{\prime\prime})\right|$
$\displaystyle=\sum_{\rho\in\Omega}\left|\sum_{r\in\rho}P(r,2s)-\sum_{r\in\rho}P(r,2s-1)\right|$
$\displaystyle=\sum_{\rho\in\Omega}\left|\sum_{r\in\rho}\frac{2+\varepsilon}{2n}-\sum_{r\in\rho}\frac{2-\varepsilon}{2n}\right|$
$\displaystyle=\left|\Omega\right|\cdot
2\cdot\left(\frac{2+\varepsilon}{2n}-\frac{2-\varepsilon}{2n}\right)=\frac{n}{2}\cdot
2\cdot\frac{\varepsilon}{n}=\varepsilon$

Now, consider $r,s\in\\{1,\ldots,m\\}$. Setting
$\rho=\\{2r-1,2r\\},\sigma=\\{2s-1,2s\\}$, we have

$\displaystyle\Pi(\rho,\sigma)$
$\displaystyle=\frac{1}{2}\cdot\left(P(2r-1,2s-1)+P(2r-1,2s)\right)+\frac{1}{2}\cdot\left(P(2r,2s-1)+P(2r,2s)\right)=\frac{2}{n}$

Hence, for $\rho=\\{2r-1,2r\\}\in\Omega$, it holds that

$\displaystyle\tau(\rho)$
$\displaystyle=\sum_{\sigma\in\Omega}\sum_{s\in\sigma}\left|\alpha(s)\Pi(\rho,\sigma)-\sum_{r\in\rho}\alpha(r)P(r,s)\right|$
$\displaystyle=\sum_{s=1}^{\left|\Omega\right|}\left(\left|\alpha(2s-1)\cdot\frac{2}{n}-\sum_{r\in\rho}\alpha(r)P(r,2s-1)\right|+\left|\alpha(2s)\cdot\frac{2}{n}-\sum_{r\in\rho}\alpha(r)P(r,2s)\right|\right)$
$\displaystyle=\sum_{s=1}^{\left|\Omega\right|}\left(\left|\frac{1}{2}\cdot\frac{2}{n}-\frac{2-\varepsilon}{2n}\right|+\left|\frac{1}{2}\cdot\frac{2}{n}-\frac{2+\varepsilon}{2n}\right|\right)=\sum_{s=1}^{\left|\Omega\right|}\frac{\varepsilon}{n}=\frac{\left|\Omega\right|}{n}\cdot\varepsilon=\frac{\varepsilon}{2}$

We thus get
$\sum_{\rho\in\Omega}\tau(\rho)=\left|\Omega\right|\cdot\frac{1}{2}\cdot\varepsilon$.
Therefore, we cannot drop the dependence on $\left|\Omega\right|$ in the bound
given in 19. Actually, the improved bound given in (5.4) is tight in this case
as we have
$\min_{\rho\in\Omega}\left|\rho\right|=\max_{\rho\in\Omega}\left|\rho\right|=2$.

We now develop an algorithm which finds an $\varepsilon$-almost exactly
lumpable partition as a counterpart to the SVD approach for almost
aggregatable partitions. The algorithm works for both DTMCs as well as CTMCs
(we give the DTMC version, but for CTMCs, $P$ only has to be replaced by $Q$).
For a given $\varepsilon$, the algorithm should find a partition which is as
coarse as possible and still satisfies $\varepsilon$-almost exact lumpability.
Note that in general, there is no unique coarsest $\varepsilon$-almost exactly
lumpable partition. However, there always is a unique coarsest _exactly
lumpable_ partition which may be found by successive refinement of the
partition $\Omega=\\{S\\}$. We can thus hope to get good results by using a
successive refinement algorithm for $\varepsilon$-_almost exact lumpability_
as well. For completeness and better understanding, we first show that there
always is a unique coarsest _exactly lumpable_ partition (which was already
shown in [4], where the partition refinement approach is mentioned on p. 269).
This justifies the structure of the algorithm below for practical purposes.

###### Lemma 21

Assume that the partition $\Omega$ satisfies exact lumpability. Let
$\widetilde{\Omega}$ be a different partition of the state space which
satisfies:

$\displaystyle\forall\widetilde{\sigma}\in\widetilde{\Omega}:\exists\sigma_{1},\ldots,\sigma_{k}\in\Omega:\widetilde{\sigma}=\bigcup_{i=1}^{k}\sigma_{i}$

i.e. $\widetilde{\Omega}$ is coarser than $\Omega$ in the sense that the
aggregated states in partition $\Omega$ are subsets of the aggregated states
in $\widetilde{\Omega}$.

Then, we have for DTMCs:

$\forall s,s^{\prime}\in S\textrm{ s.t.\leavevmode\nobreak\
}\omega(s)=\omega(s^{\prime}):\forall\widetilde{\sigma}\in\widetilde{\Omega}:\qquad\sum_{r\in\widetilde{\sigma}}P(r,s)=\sum_{r\in\widetilde{\sigma}}P(r,s^{\prime})$
(5.5)

The same holds for CTMCs.

###### Remark

The partition $\widetilde{\Omega}=\\{S\\}$ is always coarser than any exactly
lumpable partition $\Omega$, hence 21 implies that

$\displaystyle\forall s,s^{\prime}\in S\textrm{ s.t.\leavevmode\nobreak\
}\omega(s)$ $\displaystyle=\omega(s^{\prime}):\qquad\sum_{r\in
S}P(r,s)=\sum_{r\in S}P(r,s^{\prime})$

Also note that 21 does _not_ imply that $\widetilde{\Omega}$ is exactly
lumpable as well since in (5.5), the aggregation function $\omega$ corresponds
to the partition $\Omega$ and not to $\widetilde{\Omega}$.

######  Proof of (21)

This is easy to see. Let $\widetilde{\sigma}\in\widetilde{\Omega}$ be
arbitrary and $\sigma_{1},\ldots,\sigma_{k}\in\Omega$ s.t.
$\widetilde{\sigma}=\bigcup_{i=1}^{k}\sigma_{i}$. Then, for $s,s^{\prime}\in
S$ s.t. $\omega(s)=\omega(s^{\prime})$, we have

$\displaystyle\sum_{r\in\widetilde{\sigma}}P(r,s)$
$\displaystyle=\sum_{i=1}^{k}\sum_{r\in\sigma_{i}}P(r,s)=\sum_{i=1}^{k}\sum_{r\in\sigma_{i}}P(r,s^{\prime})=\sum_{r\in\widetilde{\sigma}}P(r,s^{\prime})$

The same calculation holds for CTMCs.

###### Proposition 22

For every DTMC or CTMC, there exists a unique coarsest exactly lumpable
partition.

###### Proof

The partition where every aggregate contains exactly one state is always
exactly lumpable, so we know that an exactly lumpable partition exists.

Let us call $\Omega^{(1)}=\\{S\\}$. We will now construct a sequence of
partitions $\Omega^{(1)},\ldots,\Omega^{(k)}$ such that $\Omega^{(i+1)}$ is
finer than $\Omega^{(i)}$ as follows. Assume we already have constructed
$\Omega^{(i)}$. We now construct $\Omega^{(i+1)}$. For $s,s^{\prime}\in S$
with $s\neq s^{\prime}$, we

  1. (i)

assign $s$ and $s^{\prime}$ to different aggregates in $\Omega^{(i+1)}$ if
they already belong to different aggregates in $\Omega^{(i)}$,

  2. (ii)

assign $s$ and $s^{\prime}$ to different aggregates in $\Omega^{(i+1)}$ if
they belong to the same aggregate in $\Omega^{(i)}$ but there exists some
$\sigma\in\Omega^{(i)}$ such that

$\displaystyle\sum_{r\in\sigma}P(r,s)\neq\sum_{r\in\sigma}P(r,s^{\prime})$

  3. (iii)

assign $s$ and $s^{\prime}$ to the same aggregate in $\Omega^{(i+1)}$
otherwise.

We stop the construction as soon as
$\Omega^{(i+1)}=\Omega^{(i)}=:\Omega^{(k)}$. Note that this needs to happen at
some point, at the latest when every aggregate in $\Omega^{(i)}$ contains only
one state.

Now,

  * •

if $\Omega^{(i+1)}=\Omega^{(i)}$, then $\Omega^{(i+1)}$ is exactly lumpable
since step ii was not applied in the last iteration of the construction.
Hence, by definition, $\Omega^{(i+1)}=\Omega^{(i)}=\Omega^{(k)}$ is exactly
lumpable.

  * •

for any exactly lumpable partition $\Omega$, we have that $\Omega$ is finer
than (or equal to) $\Omega^{(i)}$ for all $i=1,\ldots,k$. We prove this by
induction. $\Omega^{(1)}=\\{S\\}$ is coarser than any partition, so the
statement holds for $i=1$. The step from $i$ to $i+1$ is done as follows: if
$\Omega^{(i)}$ is coarser than the exactly lumpable $\Omega$, then, by 21, we
have that

$\displaystyle\forall s,s^{\prime}\in S\textrm{ s.t.\leavevmode\nobreak\
}\omega(s)$
$\displaystyle=\omega(s^{\prime}):\forall\sigma\in\Omega^{(i)}:\qquad\sum_{r\in\sigma}P(r,s)=\sum_{r\in\sigma}P(r,s^{\prime})$

Hence, when calculating $\Omega^{(i+1)}$, step ii will only assign $s$ and
$s^{\prime}$ to different aggregates in $\Omega^{(i+1)}$ if
$\omega(s)\neq\omega(s^{\prime})$, i.e. if they also belong to different
aggregates in $\Omega$. The same holds for step i because $\Omega^{(i)}$ is
coarser than $\Omega$. Therefore, states $s$ and $s^{\prime}$ belonging to the
same aggregate in $\Omega$ will also be assigned to the same aggregate in
$\Omega^{(i+1)}$. So $\Omega^{(i+1)}$ is also coarser than $\Omega$.

The two statements above imply that $\Omega^{(k)}$ is an exactly lumpable
partition which is coarser than any other exactly lumpable partition, which
concludes the proof. The same proof can be applied for CTMCs.

The following Algorithm 3 for finding $\varepsilon$-_almost exactly lumpable_
partitions resulted from the successive refinement technique shown in the
proof above. It does not necessarily find a partition with as few aggregates
as possible, but has performed well in experiments.

The idea of the algorithm is as follows: we start with the initial partition
$\Omega=\\{S\\}$ (represented in Algorithm 3 by the aggregation function
$\omega:S\to\mathbb{N}$ which maps every state to aggregate $1$). $\Omega$ is
then successively refined. At every refinement step, for every aggregate
$\sigma\in\Omega$ and for all states $s\in\sigma$, we construct vectors of
incoming probabilities

$\displaystyle\textrm{inc}(s)=\left(\sum_{r\in\Omega_{1}}P(r,s)\,,\;\;\ldots\,,\;\;\sum_{r\in\Omega_{m}}P(r,s)\right)\in\mathbb{R}^{m}$

where $m$ is the current number of aggregates in $\Omega$. For an
$\varepsilon$-almost exactly lumpable partition, it needs to hold that the
entries of the vectors $\textrm{inc}(s)$ and $\textrm{inc}(s^{\prime})$ are
close together for $s,s^{\prime}$ in the same aggregate $\sigma$. Actually, by
18, we have that the current partition $\Omega$ is $\varepsilon$-almost
exactly lumpable if, and only if,
$\left\lVert\textrm{inc}(s)-\textrm{inc}(s^{\prime})\right\rVert_{1}\leq\varepsilon$.
If this is not the case, the algorithm therefore proceeds with the refinement
by partitioning the states in $\sigma\in\Omega$ into smaller aggregates such
that
$\left\lVert\textrm{inc}(s)-\textrm{inc}(s^{\prime})\right\rVert_{1}\leq\varepsilon$
for two states $s,s^{\prime}$ in the same aggregate in the resulting refined
partition.

This procedure stops when an $\varepsilon$-almost exactly lumpable partition
is found (at the latest when every aggregate consists of a single state). The
refinement step basically amounts to clustering points in $\mathbb{R}^{m}$
such that the maximal $\left\lVert\cdot\right\rVert_{1}$-distance between any
pair of points in a cluster is at most $\varepsilon$. To do this, we fall back
on the hierarchical clustering utilities offered by the scipy Python package,
but it is of course possible to apply different algorithms here.

Input: $\textrm{a Markov chain, defined via its transition matrix }P\textrm{
on state space }S\textrm{,}$
Input: $\textrm{and the parameter }\varepsilon\textrm{ (a generator matrix
}Q\textrm{ can be used instead of }P\textrm{)}$
Output: $\textrm{an aggregation function }\omega$
Output: $\textrm{whose corresponding partition is }\varepsilon\textrm{-almost
exactly lumpable}$

1:$\omega^{(1)}\leftarrow\left((s\in S)\mapsto 1\right)$ $\triangleright$
aggregation function

2:$i\leftarrow 1$ $\triangleright$ iteration counter

3:$m\leftarrow 1$ $\triangleright$ number of aggregates

4:repeat

5: $m_{\textrm{old}}\leftarrow m$ $\triangleright$ saves number of old
aggregates

6: $m\leftarrow 0$ $\triangleright$ counts number of new aggregates

7: for all $j\in\\{1,\ldots,m_{\textrm{old}}\\}$ do $\triangleright$ loop over
old/target aggregates

8: for all $s\in\\{r\in S:\omega^{(i)}(r)=j\\}$ do $\triangleright$ loop over
states in same aggregate

9: $\textrm{inc}(s)\leftarrow\vec{0}\in\mathbb{R}^{m_{\textrm{old}}}$

10: for all $k\in\\{1,\ldots,m_{\textrm{old}}\\}$ do $\triangleright$ loop
over potential splitter aggregates

11: $\textrm{inc}(s)_{k}\leftarrow\sum_{r\in S:\omega^{(i)}(r)=k}P(r,s)$
$\triangleright$ incoming probability from agg. $k$ to state $s$

12: end for

13: end for

14: $C\leftarrow\textrm{cluster}(\\{r\in
S:\omega^{(i)}(r)=j\\},\textrm{inc},\varepsilon)$ $\triangleright$ see below

15: for all $\sigma\in C$ do $\triangleright$ loop over clusters

16: for all $s\in\sigma$ do

17: $\omega^{(i+1)}(s)\leftarrow m+1$ $\triangleright$ states in $\sigma$ are
assigned to the same aggregate

18: end for

19: $m\leftarrow m+1$ $\triangleright$ increment aggregate number

20: end for

21: end for

22: $i\leftarrow i+1$

23:until $m_{\textrm{old}}=m$ $\triangleright$ stop when no aggregates were
split

24:return $\omega^{(i)}$

Algorithm 3 Calculating almost exactly lumpable partitions

The method $\textrm{cluster}(T,f,\varepsilon)$ takes a subset of states
$T\subseteq S$, a function $f:T\to\mathbb{R}^{k}$ and a parameter
$\varepsilon>0$ as input. The output is a partition $C$ of $T$ such that for
any cluster $\sigma\in C$ and any two states $s,s^{\prime}\in\sigma$, we have
that $\left\lVert f(s)-f(s^{\prime})\right\rVert_{1}\leq\varepsilon$. Of
course, the method should try to return as few clusters as possible, but our
Python implementation does not guarantee an optimal solution. The method
scipy.cluster.hierarchy.fclusterdata is used to calculate the clustering. We
will see that this method performs well in our experiments.

To conclude this section, we will briefly discuss the runtime of Algorithm 3.
Denote by $m$ the number of aggregates returned by the algorithm (which is not
known in advance). The outer loop (lines 4 to 23) runs through at most $m$
iterations. The loops in lines 7 to 8 lead to $n$ executions of the inner loop
on lines 10 to 12\. Lines 10 to 12, in turn, run in time $\mathcal{O}(n)$
since the loop in line 10 iterates over all aggregates, and line 11 then
calculates a sum over all states in the respective aggregate. Therefore, the
loops in lines 7 to 13 contribute a runtime of $\mathcal{O}(n^{2})$ per
iteration of the outer loop.

The runtime of line 14 depends on the clustering algorithm which is used. In
our implementation, scipy.cluster.hierarchy.fclusterdata runs in time
$\mathcal{O}(mn^{2})$ because it gets at most $n$ vectors, one per state, as
input111for the runtime, see the SciPy documentation at
https://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.hierarchy.linkage.html
and because the vectors are of dimension at most $m$. Since line 14 is
executed within the loop on line 7 (which does at most $m$ iterations), this
contributes a runtime of $\mathcal{O}(m^{2}n^{2})$ per iteration of the outer
loop (larger than the $\mathcal{O}(n^{2})$ of lines 7 to 13). Lines 15 to 20
run in $\mathcal{O}(n)$, so these are faster than line 14 and do not add to
the runtime. As the outer loop runs $m$ times, we arrive at a total runtime of
$\mathcal{O}(m^{3}n^{2})$.

Comparing with the runtimes for exact calculation of transient distributions
$p_{k}$ (respectively $p_{t}$), we see that, asymptotically for large $k$ and
$n$, applying Algorithm 3 makes sense if $m^{3}\ll k$ (respectively $m^{3}\ll
qt$ where $q$ is the maximal exit rate of all states in the CTMC). Loosely
speaking, this is usually the case if we can expect a significant state space
reduction by applying Algorithm 3. As for the SVD approach, random sampling
could potentially be used for a speed-up under some assumptions, but this is
beyond the scope of this paper.

###  Experiments

In this section, we will compare the performance of SVD sgn, SVD seba, SVD dir
(Algorithm 2), and Algorithm 3 on a selection of Markov chains. The focus is
still more on theoretical aspects than on real-world examples, which are the
subject of ongoing research. By performance comparison, we mean comparing the
error bounds given by the $\tau$ factors resulting from the aggregations
returned by the different algorithms – the lower, the better.

We first look at the setting for which the SVD algorithm variants were
designed: almost aggregatable Markov chains. These are easy to generate
randomly and we can compare the performance of the different algorithms.
Afterwards, we will see an example derived from a stochastic process algebra
model which allows for an exactly lumpable partition, so we will also see a
setting for which Algorithm 3 was designed. By default, we will calculate the
$\alpha$ distributions as in (2.5) and (2.7), and $\Pi$ (or $\Theta$) will be
set as in (2.2).

$0$$20$$40$$60$$80$$100$$120$$140$$160$$180$$200$$0$$0.2$$0.4$$0.6$$0.8$number
of aggregates$\max_{\rho\in\Omega}\tau(\rho)$SVD sgn (varying
$\varepsilon$)SVD dir / Algo. 2 (varying $\varepsilon$)SVD seba (varying
$\varepsilon$)Algo. 3 (varying $\varepsilon$) Figure 1: SVD sgn, SVD dir, SVD
seba, and Algorithm 3 executed on 100 randomly generated almost aggregatable
DTMCs with 200 states, 20 aggregates and a probability of 0.5 to have no
transition between a particular pair of aggregates. The almost aggregatable
DTMCs were obtained by random perturbation of the transition matrix of an
aggregatable DTMC with a magnitude of $0.002$. Each plotted point is an
average resulting from running the algorithms with a particular fixed input
parameter $\varepsilon$ on the 100 DTMCs.

In Figure 1, the maximal $\tau$ factor of the aggregation returned by the
algorithms (run with different input parameters) is plotted against the number
of aggregates which are found (which depends on the input parameter
$\varepsilon$). We can see that the SVD variants (except for SVD seba) perform
better than Algorithm 3 for almost aggregatable chains, which is no surprise.
In addition, the improved stability of SVD dir clearly pays off in comparison
to SVD sgn: We can see a sharp drop in the error bounds around 20 aggregates,
which was the number of aggregates in the almost aggregatable partition. For
SVD sgn, the drop is more a gradual decrease in the error bounds. Algorithm 3
does not identify the almost aggregatable partition and only reaches a similar
error bound level for around 120 aggregates.

SVD seba performs similarly to SVD dir for a low number of aggregates, but
there is a sudden change around 20 aggregates when SVD seba starts to perform
worse than all other algorithms. This is due to the fact that we limited the
maximum number of iterations of the SEBA algorithm (see [8, Algorithm 3.1], we
took the MATLAB code given in [8] and translated it into Python) to 300
iterations because of its high runtime. Regardless of the number of maximum
iterations, we could never observe SVD seba performing significantly better
than SVD dir in all our experiments. The latter is therefore a good
alternative. The better performance of SVD dir might be due to its specificity
for the given problem. The SEBA algorithm only tries to find a sparse basis
for the row space of the first $m$ rows of the matrix $V$ in the singular
value decomposition. It was designed for general applications, and does not
exploit the fact that the vectors $v(r)$ and $v(s)$ are approximate multiples
of one another for almost aggregatable DTMCs with $r$ and $s$ belonging to the
same aggregate.

$0$$100$$200$$300$$400$$500$$600$$700$$800$$900$$0$$0.5$$1$$1.5$$2$number of
aggregates$\max_{\rho\in\Omega}\tau(\rho)$SVD sgn (varying $\varepsilon$)SVD
dir / Algo. 2 (varying $\varepsilon$)Algo. 3 (varying $\varepsilon$)SVD sgn
(varying $\varepsilon$, uniform $\alpha$)SVD dir (varying $\varepsilon$,
uniform $\alpha$)Algo. 3 (varying $\varepsilon$, uniform $\alpha$) Figure 2:
SVD sgn, SVD dir, and Algorithm 3 executed on the uniformisation of the RSVP
model with $M=7$, $N=5$ and $3$ mobile nodes, resulting in a total of $842$
states. By symmetry of the mobile nodes, a lossless reduction to $234$ states
is possible.

We also considered a compositional stochastic process algebra model, the RSVP
model from [15]. It comprises a lower network channel with capacity for $M$
calls, an upper network channel with capacity for $N$ calls, and a number of
identical mobile nodes which request resources for calls at a constant rate.
Due to the mobile node symmetry in the model definition, a lossless reduction
is possible for this model. Comparing the different algorithms in Figure 2 for
a uniformisation of the RSVP model (in order to be able to apply the SVD
algorithms), we see that only Algorithm 3 identifies the partition which
results in a lossless reduction and which is exactly lumpable: the error bound
$\tau$ is equal to $0$ for $234$ aggregates. The SVD variants perform much
worse.

In addition, Figure 2 also compares the two different ways to calculate
$\alpha$: proportional $\alpha$ as in (2.5) and uniform $\alpha$ as in (2.6).
As expected, the SVD variants perform worse for uniform $\alpha$, since they
are designed to find aggregates where the $\alpha$ distributions are not
necessarily uniform. For SVD sgn, this is not visible in Figure 2, probably
due to the fact that the aggregates found in this case are anyway far from an
almost aggregatable partition and the choice of $\alpha$ does not matter much
anymore. However, we observed that uniform $\alpha$ is worse than proportional
$\alpha$ for SVD sgn in other experiments. On the other hand, Algorithm 3
actually performs better for uniform $\alpha$.

$0$$100$$200$$300$$400$$500$$600$$700$$800$$0$$0.2$$0.4$$0.6$$0.8$number of
aggregates$\max_{\rho\in\Omega}\tau(\rho)$SVD dir / Algo. 2 (varying
$\varepsilon$)Algo. 3 (varying $\varepsilon$)SVD dir / Algo. 2 (varying
$\varepsilon$, median)Algo. 3 (varying $\varepsilon$, median) Figure 3: SVD
dir, and Algorithm 3 executed on the uniformisation of the RSVP model with
$M=7$, $N=5$ and $3$ mobile nodes, resulting in a total of $842$ states.

Figure 3 compares the error bounds on the RSVP model arising when using the
median-based scheme to determine $\Pi$ as indicated in (2.3) and (2.4). Since
the median-based scheme minimizes the $\tau$ factors for fixed distributions
$\alpha$, the error bounds are lower than when using (2.2). However, we can
also see that the difference in magnitude is not too big. The additional
benefit of having an aggregated process which is also a Markov chain (because
$\Pi$ is stochastic (respectively $\Theta$ is a generator) when set as in
(2.2)) might therefore often offer a bigger advantage than the slightly
improved error bounds.

$20$$40$$60$$80$$100$$120$$140$$160$$180$$200$$220$$240$$0$$5$$10$$15$$20$number
of aggregates$\max_{\rho\in\Omega}\tau(\rho)$Algo. 3 (varying
$\varepsilon$)Algo. 3 (varying $\varepsilon$, median)Algo. 3 (varying
$\varepsilon$, uniform $\alpha$)Algo. 3 (varying $\varepsilon$, uniform
$\alpha$, median) Figure 4: Algorithm 3 executed on the CTMC arising from the
RSVP model with $M=7$, $N=5$ and $3$ mobile nodes, resulting in a total of
$842$ states.

We get a similar picture in Figure 4 if we apply Algorithm 3 directly to the
CTMC corresponding to the RSVP model: the error bounds are similar in
magnitude for the usual way of calculating $\Theta$ as in (2.2) and for the
median-based scheme. But we also see that the error bounds are much higher
than for the uniformised version. This is because the $\tau$ factors of the
uniformised chain correspond to $\frac{1}{q}$ times the $\tau$ factors of the
CTMC where $q$ is the uniformisation rate, which is set to the maximal exit
rate $30.01$ in this case. In addition, Figure 4 demonstrates again that
uniform $\alpha$ works better than proportional $\alpha$ for Algorithm 3.

To get some insights into which states are put into a common aggregate in the
RSVP model, we take a closer look at the aggregates found by Algorithm 3 in
Figure 4. If we set $\varepsilon=0.01$, then our implementation identifies the
aggregation with $234$ aggregates which is dynamic-exact. There are $68$
aggregates of size $6$ corresponding to the states for which each of the three
mobile nodes is in a distinct state, $134$ aggregates of size $3$
corresponding to the states for which two mobile nodes are in the same state,
and $32$ aggregates of size $1$ corresponding to the states for which all
mobile nodes are in the same state.

For $\varepsilon=0.1$, the identified aggregation consists of $123$ aggregates
and the error bound is still relatively low with
$\max_{\rho\in\Omega}\tau(\rho)\approx 0.057$ for uniform $\alpha$ and
$\Theta$ as in (2.2) (the median-based scheme with uniform $\alpha$ results in
$\max_{\rho\in\Omega}\tau(\rho)\approx 0.034$). We further have
$\sum_{\rho\in\Omega}\tau(\rho)\approx 1.5$. For a comparison, note that the
bound for $\sum_{\rho\in\Omega}\tau(\rho)$ given in (5.4) would yield
$123\cdot\frac{78-1}{1}\cdot 0.1=947.1$ in this case, which is more than
$10^{2}$ times higher than the actual value. The unique maximal aggregate of
size $78$ consists of a subset of the states in the original RSVP model where
one mobile node is in state $MN_{0}$ (the idle state), and the other two
mobile nodes are in states $MN_{1}$ (waiting for network resources), $MN_{3}$