2. Mathematical preliminaries

Notational conventions. We use ${\mathbb{N}}$ and ${\mathbb{R}}$ to denote the set of natural numbers and the set of real numbers. Also, we define ${\mathbb{N}_0} \,{:\!=}\, {\mathbb{N}} \cup \{0\}$ and ${{\mathbb{R}}_{+}} \,{:\!=}\, {\mathbb{R}} \setminus (-\infty,0]$ . Additionally, we denote the set $\{1,2,\ldots, N\} $ by [N]. For a set A, we denote its cardinality by | A |, and the class of all subsets of A by $2^A$ . Given $N, K \in {\mathbb{N}}$ , the set of all non-negative integer solutions to the Diophantine equation $x_1+x_2+\cdots+x_K=N$ is denoted by ${\Lambda (N,K)}$ , that is,

$$ {\Lambda (N,K)} \,{:\!=}\, \{ x = (x_1,x_2,\ldots,x_K) \in {\mathbb{N}_0}^K {\mid} x_1+x_2+\cdots+x_K=N \} .$$

We use $\mathbb{1}(.)$ to denote the indicator function. The symmetric group on a set A is denoted by ${{{\text{Sym}}} ( {A} ) }$ .

2.1. Lumpability

We first define (strong) lumpability for a discrete time Markov chain (DTMC) for ease of understanding. Standard references on this topic are [Reference Buchholz10, Reference Kemeny and Snell25, Reference Rubino and Sericola37, Reference Rubino and Sericola38].

Let $\{Y(t)\}_{t \in {\mathbb{N}}}$ be a DTMC on a state space $\mathcal{Y}= [K]$ with transition probability matrix $T = ((t_{i,j}))_{K\times K}$ , where $t_{i,j} \,{:\!=}\, {{\textsf{P}}}({ Y(2) =j {\mid} Y(1)=i })$ . Given a partition $\{ \mathcal{Y}_1, \mathcal{Y}_2,\ldots, \mathcal{Y}_M \}$ of $\mathcal{Y}$ , we define a process $\{Z(t) \}_{t \in {\mathbb{N}}}$ on [M] as follows: $Z(t) = i \in [M]$ if and only if $ Y(t) \in \mathcal{Y}_i $ , for each $t \in {\mathbb{N}}$ . The process Z is called the lumped or the aggregated process. The sets $\mathcal{Y}_i$ are often called lumping classes.

Definition 1. (Lumpability of a DTMC.) A DTMC Y on a state space $\mathcal{Y}$ is lumpable with respect to the partition $\{ \mathcal{Y}_1, \mathcal{Y}_2,\ldots, \mathcal{Y}_M \}$ of $\mathcal{Y}$ , if the lumped process Z is itself a DTMC for every choice of the initial distribution of Y [Reference Kemeny and Snell25, Chapter VI, p. 124].

A necessary and sufficient condition for lumpability, known as Dynkin’s criterion in the literature, is the following: for any two pairs of lumping classes $\mathcal{Y}_i$ and $\mathcal{Y}_j$ with $i \neq j$ , the transition probabilities of moving into $\mathcal{Y}_j$ from any two states in $\mathcal{Y}_i$ are the same, i.e. $t_{u, \mathcal{Y}_j } = t_{v, \mathcal{Y}_j }$ for all $u, v \in \mathcal{Y}_i$ , where we have used the shorthand notation $t_{u, A} = \smash{\sum_{j \in A }} t_{u,j}$ for $A \subseteq \mathcal{Y}$ . The common values, i.e. $\tilde{t}_{i,j} = t_{u, \mathcal{Y}_{j} }$ , for some $u \in \mathcal{Y}_i$ , and $ i, j \in [M]$ , form the transition probabilities of the lumped process Z. Let $\tilde{T} = (( \tilde{t}_{i,j}))_{M\times M}$ . Since Dynkin’s criterion is both necessary and sufficient, some authors alternatively define lumpability in terms of Dynkin’s criterion. In the literature, the process Z is sometimes denoted as $Z ={{{\textsf{agg}}} ( {Y} ) }$ .

Now, we move to the continuous-time case. The lumped process and the notion of lumpability for the continuous-time case are defined similarly.

Definition 2. (Lumpability of a CTMC.) A CTMC Y on a state space $\mathcal{Y}$ is lumpable with respect to the partition $\{ \mathcal{Y}_1, \mathcal{Y}_2,\ldots, \mathcal{Y}_M \}$ of $\mathcal{Y}$ , if the lumped process Z is itself a CTMC for every choice of the initial distribution of Y.

The lumpability of a CTMC can be equivalently described in terms of lumpability of a linear system of ordinary differential equations (ODEs). Consider the linear system $ \skew1\dot{y} = y A $ , where $A = ((a_{i,j}))$ is a $K\times K$ matrix (representing the transition rate matrix of the corresponding continuous-time Markov chain).

Definition 3. (Lumpability of a linear system.) The linear system $ \skew1\dot{y} = y A $ is said to be lumpable with respect to a partition $\{ \mathcal{Y}_1, \mathcal{Y}_2,\ldots, \mathcal{Y}_M \}$ of $\mathcal{Y}$ , if there exists an $M\times M$ matrix $B = ((b_{i,j}))$ satisfying Dynkin’s criterion, i.e. if $b_{i,j} = \smash{\sum_{l \in \mathcal{Y}_j }} a_{u,l} = \smash{\sum_{l \in \mathcal{Y}_j }} a_{v,l} $ for all $u, v \in \mathcal{Y}_i$ .

An alternative approach to study the lumpability of a CTMC is via its uniformization. This approach will be particularly useful when we discuss lumpability using local symmetries later. Let us consider a CTMC $\{Y(t)\}_{t \in {\mathcal{T}}}$ with transition rate matrix $A= ((a_{i,j}))$ . It is known that the original CTMC is lumpable with respect to a given partition if and only if the uniformized DTMC is lumpable with respect to the same partition. The uniformized DTMC $\tilde{Y}$ is often denoted by ${{{\textsf{unif}}} ({Y}) }$ , i.e. $ \tilde{Y} ={{{\textsf{unif}}} ({Y}) } $ . It was proved in [Reference Ganguly, Petrov and Koeppl19, Reference Rubino and Sericola38] that

\begin{align*} {{{\textsf{agg}}} ( {{{{\textsf{unif}}} ({Y}) }} ) } ={{{\textsf{unif}}} ( {{{{\textsf{agg}}} ( {Y} ) }}) } . \end{align*}

Another useful observation that will be helpful later is regarding permutation of the states. It is intuitive that permutation of elements of the state space does not destroy the lumpability property of a process. The proof of the following remark is straightforward, but is provided in Appendix A for the sake of completeness.

Remark 1. Let Y be a CTMC on $\mathcal{Y}$ with transition rate matrix $A =((a_{i,j})) $ . Let $f \in {{{\textsf{Sym}}} ( { \mathcal{Y} } ) } $ be used to permute the states. If Y (or the linear system $ \skew1\dot{y}=y A $ ) is lumpable with respect to a partition $\{ \mathcal{Y}_1, \mathcal{Y}_2,\ldots, \mathcal{Y}_M\}$ , then the process Z = f(Y) is lumpable with respect to the partition $\{ \mathcal{\tilde{Y}}_1, \mathcal{\tilde{Y}}_2,\ldots, \mathcal{\tilde{Y}}_M \} $ , where $ \mathcal{\tilde{Y}}_i = \{ f(u) {\mid} u \in \mathcal{Y}_i \} $ .

The notion of lumpability considered here is often referred to as strong lumpability in the literature. There is also a notion of weak lumpability [Reference Rubino and Sericola37, Reference Rubino and Sericola38], which we do not consider. Therefore, the term ‘lumpability’ always refers to strong lumpability in this paper. Also, note that both the DTMC and the CTMC, and the lumped processes in Definitions 1, and 3 are all assumed time-homogeneous. However, in general, it is possible to extend the notion of lumpability to allow the lumped process to be Markovian, but not necessarily time-homogeneous.

3. Markovian agent-based model

3.1. Interaction rules and the transition intensities

The most important ingredients of an MABM are the interaction rules of the agent-based local processes $X_i$ . These rules of interaction determine the dynamics of the process. Note that an MABM can also be viewed as a collection of local CTMCs that are connected to each other via the graph G. In other words, each $X_i$ can be seen as a local CTMC, conditioned on the rest. In this work, we assume the intensities of the local CTMC $X_i$ depend on the local states $X_j$ of the neighbours of the vertex $i \in V$ (such that ${(i,j) \in E}$ ). Let $d_i = |\{ j \in V {\mid} (i,j)\in E \} |$ denote the number of neighbours of vertex i. Additionally, we assume the intensities depend only on the counts of neighbours for each local state $a \in \mathcal{X}$ . Therefore, we define the following summary function c that returns population counts for different configurations of local states,

$$\begin{align*} c \colon \{ \emptyset \} \cup \bigg( \bigcup\limits_{l=1 }^{N} \mathcal{X}^l \bigg) \longrightarrow \{ \emptyset \} \cup \bigg( \bigcup\limits_{l=1 }^{N} {\Lambda (l,K)} \bigg) \end{align*}$$

such that, for $x = (x_1,x_2,\ldots, x_l) \in \mathcal{X}^l$ , and $l \in [N]$ ,

$$ \begin{align*} c( x) = ( y_1,y_2,\ldots, y_K) \in {\Lambda (l,K)} \quad \text{where } y_i = | \{ x_j = i \in \mathcal{X} \colon j =1,2, \ldots, l \} | , \end{align*} $$

and set $c(\emptyset) = \emptyset$ . Here, x is a given configuration of states of l vertices and $y_i$ is the count of vertices (within the given configuration x) that are in the ith state. The empty set $\emptyset$ is used to denote the neighbourhood of an isolated vertex. An important feature of the set-valued function c is that it is permutation-invariant in the sense that $c(x)=c(x')$ if the elements of x’ are permutations of the elements of x. In order to extract the neighbourhood information out of the global configuration, we define a family of set-valued functions $n_i$ in the following way:

$$ \begin{align*} n_i \colon \mathcal{X}^N \longrightarrow \{ \emptyset \} \cup \bigg( \bigcup\limits_{l=1 }^{N-1} \mathcal{X}^l \bigg) \quad \text{for } i \in [N], \end{align*} $$

such that, for $x=(x_1,x_2,\ldots, x_N) \in \mathcal{X}^N$ ,

$$ \begin{align*} n_i (x) = \begin{cases} (x_{i_1},x_{i_2}, \ldots, x_{i_l}) & \text{if $ (i,i_j) \in E$ for all $ j = 1,2,\ldots, l $ and $ l=d_i$,} \\ \emptyset & \text{otherwise.} \end{cases} \end{align*} $$

The neighbourhood extraction function $n_i$ for the ith vertex takes as input x, the (global) configuration of the states of all N vertices, and returns the tuple of states of the neighbours of the vertex i. Having defined these two important functions, we now define the interaction rules by means of local transition intensities. We assume the intensities depend only on the type of local transition and the summary of the neighbourhood configuration of a vertex. Therefore, we define the local intensity function

(3.1) $$ \begin{align} \gamma {\colon}\, \mathcal{X}\times \mathcal{X} \times \bigg( \{ \emptyset \} \cup \bigg( \bigcup\limits_{l=1 }^{N-1} {\Lambda (l,K)} \bigg) \bigg) \longrightarrow {{\mathbb{R}}_{+}} , \label{eqn1} \end{align} $$

where we interpret $ \gamma(a,b, y) $ as the local intensity of making a transition from local state a to b by a vertex when the summary of its neighbourhood configuration is y.

We are now in a position to specify the transition rate or the infinitesimal generator matrix for our MABM X. Note that the process X jumps from a state x to y whenever one of the local processes $X_i$ jumps. Therefore, only one of the coordinates of the states x and y differ. Let the $K^N \times K^N$ matrix $Q = ((q_{x,y}))$ denote the transition rate matrix of X. The elements of the matrix Q are given by

$$ \begin{align*} q_{x,y}= \begin{cases} \sum_{i \in [N]} {\mathbb{1}}\;(x_{i}\neq y_i,x_j=y_j \ \forall j \in V \setminus \{i\})\; \gamma (x_i, y_i, c(n_i (x)) ) & \text{if $ x \neq y$,} \\[3pt] - \sum_{y \neq x} q_{x,y} & \text{if $ x=y$.} \end{cases} \end{align*} $$

We interpret $q_{x,y}$ as the rate of transition from x to y, where $x,y \in \mathcal{X}^N$ . For ease of understanding, we have suffixed the entries of Q by the different configurations $x,y \in \mathcal{X}^N$ and interpret them as functions on $\mathcal{X}^N\times \mathcal{X}^N$ , instead of introducing a bijection between $\mathcal{X}^N$ and $[K^N]$ to label the states in a linear order so that the suffixes range over the integers from 1 to $K^N$ . Note that the particular choice of bijection to label the states is immaterial for our purposes, because such a bijection essentially yields a permutation of $[K^N]$ , and in the light of Proposition 1, does not alter lumpability properties of Q. Finally, we study the dynamics of X via the linear system

$$ \begin{align*} \dot{p} = p Q . \end{align*} $$

The vector-valued function p gives the probability distribution of X.

3.2. Examples

Susceptible-infected-susceptible (SIS) epidemics. The SIS epidemic model [Reference Simon, Taylor and Kiss40] captures the dynamics of an epidemic spread over a human or an animal population. It encapsulates binary dynamics in the sense that the local state space is written as $\mathcal{X} \,{:\!=}\, \{1,2\}$ , where 1 indicates susceptibility and 2 indicates an infected status. Infected vertices infect one of its randomly chosen neighbours at each ticking of a Poisson clock with a fixed rate a > 0. Infected vertices themselves recover to susceptibility at a rate $b\geq 0$ , independent of the neighbours’ statuses. When b = 0, the model is called a susceptible-infected (SI) model. The local transition intensities are given by

$$ \begin{align*} \gamma(1,2, (x_1,x_2)) = x_2 a \quad \text{and} \quad \gamma(2,1, (x_1,x_2)) =b . \end{align*} $$

We set $\gamma$ to zero in every other case. This fully describes the dynamics of the system.

Peer-to-peer live media streaming systems. Peer-to-peer networks are engineered networks where the vertices, called peers, communicate with each other to perform certain tasks in a distributed fashion. In particular, content delivery platforms such as BitTorrent, file-sharing platforms such as Gnutella, and (live) media (audio/video) streaming platforms use peer-to-peer networks. The main advantage of peer-to-peer systems over centralized systems is their ‘scalability’. That is, the performance of a peer-to-peer system scales well with the number of peers by virtue of the distributed task-sharing aspect, even though the global demand increases as more and more peers join the network. For the purposes of performance analysis, Markov chain models are often used for such systems.

In a peer-to-peer live-streaming system, each peer maintains a buffer of length L. The availability of a media chunk at buffer index $i\in [L]$ is indicated by 1, and likewise unavailability is indicated by 0 (see [Reference KhudaBukhsh, Rückert, Wulfheide, Hausheer and Koeppl26]). Therefore, the local state space is given by $\mathcal{X}=\{0,1\}^L$ . Set $K=2^L$ so that $\{0,1\}^L$ can be put in one-to-one correspondence with [K]. The chunk at buffer index L, if available, is played back at rate unity and then removed. After playback, all other chunks are moved one index to the right, that is, the chunk at buffer index i is shifted to buffer index $i+1$ . The central server selects a finite number of peers at random and uploads chunks at buffer index 1. All other peers (not receiving chunks from the server) download chunks from their neighbours, following a pull mechanism. Note that there are also systems where the peers push chunks into their neighbours’ buffers instead of pulling. However, we do not consider them here. The peers contact their neighbours to download missing chunks. The neighbours fulfil the request if the requested chunk is available. When multiple chunks are missing, the peers prioritize the chunks in some way, giving rise to different chunk selection strategies, such as the latest deadline first (LDF) and the earliest deadline first (EDF) strategies. Let us introduce a function, called the chunk selection function, that captures this prioritization, usually represented as probabilities. Let $s \colon [L]\times \mathcal{X} \times \mathcal{X}$ be the chunk selection function. We interpret s(i, u, v) as the probability of a vertex with buffer configuration u selecting to fill buffer index i when it contacts a neighbour with buffer configuration v.

We assume the peer-to-peer live-streaming system described above is Markovian. That is, the peers are assumed to maintain their private Poisson clocks at the tickings of which they contact their neighbours for missing chunks. Let the rate of these Poisson clocks be a > 0. Let $y_1,y_2,\ldots,y_K$ be a linear arrangement of the states in $\mathcal{X}$ . Denote the jth component of $y_i$ by $y_{i,j}$ , i.e. $y_i=( y_{i,1}, y_{i,2}, \ldots, y_{i,L})$ . The local intensity function is then given by [Reference KhudaBukhsh, Rückert, Wulfheide, Hausheer and Koeppl26],

$$ \begin{align*} \gamma(u,u+e_j, (x_1,x_2,\ldots,x_K)) = a \sum_{ i \in [K]} {\mathbb{1}}(y_{i,\:j}=1)\ x_i s( j,u,y_i) \quad \text{if } j >1 , \end{align*} $$

where $e_j$ is the jth unit vector in the L-dimensional Euclidean space, and $ (x_1,x_2,\ldots,x_K) $ is the population count vector of the neighbours of a vertex with different buffer configurations. Apart from the above transitions due to download of a chunk from a neighbour, there are two other transitions, namely, the transition due to the shifting after playback that takes place at rate unity irrespective of the buffer configurations of the neighbours, and the transition due to being directly served by the server. The latter event also takes place irrespective of the buffer configurations of the neighbours, but a rate that depends on the exact implementation set-up of the peer-to-peer system. See [Reference KhudaBukhsh, Rückert, Wulfheide, Hausheer and Koeppl26] for a detailed account.

4. Automorphism-based lumping of an MABM

Now we discuss how graph automorphisms can be used to lump states of X. The idea was introduced by Simon, Taylor, and Kiss [Reference Simon, Taylor and Kiss40] for SIS epidemics on graphs. The purpose of lumping states is to generate a Markov chain on a smaller state space. However, we should make sure that the loss of information is not too much. For instance, X is always lumpable with respect to the partition $\{\mathcal{X}^N\} $ , but if all states are lumped together, all information about the dynamics of X is lost except for the fact that total probability is conserved at all times. On the other hand, X is also lumpable with respect to the partition $ \{ \{x\} {\mid} x \in \mathcal{X}^N \} $ , which retains all the information but does not yield any state space reduction. Therefore, one needs to find a meaningful partition that yields as much state space reduction as possible with minimal loss of information. For an MABM, population counts are very useful quantities. Therefore, in order to retain information about the population counts, we first partition $\mathcal{X}^N$ into $\{ \mathcal{X}_a {\mid} a \in {\Lambda (N,K)} \}$ , that is,

(4.1) $$ \begin{align} \mathcal{X}^N = \cup_{a \in {\Lambda (N,K)} } \mathcal{X}_a \quad \text{where } \mathcal{X}_a \,{:\!=}\, \{ b \in \mathcal{X}^N {\mid} c(b)=a \} , \label{eqn2} \end{align} $$

and then seek a lumpable partition that is ideally minimally finer than this. The partition in (4.1) lumps together states that produce the same population counts. The size of this partition, i.e. $ \smash{| \{ \mathcal{X}_a {\mid} a \in {\Lambda (N,K)} \} | } $ , is $\binom{N+K-1}{K-1} $ . Note that, in the standard mean-field approach, one assumes that X is lumpable with respect to the partition in (4.1) and studies (approximate) master equations (Kolmogorov forward equations) corresponding to the different population counts. Next, we refine this partition using automorphisms.

A bijection $f \colon V \longrightarrow V$ is an automorphism on G if ${(i,j)\in E}$ if and only if ${( f(i), f( j)) \in E}$ , for all $i,j \in V$ (see [Reference Godsil and Royle21]). The collection of all automorphisms forms a group under the composition of maps. This group is denoted by ${{{\textsf{Aut}}} ( {G} ) }$ . Clearly, ${{{\textsf{Aut}}} ( {G} ) }$ is a subgroup of ${{{\textsf{Sym}}} ( { V } ) }$ . In order to use automorphisms to produce a partition of $\mathcal{X}^N$ , we shall let ${{{\textsf{Aut}}} ( {G} ) }$ act on $\mathcal{X}^N$ . We define the following group action (a map from $ {{{\textsf{Aut}}} ( {G} ) }\times \mathcal{X}^N $ to $ \mathcal{X}^N$ ):

$$ \begin{align*} f\cdot x = y \in \mathcal{X}^N \quad \text{{if and only if}} \ x_{f(i) } = y_{ i }\ \text{{for all} } i \in [N], \, \text{ for } f \in {{{\textsf{Aut}}}( {G} )},\ x \in \mathcal{X}^N . \end{align*} $$

The rationale is that, for our purpose, an automorphism needs to preserve the local states of vertices as well. Note that the action of the group ${{{\textsf{Aut}}} ( {G} ) }$ defined above can be used to introduce an equivalence relation on $ \mathcal{X}^N$ as follows: we say x and y are equivalent with respect to the action of ${{{\textsf{Aut}}} ( {G} ) }$ , denoted as $x\sim y$ , if and only if there exists an $f \in {{{\textsf{Aut}}} ( {G} ) }$ such that $ f\cdot x = y $ . The equivalence classes $\{ \tilde{\mathcal{X}}_1, \tilde{\mathcal{X}}_2, \ldots, \tilde{\mathcal{X}}_M \}$ of the relation $\sim$ yield a lumpable partition of $\mathcal{X}^N$ . Moreover, the partition thus obtained is finer than $ \{\mathcal{X}_a {\mid} a \in {\Lambda (N,K)} \}$ . We prove this in the following.

Proof. Pick any $ \tilde{\mathcal{X}}_{i} $ and $ x \in \tilde{\mathcal{X}}_{i} $ . Then $a=c(x) \in {\Lambda (N,K)}$ , and therefore $x \in \mathcal{X}_{a}$ . The proof is completed when we show that every other y in $ \tilde{\mathcal{X}}_{i} $ is also in $\mathcal{X}_{a}$ . Now $y \in \tilde{\mathcal{X}}_{i} $ implies $x \sim y$ , and therefore there exists an $f \in {{{\textsf{Aut}}} ( {G} ) }$ such that $f\cdot x=y$ . From the permutation invariance of c, we get $c( y)= c( f\cdot x ) = c(x)=a$ , implying $y \in \mathcal{X}_{a}$ .

Before proving Theorem 1, we prove the following useful lemma regarding the neighbourhood function and the action of the group ${{{\textsf{Aut}}} ( {G} ) }$ .

Lemma 1. For all $i \in [N]$ and for any $z \in \mathcal{X}^N$ , the following is true for all $f \in$ Aut(G) :

(4.2) $$ \begin{align} n_{f^{-1}(i)}( f \cdot z) = n_i(z) . \label{eqn3} \end{align} $$

Proof of Lemma 1. Let us put $ f \cdot z=x$ and $f^{-1}(i)=k $ . If $d_k=0$ , the assertion follows immediately because both sides of (4.2) are the empty set. Therefore, we assume $d_k=l>0$ . Then,

$$ \begin{align*} n_{f^{-1}(i)}( f\cdot z) & = ( x_{i_1}, x_{i_2}, \ldots, x_{i_l} )\quad \text{if } (i_j, k) \in E\ \text{{for all} } j \in [i_l] \\ & = ( z_{f(i_1)}, z_{f(i_2)}, \ldots, z_{f(i_l)} ) \quad \text{if } (i_j, k) \in E\ \text{{for all} } j \in [i_l] \\ & = n_{f(k)}(z) , \end{align*} $$

but $f(k)=i$ , implying $ n_{f^{-1}(i)}( f\cdot z) = n_i (z)$ .

Now we present the proof of Theorem 1.

Proof of Theorem 1. We check Dynkin’s criterion to establish lumpability. For any two distinct $i,j \in [M]$ , we check if $\,\vphantom{\int_{{\int}_{\int}}}\tilde{\!q}_{i,j}= \smash{\sum_{y \in \tilde{\mathcal{X}}_j }} q_{x,y}= \smash{\sum_{y \in \tilde{\mathcal{X}}_j }} q_{z,y}$ for each distinct pair $x, z \in \tilde{\mathcal{X}}_i $ . Since $z\sim x$ , there exists an $f \in {{{\textsf{Aut}}} ( {G} ) }$ such that $f\cdot z=x$ . The idea is to apply f on the states of $ \tilde{\mathcal{X}}_j $ and then show that, for any two states $x, z \in \tilde{\mathcal{X}}_i $ , there are two states $y, f\cdot y \in \tilde{\mathcal{X}}_j $ such that the neighbourhood information is preserved, that is,

$$ \begin{align*} \sum_{y \in \tilde{\mathcal{X}}_j } q_{x,y} & = \sum_{y \in \tilde{\mathcal{X}}_j } \sum_{i \in [N]} {\mathbb{1}}\;(x_{i}\neq y_i,x_j=y_j \ \forall j \neq i)\; \gamma (x_i, y_i, c(n_i (x)) ) \\ & = \sum_{f\cdot y \in \tilde{\mathcal{X}}_j } \sum_{i \in [N]} {\mathbb{1}}\;(x_{i}\neq y_{f(i)},x_j=y_{f(j)} \ \forall j \neq i)\; \gamma (x_i, y_{f(i)}, c(n_i (x)) ) \\ & = \sum_{ f\cdot y \in \tilde{\mathcal{X}}_j } \sum_{i \in [N]} {\mathbb{1}} \;(z_{f(i)}\neq y_{f(i)},z_{f(j)}=y_{f(j)} \ \forall j \neq i)\; \gamma (z_{f(i)}, y_{f(i)}, c(n_i ( f\cdot z )) ) \\ & = \sum_{ f\cdot y \in \tilde{\mathcal{X}}_j } \sum_{ f^{-1}( i) \in [N]} {\mathbb{1}} \;(z_{i}\neq y_i,z_j=y_j \ \forall j \neq i)\; \gamma (z_{i}, y_{i}, c(n_i ( z )) ) \\ & = \sum_{ y \in \tilde{\mathcal{X}}_j } \sum_{i \in [N]} {\mathbb{1}} \;(z_{i}\neq y_i,z_j=y_j \ \forall j \neq i)\; \gamma (z_{i}, y_{i}, c(n_i ( z )) )\\ & = \sum_{y \in \tilde{\mathcal{X}}_j } q_{z,y} , \end{align*} $$

where we have used $ n_{f^{-1}(i) }( f \cdot z ) = n_i (z)$ from Lemma 1. Denoting the common value by $\ \tilde{\!q}_{i,j}= \smash{\sum_{y \in \tilde{\mathcal{X}}_j }} q_{x,y} $ , the matrix $\,\tilde{\!Q} = ((\, \tilde{\!q}_{i,j} ))$ is the transition rate matrix of ${{{\textsf{agg}}} ( {X} ) }$ .

Example 2. (Star graph.) An automorphism on a star graph leaves the central node (root) unchanged and permutes the rest of the nodes (leaf nodes) in any possible manner. Without loss of generality, let us assume the central node is labelled N. Then, the automorphism group ${{{\textsf{Aut}}} ( {G} ) }$ is given by

$$ \begin{align*} &{{{\textsf{Aut}}} ( {G} ) } =\\ &\quad \{g \in {{{\textsf{Sym}}} ( {[N]} ) } {\mid} g(N)= N, g(i)=f(i)\ \text{{for all} } i\in [N-1], \text{ for some } f \in {{{\textsf{Sym}}} ( {[N-1]} ) } \} . \end{align*} $$

5. Lumping states using local symmetry

In this section, we discuss lumping ideas based on a local notion of automorphism. In many cases, the number of automorphisms decrease drastically as the graph grows arbitrarily large. For instance, it is known that Erdős–Rényi random graphs tend to be asymmetric with probability approaching unity as the size of the graph N grows to infinity [Reference Łuczak31]. Similar statements are true for d-regular random graphs under various sets of conditions on d relative to the number of vertices N [Reference Kim, Sudakov and Vu27], and random graphs with specified degree distributions [Reference McKay and Wormald32]. As a consequence, the automorphism-based lumping tends to be ineffective in state space reduction as the size of the graphs grows arbitrarily. Therefore, it is desirable to bring in a notion of local automorphism or local symmetry that would allow swapping vertices that are locally indistinguishable (i.e. have similar neighbourhoods), but are not so globally. This notion of symmetry is weaker than an automorphism, which endows global symmetry on a graph. However, the potential gain is in the ability to engender state space reduction when the graph grows arbitrarily large, rendering automorphism-based lumping virtually ineffective. In the following, we make these ideas precise.

5.1. Local symmetry

There have been several attempts to formulate a more flexible notion of local symmetry. However, the literature seems divided on this and there is not a single universally accepted concept. In our set-up, it seems intuitive that two vertices that are locally indistinguishable in a large graph would also behave indistinguishably and could therefore be swapped. A notion of local symmetry identifying such vertices was proposed in [Reference Elbert Simões, Figueiredo and Barbosa16], which we adopt in this paper. We need a few definitions to make precise what we mean by two vertices being locally indistinguishable.

In order to define locality, we need some notion of distance between vertices of G. Let d(u, v) denote the smallest distance (length of the minimal path) between two vertices $u,v \in V$ . If u, and v are not connected, i.e. there is no path between them, we simply set $d(u,v)=\infty$ .

Definition 4. Given a vertex $u \in V$ , define its k-neighbourhood in G, denoted by $N_{k}(u)$ , as follows:

$$ \begin{align*} N_{k} \colon V \longrightarrow 2^V \quad \text{such that } N_{k}(u) \,{:\!=}\, \{ v \in V {\mid} d(u,v) \leq k\} . \end{align*} $$

Let $G[ N_{k}(u)]$ denote the subgraph of G induced by $ N_{k}(u) $ . The notion of locality we adopt in this paper hinges on these k-neighbourhoods and their induced subgraphs. If two vertices induce isomorphic subgraphs, they are indistinguishable locally and we say they are k-locally symmetric [Reference Elbert Simões, Figueiredo and Barbosa16].

Definition 5. Two vertices $u,v\in V$ are defined to be k-locally symmetric if there exists an isomorphism f between $G[N_k(u)]$ and $G[N_k(v)]$ such that f(u) = v.

Therefore, two vertices $u,v \in V$ are k-locally symmetric if their kth-order local structures (k-hop neighbourhoods) are equivalent in the sense that there is a structure-preserving (edge-preserving in this case) bijection between them. When k = 1, we simply say the vertices are locally symmetric.

As with automorphisms, local symmetries also induce an equivalence relation on the set of vertices V. We say two vertices $u,v\in V$ are equivalent with respect to k-local symmetry, denoted by $u \smash{\overset{k} \sim} v$ , if there exists an isomorphism f between $G[N_k(u)]$ and $G[N_k(v)]$ such that f(u) = v. The notion of local symmetry is related to the concept of views in the discrete mathematics literature [Reference Hendrickx23, Reference Yamashita and Kameda43]. The view of depth k of a vertex is a tree containing all walks of length k leaving that vertex. However, note that, in our context, the induced subgraphs $G[N_k(u)]$ need not be trees. The following facts about local symmetry are useful for our study of lumpability [Reference Elbert Simões, Figueiredo and Barbosa16, Reference Norris35].

Proposition 2. The following properties are satisfied by k-local symmetry.

(i) For $u,v \in V$ , $ u {\overset{k+1} \sim} v \Rightarrow u {\overset{k} \sim} v $ . Consequently, ${V}/\smash{\overset{k+1} \sim},$ the equivalence classes of $\smash{\overset{k+1} \sim}$ , form a refinement of ${V}/\smash{\overset{k} \sim}$ , the equivalence classes of $\smash{\overset{k} \sim}$ .

(ii) If the equivalence classes of $ \smash{\overset{k+1} \sim}$ are the same as those of $ \smash{\overset{k} \sim}$ , the equivalence classes of all $ {\overset{k+j} \sim}$ are the same as those of $ \smash{\overset{k} \sim}$ , for $j \in {\mathbb{N}}$ .

(iii) If $k \geq \textrm{diam}(G)$ , the diameter of G, then, for two vertices $u,v \in V$ , we have $ u \smash{\overset{k} \sim} v $ if and only if there exists an $f \in$ Aut(G) such that f(u) = v. That is, k-local symmetry is equivalent to automorphism if k is at least as large as the diameter of G.

In addition to the above, it can be verified that the quotient spaces ${V}/\smash{\overset{k} \sim}$ are equitable partitions [Reference Godsil and Royle21, Chapter 9] for each $k\geq 1$ . We use these properties to lump states of $\mathcal{X}^N$ in the next section.

5.2. Lumping states using local symmetry

The procedure to lump states in $\mathcal{X}^N$ using local symmetry is similar to the procedure used to lump states using automorphism. However, unlike the case with automorphism, we now allow permutations that only need to ensure symmetry locally. That is, in order to lump states using k-local symmetry, we allow permuting two vertices u and v in V if and only if u and v are k-locally symmetric. Therefore, define

$$ \begin{align*} \Psi_k (G) \,{:\!=}\, \{ f \in {{{\textsf{Sym}}} ( {V} ) } {\mid} f(u)=v \ \text{{if and only if}} \ u \smash{\overset{k}\sim} v, \text{ for } u, v \in V \} . \end{align*} $$

We refer to $| \Psi_k (G) |$ as the number of k-local symmetries. It can be verified that $ \Psi_k (G)$ , for each $k \geq 1$ , forms a group under the composition of maps. Therefore, we can let the group $ \Psi_k(G) $ act on $\mathcal{X}^N$ . We define the action of $ \Psi_k(G) $ as follows:

$$ \begin{align*} f\cdot x = y \in \mathcal{X}^N \quad \text{{if and only if}} \ x_{f(i) } = y_{ i }\ \text{{for all} } i \in [N]{,} \text{ for } f \in \Psi_k (G) , x \in \mathcal{X}^N . \end{align*} $$

Note that a state x in $\mathcal{X}^N$ is taken to y if and only if the local states of all vertices are preserved and two vertices are swapped only when they are k-locally symmetric. The above action induces the following partition of the state space: two states $ x,y \in \mathcal{X}^N$ are said to be equivalent with respect to k-local symmetry, denoted as $x \smash{\overset{k}\sim} y$ , if there exists an $f \in \Psi_k(G)$ such that $f\cdot x=y$ . We use the same symbol $\smash{\overset{k}\sim}$ since there is no scope for confusion. The equivalence classes of ${\smash{\overset{k}\sim}} $ are obtained, as before, by determining the orbits of states in $\mathcal{X}^N$ . The orbit of a state $x \in \mathcal{X}^N$ is given by $ \Psi_k(G) \cdot x \,{:\!=}\, \{ f\cdot x \in \mathcal{X}^N {\mid} f \in \Psi_k (G) \} $ .

The partition thus obtained (based on k-local symmetry) does not, in general, guarantee lumpability, i.e. the X need to be lumpable with respect to $\mathcal{X}^N/\smash{\overset{k}\sim}$ . Therefore, we construct another Markov chain that approximates the lumped process and seek to quantify the approximation error in the next section. The following observation is integral to the quantification of the approximation error incurred when states of $\mathcal{X}^N$ are lumped according to k-local symmetry instead of automorphism.

Proposition 3. The quotient space $\mathcal{X}^N/\smash{\overset{k+1}\sim}$ is a refinement of $\mathcal{X}^N/\smash{\overset{k}\sim}$ .

Proof of Proposition 3. Let $ \mathcal{X}_{1}^{(k+1)} , \mathcal{X}_{2}^{(k+1)} ,\ldots, \mathcal{X}_{M_{k+1}}^{(k+1)} $ be the equivalence classes of $\smash{\overset{k+1} \sim}$ . Also, denote the equivalence classes of $\smash{\overset{k} \sim}$ by $ \mathcal{X}_{1}^{(k)} , \mathcal{X}_{2}^{(k)} ,\ldots, \mathcal{X}_{M_k}^{(k)} $ . Let $ i \in [M_{k+1}]$ and $x \in \mathcal{X}_{i}^{(k+1)} $ . If $ \mathcal{X}_{i}^{(k+1)} $ is a singleton, the identity map is the only map in $ \Psi_{k+1} (G) $ , but it is also in $\Psi_{k} (G)$ . Therefore, $x \in \mathcal{X}_{j}^{(k)} $ for some $j \in [M_k]$ , and the assertion follows. If $ \mathcal{X}_{i}^{(k+1)} $ has at least two elements, say, x, y, then $y \smash{\overset{k+1}\sim} x$ . By Proposition 2, we must have $y \smash{\overset{k}\sim} x$ . Therefore, there exists a $j \in [M_k]$ such that $x,y \in \mathcal{X}_{j}^{(k)} $ . Since the choice of x, y is arbitrary, the assertion follows.

For practical applications, one would start with $\mathcal{X}^N/ \smash{\overset{1}\sim}$ and then iteratively obtain further refinements $\,\mathcal{X}^N\,{/}\, \smash{\overset{2}\sim}, \mathcal{X}^N/ \smash{\overset{3}\sim} $ , and so on until satisfactory accuracy is achieved (assuming we can quantify accuracy for the time being). In the light of Proposition 2, two important remarks are in order. They emphasize the benefits of local symmetry-driven lumping over automorphism-driven lumping.

Remark 3. In an algorithmic implementation, item (ii) in Proposition 2 provides a stopping rule for an iterative procedure to obtain local symmetry-driven partitions. That is, we can stop at the first instance of no improvement (the equivalence classes of $\smash{\overset{k+1}\sim}$ and $\smash{\overset{k}\sim}$ are the same). For the sake of illustration, consider the chain graph in Figure 1. For this trivial example, notice that going from k = 1 to k = 2 does not cause any improvement. Therefore, we can stop at k = 1 even though the diameter is 4.

Figure 1: A chain graph with $1\smash{\overset{1}\sim} 4$ and $2\smash{\overset{1}\sim} 3$ as well as $1\smash{\overset{2}\sim} 4$ and $2\smash{\overset{2}\sim} 3$ .

Remark 4. The diameters in many random graphs grow slowly as the number of vertices goes to infinity. For instance, the diameter of Erdős–Rényi random graphs with N vertices and edge probability $\lambda/N$ , for some fixed $\lambda>1$ , grows as log N (see [Reference Riordan and Wormald36]). A similar result holds for Newman, Strogatz, and Watts graphs too (see [Reference Chatterjee and Durrett12, Lemma 1.2]). In the light of item (iii) of Proposition 2, our approach needs (at most) as many steps as the diameter of G to produce an exactly lumpable partition of $\mathcal{X}^N$ . Note that $k \geq \textrm{diam}(G)$ is only a sufficient condition for $\mathcal{X}^N/ \smash{\overset{k}\sim}$ to be an exactly lumpable partition. For practical purposes, we actually achieve sufficient accuracy (including exact lumpability) even for small values of $k < \textrm{diam}(G) $ . Supporting numerical results are provided later in Section 7.

In the next example, we show the amount of state space reduction that can be achieved by means of local symmetry-driven lumping.

Example 5. Consider an SIS process on the graph with 20 vertices in Figure 2. The number of states in this case is $2^{20}=1048576$ . However, there are plenty of local symmetries that we can exploit to reduce the size of the state space. In Table 1, we summarize the state space reduction achieved for increasing k. The compression level (used to quantify the reduction achieved) is calculated as one minus the ratio of the size of the reduced state space to the size of the original state space, i.e. 1048576.

Figure 2: A graph with 20 vertices for Example 5.

Table 1. State space reduction with k.

Our local symmetry-driven lumping approach shares a close relationship with what are known as fibrations in algebraic graph theory. We briefly describe the relationship in the following.

6. Graph fibrations

Fibrations of graphs were first inspired by fibrations between a pair of categories [Reference Boldi and Vigna9]. Although the idea of fibrations originated from category theory, it has deep implications for graph theory, theoretical computer science, and other mathematical disciplines. For instance, Boldi, Lonati, Santini, and Vigna [Reference Boldi, Lonati, Santini and Vigna8] discussed its interesting connections to the PageRank citation ranking algorithm. Nijholt, Rink, and Sanders [Reference Nijholt, Rink and Sanders33] explored the similarities between dynamical systems with a network structure and dynamical systems with symmetry by means of fibrations of graphs. Let us now define the necessary graph-theoretic concepts.

Given the graph G = (V, E), we first define the source and target maps ${{\textsf{s}}}_G, {{\textsf{t}}}_G \colon E \longrightarrow V$ on G such that ${{\textsf{s}}}_G (u,v) =u$ and ${{\textsf{t}}}_G(u,v)=v$ for each $(u,v)\in E$ . Let $H =(V', E')$ be another graph. The source and the target maps ${{\textsf{s}}}_H, {{\textsf{t}}}_H$ are defined analogously. A map $f \,{:\!=}\, ( f_v, f_e)$ , where $f_v\,:\,V \longrightarrow V'$ and $f_e \colon E \longrightarrow E'$ , is called a graph morphism between G and H (from G to H, to be precise) if $f_v$ and $f_e$ commute with the source and the target maps of G and H, i.e. if ${{\textsf{s}}}_H\, f_e = f_v {{\textsf{s}}}_G $ and $ {{\textsf{t}}}_H\, f_e = f_v {{\textsf{t}}}_G $ . A morphism is called an epimorphism if both $f_v$ and $f_e$ are surjective. Finally, we define a graph fibration as follows [Reference Boldi and Vigna9].

In the original paper [Reference Boldi and Vigna9], Boldi and Vigna defined colour-preserving graph morphisms when graphs are endowed with a colouring function. In that case, $f_e$ also commutes with the colouring function. For our present purposes, we do not require this generality and only consider uncoloured graphs. Boldi and Vigna [Reference Boldi and Vigna9] showed that a left action of a group on G can be used to induce fibrations. They also showed that fibrations and epimorphisms satisfying certain local in-isomorphism property are equivalent [Reference Boldi and Vigna9, Theorem 2]. Indeed, fibrations have a close relationship with the notion of local symmetry described in Section 5. The proof of the following proposition follows analogously from [Reference Boldi and Vigna9, Theorem 2]. However, for the sake of completeness, we also provide it in appendix A.

The above proposition essentially shows that the equivalence classes of local symmetry (with k = 1) and fibres induced by a graph fibration are the same. Therefore, the fibres can also be used to aggregate the states of $\mathcal{X}^N$ to achieve approximate lumpability in the same fashion as we did with local symmetry. See Figure 3.

Figure 3: Fibrations map vertices to vertices and edges to edges. When three vertices form a triangle, fibrations also preserve the triangle structure. Therefore, one can define an isomorphism between local neighbourhoods using fibrations.

7. Approximation error

As the lumping based on local symmetry does not ensure Markovianness of the lumped process, we need to quantify the approximation error. In order to do so, we work with the uniformization of X. Then we lump ${{{\textsf{unif}}} ( {X}) }$ to produce ${{{\textsf{agg}}} ( {{{{\textsf{unif}}} ( {X}) }} ) }$ according to k-local symmetry. A direct assessment of the quality of aggregation is cumbersome. Therefore, it is suggested [Reference Deng, Mehta and Meyn15, 20] that we lift the aggregated process ${{{\textsf{agg}}} ( {{{{\textsf{unif}}} ( {X}) }} ) }$ to a Markov chain on the same state space $\mathcal{X}^N$ as $ {{{\textsf{unif}}} ( {X}) }$ and then compare their transition probability matrices. The lifting allows us to use known metrics of divergence such as the Kullback–Leibler divergence to quantify the approximation error. We follow the scheme depicted in Figure 4.

Figure 4: Lifting procedure used to assess the quality of the approximation.

In order to fix ideas, let us lump $ {{{\textsf{unif}}} ( {X}) }$ according to k-local symmetry, i.e. according to the partition $\{ \mathcal{X}_{1}^{(k)} , \mathcal{X}_{2}^{(k)} ,\ldots, \mathcal{X}_{M_k}^{(k)} \} $ of $\mathcal{X}^N$ obtained as the equivalence classes of $\smash{\overset{k}\sim}$ . We introduce two notations in this connection. Let $\eta_k \colon \mathcal{X}^N \longrightarrow [M_k] $ be the partition function associated with $\smash{\overset{k}\sim}$ , i.e. $\eta_k(x)\,{:\!=}\, i $ if and only if $ x \in \mathcal{X}_{i}^{(k)} $ . For $u\in \mathcal{X}^N$ , let us denote the equivalence class containing u by $\langle {x} \rangle_{k}$ , i.e. $\langle {x} \rangle_{k} \,{:\!=}\, \mathcal{X}_{i}^{(k)} $ if and only if $ x \in \mathcal{X}_{i}^{(k)} $ . Note that $\langle {x} \rangle_{k} = \eta_k^{-1}(\eta_k (x)) $ .

Let $T =((t_{i,j}))$ be the transition probability matrix associated with $ {{{\textsf{unif}}} ( {X}) }$ . Now, since X is not necessarily lumpable with respect to the partition $\{ \mathcal{X}_{1}^{(k)} , \mathcal{X}_{2}^{(k)} ,\ldots, \mathcal{X}_{M_k}^{(k)} \} $ , for $i \neq j \in [M_k] $ and two distinct $\smash{x,y \in \mathcal{X}_{i}^{(k)}}$ , the quantity ${\sum_{z \in \mathcal{X}_{j}^{(k)}}} t_{x,z}$ may not equal ${\sum_{z \in \mathcal{X}_{j}^{(k)}}} t_{y,z}$ . In other words, ${{{\textsf{agg}}} ( { {{{\textsf{unif}}} ( {X}) }} ) }$ may not be Markovian. The idea is to approximate ${{{\textsf{agg}}} ( { {{{\textsf{unif}}} ( {X}) }} ) }$ with a Markov chain and then quantify the approximation error. Therefore, in order to avoid additional notations, we proceed as if ${{{\textsf{agg}}} ( { {{{\textsf{unif}}} ( {X}) }} ) }$ were Markovian and estimate its transition probabilities. If ${{{\textsf{unif}}} ( {X}) }$ is stationary with distribution $\pi$ , i.e. if $\pi$ is the solution to $\pi T= \pi$ and $p(0)=\pi$ , a natural estimate of the transition probability of the lumped process is as follows:

$$ \begin{align*} \tilde{t}_{i,j}^{(k)} \,{:\!=}\, \dfrac{ \sum_{u \in \mathcal{X}_{i}^{(k)} } \pi_u \sum_{v \in \mathcal{X}_{j}^{(k)} } t_{u,v} }{ \sum_{u \in \mathcal{X}_{i}^{(k)} } \pi_u } \quad \text{for } i, j \in [M_k] . \end{align*} $$

That is, we estimate the transition probabilities of the (Markov chain that approximates the) lumped process ${{{\textsf{agg}}} ( { {{{\textsf{unif}}} ( {X}) }} ) }$ by averaging the different values $\smash{\sum_{z \in \mathcal{X}_{j}^{(k)}}} t_{x,z}$ and $\smash{\sum_{z \in \mathcal{X}_{j}^{(k)}}} t_{y,z}$ , weighted by the stationary probabilities [Reference Geiger, Petrov, Kubin and Koeppl20]. Let $\tilde{T}^{(k)} \,{:\!=}\, (( \tilde{t}_{i,j}^{(k)} ))$ . Before we proceed to explain the lifting procedure, let us consider an example.

Figure 5: (a) Two graphs, and (b) comparisons of the expected number of infected nodes for SIS dynamics on the graphs.

Now, we describe how the transition probabilities of the lifted Markov chain are calculated. There are two common ways of lifting ${{{\textsf{agg}}} ( {{{{\textsf{unif}}} ( {X}) }} ) }$ to a Markov chain on $\mathcal{X}^N$ , one using a probability vector, called $\pi$ -lifting, and the other using the transition probabilities, called P-lifting. Let us discuss $\pi$ -lifting first.

Definition 7. ( $\pi$ -lifting.) The $\pi$ -lifting of $\eta_k({{{\textsf{unif}}} ( {X}) } )$ is a DTMC with transition probability matrix $T^{\pi}_{k} \,{:\!=}\, ((t_{u,v}^{\pi,k} ))$ , where

$$ \begin{align*} t_{u,v}^{\pi,k} \,{:\!=}\, \dfrac{ \pi_v }{ \sum_{ x \in \langle {v} \rangle_{k} } \pi_x } \tilde{t}_{ \eta_k(u) ,\eta_k (v)}^{(k)} , \quad \text{where } u, v \in \mathcal{X}^N . \end{align*} $$

Note that, in principle, $\pi$ -lifting can be done using any probability vector as long as the denominator remains non-zero for the choice of the candidate probability vector. Nevertheless, the most common choice is the stationary probability vector. The reason for this choice is the fact that the stationary probability vector achieves the minimum KL divergence rate [Reference Deng, Mehta and Meyn15]. For this reason, in this paper we consider $\pi$ -lifting with the stationary distribution for numerical computations. Another immediate consequence of $\pi$ -lifting is that the lifted Markov chain with transition probability matrix $T_k^{\pi}$ given in (7) is lumpable with respect to the partition $\mathcal{X}^N\!/\smash{\overset{k}\sim}$ and has $\pi$ as the stationary probability. In many situations, computing the stationary distribution $\pi$ may not be straightforward. Many numerical methods can be adopted in such situations. See Section 8 for a more detailed discussion. Now, we define the approximation error.

Definition 8. We define the approximation error to be the KL divergence rate between ${{{\textsf{unif}}} ( {X}) }$ and the lifted DTMCs. Therefore, for $\pi$ -lifting, the approximation error is given by

$$ \begin{align*} {D_\textrm{KL} ( T {\mid\!\!\!\!\mid} T^{\pi}_{k} ) } \,{:\!=}\, \sum_{u \in \mathcal{X}^N } \sum_{v \in \mathcal{X}^N } \pi_u t_{u,v} \log \bigg( \dfrac{ t_{u,v} }{ t_{u,v}^{\pi,k} } \bigg) . \end{align*} $$

Having defined the approximation error, we show that it indeed decreases monotonically with the order of local symmetry. This is precisely the assertion of Theorem 2. However, in order to prove Theorem 2, we need to make use of the following calculation, which we present in the form of a lemma.

Lemma 2. For any two states $u,v \in \mathcal{X}^N$ , and for any k, define the ratio

$$ \begin{align*} \rho_k (u,v) \,{:\!=}\, \dfrac{ \sum_{t \in \langle {v} \rangle_{k}} \pi_t }{ \tilde{t}_{\eta_k(u), \eta_k(v) }^{(k)} } = \dfrac{ \sum_{p \in \langle {u} \rangle_{k}} \pi_p \sum_{q \in \langle {v} \rangle_{k} } \pi_q }{ \sum_{p \in \langle {u} \rangle_{k} } \sum_{q \in \langle {v} \rangle_{k}} \pi_p t_{p,q} } . \end{align*} $$

Then, the following recursion relation holds:

$$ \begin{align*} \sum_{x \in \langle {u} \rangle_{k} } \sum_{y \in \langle {v} \rangle_{k} } \pi_x t_{x,y} \rho_{k+1}(x,y) & = \rho_k(u,v) \sum_{x \in \langle {u} \rangle_{k}} \sum_{y \in \langle {v} \rangle_{k} } \pi_x t_{x,y} . \end{align*} $$

Proof of Lemma 2. By the refinement property of local symmetry in Proposition 3, we can find distinct integers $i_1,i_2,\ldots, i_m$ and $j_1,j_2, \ldots, j_n$ in $[M_{k+1}]$ such that

(7.1) $$ \begin{align} \langle {u} \rangle_{k} =\cup_{ l \in [m] } \mathcal{X}_{ i_l }^{(k+1)} \quad\text{and}\quad \langle {v} \rangle_{k} =\cup_{ l \in [n] } \mathcal{X}_{ j_l }^{(k+1)} . \label{eqn4} \end{align} $$

Therefore, we can split the summation over $ \smash{\langle {u} \rangle_{k}} $ , $ \langle {v} \rangle_{k} $ into disjoint equivalence classes of $\overset{k+1}\sim$ . Within each of these equivalence classes of $\overset{\kern-1ptk+1\phantom{0}}\sim,$ the quantity $\tilde{t}_{ \eta_{k+1}(x), \eta_{k+1}( y) }$ is constant and can therefore be pulled out of the summation. Therefore,

$$ \begin{align*} &\sum_{x \in \langle {u} \rangle_{k} } \sum_{y \in \langle {v} \rangle_{k} } \pi_x t_{x,y} \rho_{k+1}(x,y) \\ & \qquad = \sum_{p \in [m]} \sum_{q \in [n]} \sum_{x \in \mathcal{X}_{ i_p }^{(k+1)} } \sum_{y \in \mathcal{X}_{ j_q }^{(k+1)} } \pi_x t_{x,y} \bigg( \dfrac{ \sum_{s \in \langle {x} \rangle_{k+1}} \pi_s \sum_{t \in \langle {y} \rangle_{k+1}} \pi_t }{ \sum_{s \in \langle {x} \rangle_{k+1} } \sum_{t \in \langle {y} \rangle_{k+1} } \pi_s t_{s,t} } \bigg) \\ & \qquad = \sum_{p \in [m]} \sum_{q \in [n]} \bigg( \dfrac{ \sum_{s \in \mathcal{X}_{ i_p }^{(k+1)} } \pi_s \sum_{t \in \mathcal{X}_{ j_q }^{(k+1)} } \pi_t }{ \sum_{s \in \mathcal{X}_{ i_p }^{(k+1)} } \sum_{t \in \mathcal{X}_{ j_q }^{(k+1)} } \pi_s t_{s,t} } \bigg) \sum_{x \in \mathcal{X}_{ i_p }^{(k+1)} } \sum_{y \in \mathcal{X}_{ j_q }^{(k+1)} } \pi_x t_{x,y}\\ & \qquad = \sum_{x \in \langle {u} \rangle_{k} } \sum_{y \in \langle {v} \rangle_{k} } \pi_x \pi_y \\ & \qquad = \rho_k(u,v) \sum_{x \in \langle {u} \rangle_{k} } \sum_{y \in \langle {v} \rangle_{k} } \pi_x t_{x,y} . \end{align*} $$

This completes the proof.

Note that $\rho_k(u,v)=\rho_k(x,y)$ for any $u \smash{\overset{k}\sim} x$ and $v \smash{\overset{k}\sim} y$ . Therefore, we can use the shorthand notation $\rho_k( \mathcal{X}_i^{(k)} , \mathcal{X}_j^{(k)} ) $ to mean $\rho_k(u,v)$ for any $u \in \mathcal{X}_i^{(k)} , v \in \mathcal{X}_j^{(k)} $ .

Remark 5. (Averaging argument.) The main implication of Lemma 2 is that the quantity $\rho_k(u,v) $ can be seen as a weighted average of $\rho_{k+1} (x,y)$ where the x, y are in the equivalence classes of $\smash{\overset{k+1}\sim}$ . The weights are precisely

$$ \begin{align*} W_{ \langle {u} \rangle_{k}, {\langle {v} \rangle_{k}} } ( \mathcal{X}_{ i_p }^{(k+1)} , \mathcal{X}_{ j_q }^{(k+1)} ) \,{:\!=}\, \dfrac{ \sum_{x \in \mathcal{X}_{ i_p }^{(k+1)} } \sum_{y \in \mathcal{X}_{ j_q }^{(k+1)} } \pi_x t_{x,y} }{ \sum_{x \in {\langle {u} \rangle_{k}} } \sum_{y \in {\langle {v} \rangle_{k}} } \pi_x t_{x,y}} , \end{align*} $$

where we have partitioned ${\langle {u}\rangle_{k}}$ and ${\langle {v} \rangle_{k}}$ into $ \mathcal{X}_{ i_p }^{(k+1)} $ and $ \mathcal{X}_{ j_q }^{(k+1)} $ respectively, as shown in (7.1). We interpret $W_{ {\langle {u}\rangle_{k}}, {\langle {v} \rangle_{k}} } ( \mathcal{X}_{ i_p }^{(k+1)} , \mathcal{X}_{ j_q }^{(k+1)} ) $ as the weight for the cross-section $ \mathcal{X}_{ i_p }^{(k+1)} \times \mathcal{X}_{ j_q }^{(k+1)} $ with regard to the partition of ${\langle {u}\rangle_{k}}$ and ${\langle {y} \rangle_{k}}$ given in (7.1). Therefore, it follows from Lemma 2 that

$$ \begin{align*} \rho_k({\langle {u}\rangle_{k}} ,{\langle {v} \rangle_{k}}) = \sum_{p \in [m]} \sum_{q \in [n]} \rho_{k+1}( \mathcal{X}_{ i_p }^{(k+1)} , \mathcal{X}_{ j_q }^{(k+1)} ) W_{ {\langle {u}\rangle_{k}}, {\langle {v} \rangle_{k}} } ( \mathcal{X}_{ i_p }^{(k+1)} , \mathcal{X}_{ j_q }^{(k+1)} ) . \end{align*} $$

Since the weights sum up to unity, $\rho_k(u,v) $ can indeed be seen as an average. Keeping this remark in mind, we now proceed to state and prove Theorem 2 about the monotonicity of the approximation error.

Theorem 2. For $\pi$ -lifting, the aggregation of states in $\mathcal{X}^N$ using local symmetry ensures monotonically decreasing approximation error with increasing order of local symmetry, that is,

$$ \begin{align*} {D_\textrm{KL}( {T} {\mid\!\!\!\!\mid} {T_{k+1}^{(\pi)} } ) } \leq {D_\textrm{KL}( {T} {\mid\!\!\!\!\mid} {T_{k}^{(\pi)}} ) } \quad \text{for all } k \geq 1 . \end{align*} $$

Proof of Theorem 2. By the refinement property of local symmetry proved in Proposition 3, we can partition $[M_{k+1}] =\{1,2,\ldots,M_{k+1}\} $ into $\{\Lambda_1, \Lambda_2, \ldots, \Lambda_{M_k} \}$ such that

$$ \begin{align*} \mathcal{X}_{i}^{(k)} = \cup_{l \in \Lambda_i } \mathcal{X}_{l }^{k+1} . \end{align*} $$

Note that

$$ \begin{align*} &{D_\textrm{KL} ( {T} {\mid\!\!\!\!\mid} {T_{k}^{(\pi)}} ) }- {D_\textrm{KL}( {T} {\mid\!\!\!\!\mid} {T_{k+1}^{(\pi)} } ) } \\ & \qquad = \sum_{i, j \in [M_k] } \sum_{u \in \mathcal{X}_{i}^{(k)} } \sum_{v \in \mathcal{X}_{j}^{(k)} } \pi_u t_{u,v} \log \bigg( \dfrac{ \rho_k(u,v) }{ \rho_{k+1}(u,v) } \bigg) \\ & \qquad = \sum_{i, j \in [M_k] } (\!\log ( \rho_k(\mathcal{X}_{i}^{(k)},\mathcal{X}_{j}^{(k)} ) ) \sum_{u \in \mathcal{X}_{i}^{(k)} } \sum_{v \in \mathcal{X}_{j}^{(k)} } \pi_u t_{u,v} - \sum_{u \in \mathcal{X}_{i}^{(k)} } \sum_{v \in \mathcal{X}_{j}^{(k)} } \pi_u t_{u,v} \log ( \rho_{k+1}(u,v) ) ) \\ & \qquad = \sum_{i, j \in [M_k] } \Theta_{i,j} , \end{align*} $$

where

$$ \begin{align*} \Theta_{i,j} & \,{:\!=}\, \log ( \rho_k(\mathcal{X}_{i}^{(k)},\mathcal{X}_{j}^{(k)} ) ) \sum_{u \in \mathcal{X}_{i}^{(k)} } \sum_{v \in \mathcal{X}_{j}^{(k)} } \pi_u t_{u,v} \\ & \quad\, - \sum_{p \in \Lambda_i} \sum_{q \in \Lambda_j} \sum_{u \in \mathcal{X}_{p}^{(k+1)} } \sum_{v \in \mathcal{X}_{q}^{(k+1)} } \pi_u t_{u,v} \log ( \rho_{k+1}(u,v) ) \\ & = \bigg( \sum_{u \in \mathcal{X}_{i}^{(k)} } \sum_{v \in \mathcal{X}_{j}^{(k)} } \pi_u t_{u,v} \bigg) \times (\!\log ( \rho_k(\mathcal{X}_{i}^{(k)},\mathcal{X}_{j}^{(k)} ) ) \\ &\quad\, - \sum_{p \in \Lambda_i} \sum_{q \in \Lambda_j} W_{ \mathcal{X}_{i}^{(k)},\mathcal{X}_{j}^{(k)} }(\mathcal{X}_{p}^{(k+1)},\mathcal{X}_{q}^{(k+1)} ) \log ( \rho_{k+1}(\mathcal{X}_{p}^{(k+1)},\mathcal{X}_{q}^{(k+1)} ) ) ) \\ & \geq 0 , \end{align*} $$

by Jensen’s inequality and the averaging argument given in Remark 5 and Lemma 2. This completes the proof.

Note that ${D_\textrm{KL} ( {T} {\mid\!\!\!\!\mid} {T_{k}^{(\pi)}} ) }- {D_\textrm{KL} ( {T} {\mid\!\!\!\!\mid} {T_{k+1}^{(\pi)} } ) }=0$ is achieved if (and only if) equality occurs in Jensen’s inequality, forcing the individual $\Theta_{i,j}$ to be zeros. This is the case when the $\rho_k$ and $\rho_{k+1}$ are equal. There are two possibilities. First, the equivalence classes of $\smash{\overset{k}{\sim}}$ and $\smash{\overset{k+1}{\sim}}$ are the same. In this case, by Proposition 2, the equivalence classes of all $\smash{{\overset{k+j}{\sim}}}$ , for $j\geq 2$ , will remain the same. Therefore, we have already reached automorphism, and hence exact lumpability. Second, the equivalence classes of $\smash{\overset{k}{\sim}}$ and $\smash{\overset{k+1}{\sim}}$ are different (so we are not yet at automorphism), but exact lumpability has already been achieved at the order of local symmetry k. In both cases, we need not refine our partition further because exact lumpability has been achieved. Therefore, ${D_\textrm{KL}( {T} {\mid\!\!\!\!\mid} {T_{k}^{(\pi)}} ) }- {D_\textrm{KL}( {T} {\mid\!\!\!\!\mid} {T_{k+1}^{(\pi)} } ) }=0$ serves as a definite stopping rule for any iterative algorithmic implementation of local symmetry-driven lumping.

Now, we discuss the second type of lifting, which makes use of the transition probabilities and is called P-lifting. The following is the definition.

Definition 9. (P-lifting.) The P-lifting of $\eta_k({{{\textsf{unif}}}( {X})} )$ is a DTMC with transition probability matrix $T^{P}_{k} \,{:\!=}\, ((t_{u,v}^{P,k} ))$ where, for $u, v \in \mathcal{X}^N$ ,

$$ \begin{align*} t_{u,v}^{P,k} \,{:\!=}\, \begin{cases} \dfrac{ t_{u,v} }{ \sum_{ x \in {\langle {v} \rangle_{k}} } t_{u,x} } \tilde{t}_{ \eta_k(u) ,\eta_k (v)}^{(k)} & \text{if } \sum_{ x \in {\langle {v} \rangle_{k}} } t_{u,x} >0 , \\\\[-6pt] \dfrac{1 }{ | {\langle {v} \rangle_{k}} | } \tilde{t}_{ \eta_k(u) ,\eta_k (v)}^{(k)} & \text{if } \sum_{ x \in {\langle {v} \rangle_{k}} } t_{u,x} =0 . \end{cases} \end{align*} $$

The approximation error for P-lifting is defined similarly. Note that P-lifting is sharp, in the sense that if the lumping is in fact exact, then ${D_\textrm{KL} ( { T } {\mid\!\!\!\!\mid} { T_{k}^{(P)} } ) } =0$ , whereas $\pi\!$ -lifting is not [Reference Geiger, Petrov, Kubin and Koeppl20]. In Figure 6, we show numerical results pertaining to Theorem 2. We consider the Barabási–Albert preferential attachment, the Erdős–Rényi and the Watts–Strogatz small-world random graphs. The claimed monotonicity is observed in all three cases. In fact, the KL divergence rate decreases steeply in all three cases, for both $\pi\!$ - and P-lifting. The figures are particularly encouraging in that a satisfactory level of accuracy is achieved even for small orders of local symmetry. Since one of the main purposes of aggregation is to engender state space reduction, we need to evaluate the performance of local symmetry-driven aggregation in terms of some notion of compression level as well. Therefore, we define compression level C at the order of local symmetry k as follows:

$$ \begin{align*} {C}(k)=1- \dfrac{M_k}{| \mathcal{X}^N | } , \end{align*} $$

Figure 6: Monotonicity of the KL divergence with the order of the local symmetry for the SIS dynamics on different models of random graphs with 10 vertices. All graphs are undirected. The Erdős–Rényi graphs are created with edge probability 0.3, while the Watts–Strogatz small-world networks are created with rewiring probability 0.3 and each vertex having three neighbours. The infection and recovery rates are both 0.5.

where $M_k$ is the cardinality of the quotient space $ \mathcal{X}^N/ \smash{\overset{k}\sim} $ , i.e. the number of equivalence classes of $\smash{\overset{k}\sim}$ . If there is no non-trivial local symmetries, the compression level is zero because the partition is simply $\{ \{x\} {\mid} x \in \mathcal{X}^N \}$ . In Figure 7 we show how the number of local symmetries decreases drastically as we increase the order of local symmetry. Consequently, the compression level also falls steeply. This is expected because random graphs tend to become asymmetric as the number of vertices increases.

Figure 7: As we increase the order, the number of local symmetries (the cardinality of $\Psi_k$ ) decreases drastically. Therefore, the compression level also decreases. The simulation set-up is the same as in Figure 6.

8. Discussions

The idea of using Markov chain lumpability for model reduction has been discussed in the literature for some years now. For instance, Simon, Taylor, Kiss, and Miller [Reference Kiss, Miller and Simon28, Reference Simon and Kiss39, Reference Simon, Taylor and Kiss40] considered epidemiological scenarios, focusing mainly on binary dynamics. More general Markovian agent-based models were considered in [Reference Banisch6]. Lumpability abstractions were applied to rule-based systems in [Reference Feret, Henzinger, Koeppl and Petrov17] from a theoretical computer science perspective. While model reduction is one of the main purposes of lumpability, it is not the only one. In a recent paper [Reference Katehakis and Smit24], Katehakis and Smit identified a class of Markov chains, which they call successively lumpable and for which the stationary probabilities can be computed successively by computing stationary probabilities of a cleverly constructed sequence of Markov chains (typically on much smaller state spaces).