# On Kemeny’s constant and stochastic complement111This article was partially
funded by the “INdAM – GNCS Project: Metodi basati su matrici e tensori
strutturati per problemi di algebra lineare di grandi dimensioni” code
CUP_E53C22001930001, and by the PRIN “Low-rank Structures and Numerical
Methods in Matrix and Tensor Computations and their Application” code
20227PCCKZ. S. Kim is supported in part by funding from the Fields Institute
for Research in Mathematical Sciences and from the Natural Sciences and
Engineering Research Council of Canada. The second and fourth authors are
member of the INdAM GNCS group.

Dario Andrea Bini<EMAIL_ADDRESS>Fabio Durastante
<EMAIL_ADDRESS>Sooyeong Kim<EMAIL_ADDRESS>Beatrice Meini
<EMAIL_ADDRESS>
###### Abstract

Given a stochastic matrix $P$ partitioned in four blocks $P_{ij}$, $i,j=1,2$,
Kemeny’s constant $\kappa(P)$ is expressed in terms of Kemeny’s constants of
the stochastic complements $P_{1}=P_{11}+P_{12}(I-P_{22})^{-1}P_{21}$, and
$P_{2}=P_{22}+P_{21}(I-P_{11})^{-1}P_{12}$. Specific cases concerning periodic
Markov chains and Kronecker products of stochastic matrices are investigated.
Bounds to Kemeny’s constant of perturbed matrices are given. Relying on these
theoretical results, a divide-and-conquer algorithm for the efficient
computation of Kemeny’s constant of graphs is designed. Numerical experiments
performed on real world problems show the high efficiency and reliability of
this algorithm.

###### keywords:

Markov chains , Kemeny’s constant , divide-and-conquer algorithm

###### MSC:

[2010] 60J22 , 65C40 , 65F15

††journal: Linear Algebra and its Applications

[pisa]organization=Mathematics Department, University of
Pisa,addressline=Largo Bruno Pontecorvo, 5, city=Pisa, postcode=56127,
state=PI, country=Italy [york]organization=Department of Mathematics and
Statistics, York University,addressline=4700 Keele Street, city=Toronto,
postcode=M3J 1P3, state=ON, country=CA

Expression of Kemeny’s constant employing the constants of stochastic
complements.

New recursive algorithms for computing Kemeny’s constant

##  1 Introduction

Given an ergodic, finite, discrete-time, time-homogeneous Markov chain,
Kemeny’s constant is the expected time for the Markov chain to travel between
randomly chosen states, where these states are sampled according to the
stationary distribution. Originally defined in [22] as the expected time to
reach a randomly-chosen state from a fixed starting state, this quantity is
independent of the choice of initial state. Various intuitive explanations for
the constancy have been provided in [3, 14, 27]. Moreover, Kemeny’s constant
has many applications to a variety of subjects, including the study of road
traffic networks [1, 11], disease spread [24, 36], and many others.

In addition, Kemeny’s constant has recently received significant attention
within the graph theory community. Such constant finds special relevance in
the study of random walks on graphs. In the context of random walks, Kemeny’s
constant is used as a measure of connectivity of the graph: the smaller the
constant, the faster a random walker moves around the graph. As a graph
invariant, much research has been dedicated to understanding how the structure
of a graph influences Kemeny’s constant. Extremal Kemeny’s constant has been
studied in [5, 10, 19]. The impact of edge addition/removal on this quantity
has also been explored in [1, 9, 16, 23, 28, 29]. Recently, work in [25]
provides insights into the interplay between a graph and its complement
regarding Kemeny’s constant.

The pursuit of reducing the computational complexity of Kemeny’s constant of a
Markov chain with $n$ states, which is $O(n^{3})$, is interesting. Randomized
approaches for the direct approximation of Kemeny’s constant have been used in
[30, 35]. Additionally, exploring computational questions where partial
information is available without computing from scratch has been investigated.
Work in [1] investigated computation of Kemeny’s constant for a graph obtained
by removing an edge. Moreover, work from [6] provided an explicit formula for
Kemeny’s constant for graphs with bridges in terms of some quantities on the
subgraphs resulting from the deletion of the bridges.

In this article, we furnish a formula for Kemeny’s constant of a Markov chain
by utilizing two Markov chains induced from the original, known as censored
(watched) Markov chains (see [31]). Let $P$ be the transition matrix of an
irreducible Markov chain with state space $\\{1,\ldots,n\\}$. Let

$P=\begin{bmatrix}P_{11}&P_{12}\\\ P_{21}&P_{22}\end{bmatrix},$ (1)

where $P_{11}$ and $P_{22}$ are square of size $m\times m$ and
$(n-m)\times(n-m)$, respectively. The transition matrices $P_{1}$ and $P_{2}$
of the censored Markov chains are given by

$\begin{split}P_{1}=P_{11}+P_{12}(I-P_{22})^{-1}P_{21},\quad\quad
P_{2}=P_{22}+P_{21}(I-P_{11})^{-1}P_{12}.\end{split}$ (2)

In particular, we will provide expressions of the difference
$\gamma=\kappa(P)-\kappa(P_{1})-\kappa(P_{2})$, given in terms of the
stationary distribution vector of the Markov chain $P$.

The paper is organized as follows: In Section 2, we introduce some notation
and review the definitions and properties of stochastic complements and
Kemeny’s constant. Section 3 concerns the main theoretical result (Theorem
3.3), where Kemeny’s constant of a Markov chain is related to Kemeny’s
constants of censored Markov chains. Section 4 explores applications of the
main result for various structured transition matrices and subsequently
determines the minimum Kemeny’s constant of a periodic Markov chain. In
Section 5, we establish lower and upper bounds for the constant $\gamma$ and
bounds on Kemeny’s constant of perturbed matrices are given. In Section 6, we
introduce a divide-and-conquer algorithm to compute $\kappa(P)$ where $P$ is a
large and highly sparse matrix, and we present numerical experiments performed
with real world matrices to show the effectiveness and reliability of this
algorithm. Section 7 draws the conclusions.

##  2 Notation and preliminaries

A matrix $A=[a_{ij}]\in\mathbb{R}^{m\times n}$ is said to be nonnegative
(resp. positive), if $a_{ij}\geq 0$ ($a_{ij}>0$) for any $i,j$, and we write
$A\geq 0$ (resp. $A>0$). We denote by $\operatorname{\mathbf{1}}_{n}$ the all
ones vector of size $n$. If $n$ is clear from the context, we shall omit the
subscript of $\operatorname{\mathbf{1}}_{n}$. A matrix $A$ is said to be
reducible if there exists a permutation matrix $\Pi$ such that $\Pi A\Pi^{T}$
is a block upper triangular matrix with square diagonal blocks. If $A$ is not
reducible, then we say that $A$ is irreducible. A nonnegative matrix $A$ is
said to be stochastic if
$A\operatorname{\mathbf{1}}=\operatorname{\mathbf{1}}$. For
$A\in\mathbb{R}^{n\times n}$, we use $\rho(A)$ to denote the spectral radius
of $A$.

A matrix $A$ is a _non-singular (resp. singular)_ M-matrix if $A=sI-B$ with
$B\geq 0$ and $s>\rho(B)$ (resp. $s=\rho(B)$). It is known that if $A$ is a
non-singular M-matrix, then $a_{ii}>0$, $a_{ij}\leq 0$ for $i\neq j$, and
$A^{-1}\geq 0$.

Given a vector $\mathbf{v}=(v_{i})\in\mathbb{R}^{n}$, we denote by
$\|\mathbf{v}\|$ the 1-norm of $\mathbf{v}$, i.e.,
$\|\mathbf{v}\|=\sum_{i=1}^{n}|v_{i}|$, and we use $\|\mathbf{v}\|_{\infty}$
to indicate the infinity norm of $\mathbf{v}$, that is,
$\|\mathbf{v}\|_{\infty}=\max_{i}|v_{i}|$. For a matrix $A$, we denote by
$\|A\|_{\infty}$ the norm induced by the infinity norm, that is,
$\|A\|_{\infty}=\max_{j}\sum_{i=1}^{n}|a_{ij}|$.

###  2.1 Stochastic complement

Let $P$ be an $n\times n$ irreducible stochastic matrix, with stationary
distribution vector $\boldsymbol{\pi}$, i.e.,
$\boldsymbol{\pi}^{T}P=\boldsymbol{\pi}^{T}$ and
$\boldsymbol{\pi}^{T}\operatorname{\mathbf{1}}=1$. Let $P$ be partitioned into
the $2\times 2$ block matrix as in (1). Since $P$ is irreducible, $I-P_{11}$
and $I-P_{22}$ are non-singular M-matrices. It is found in [31] that the
matrices $P_{1}$ and $P_{2}$ in (2) are stochastic and irreducible. We call
$P_{1}$ (resp. $P_{2}$) the _stochastic complement_ of $P_{11}$ (resp.
$P_{22}$) in $P$. The matrix $P_{1}$ represents the transition matrix of the
Markov chain obtained by censoring the states $\\{m+1,\ldots,n\\}$, while
$P_{2}$ represents the transition matrix of the Markov chain obtained by
censoring the states $\\{1,\ldots,m\\}$. This suggests the name “censored
Markov chain”.

Let $\boldsymbol{\pi}$ be partitioned conformally with $P$ so that
$\boldsymbol{\pi}^{T}=\begin{bmatrix}\boldsymbol{\pi}_{1}^{T}&\boldsymbol{\pi}_{2}^{T}\end{bmatrix}$
where $\boldsymbol{\pi}_{1}\in\mathbb{R}^{m}$ and
$\boldsymbol{\pi}_{2}\in\mathbb{R}^{n-m}$. Denoting by
$\hat{\boldsymbol{\pi}}_{i}$ the stationary distribution vector of $P_{i}$, we
have

$\hat{}\boldsymbol{\pi}_{i}=\frac{1}{\|\boldsymbol{\pi}_{i}\|}\boldsymbol{\pi}_{i},~{}~{}i=1,2.$

Moreover, setting $\alpha_{i}=\|\boldsymbol{\pi}_{i}\|$, we see that

$\boldsymbol{\pi}^{T}=[\alpha_{1}\hat{}\boldsymbol{\pi}_{1}^{T},\alpha_{2}\hat{}\boldsymbol{\pi}_{2}^{T}],$
(3)

and the vector $\boldsymbol{\alpha}^{T}=[\alpha_{1},\alpha_{2}]$ satisfies

$\boldsymbol{\alpha}^{T}S=\boldsymbol{\alpha}^{T},~{}~{}\boldsymbol{\alpha}^{T}\operatorname{\mathbf{1}}=1,\quad
S=\begin{bmatrix}\hat{}\boldsymbol{\pi}_{1}^{T}P_{11}\operatorname{\mathbf{1}}_{m}&\hat{}\boldsymbol{\pi}_{1}^{T}P_{12}\operatorname{\mathbf{1}}_{n-m}\\\
\hat{}\boldsymbol{\pi}_{2}^{T}P_{21}\operatorname{\mathbf{1}}_{m}&\hat{}\boldsymbol{\pi}_{2}^{T}P_{22}\operatorname{\mathbf{1}}_{n-m}\end{bmatrix}.$
(4)

That is, $\boldsymbol{\alpha}$ is the stationary distribution vector of the
aggregated matrix $S$, which is stochastic and irreducible.

###  2.2 Kemeny’s constant

Let $P$ be the transition matrix of an ergodic, finite, discrete-time, time-
homogeneous Markov chain. _Kemeny’s constant_ of $P$, denoted by $\kappa(P)$,
is the expected number of time steps required by the Markov chain to go from a
given state $i$ to a random state $j$, sampled according to the stationary
distribution $\boldsymbol{\pi}$. Kemeny’s constant has an algebraic
characterization, in terms of the trace of a suitable matrix [21, 34], as
stated in the following.

######  Lemma 2.1.

Let $\mathbf{g},\mathbf{h}\in\mathbb{R}^{n}$ be vectors with
$\mathbf{h}^{T}\mathbf{g}=1$, $\mathbf{h}^{T}\operatorname{\mathbf{1}}\neq 0$,
$\boldsymbol{\pi}^{T}\mathbf{g}\neq 0$. Then, $I-P+\mathbf{g}\mathbf{h}^{T}$
is non-singular, and

$\kappa(P)=\operatorname{Tr}(Z)-\boldsymbol{\pi}^{T}Z\operatorname{\mathbf{1}}.$
(5)

where $Z=(I-P+\mathbf{g}\mathbf{h}^{T})^{-1}$.

By choosing $\mathbf{g}=\operatorname{\mathbf{1}}$, we find that

$\kappa(P)=\operatorname{Tr}(Z)-1\quad\text{where}\quad
Z=(I-P+\operatorname{\mathbf{1}}\mathbf{h}^{T})^{-1}.$ (6)

Let $\lambda_{1}=1,\lambda_{2},\ldots,\lambda_{n}$ be the eigenvalues of $P$.
Since $P$ is irreducible, $\lambda_{1}$ is simple. By the Brauer theorem [4],
the eigenvalues $1-\lambda_{i}$ of the matrix $I-P$ coincide with the
eigenvalues of $I-P+\operatorname{\mathbf{1}}\mathbf{h}^{T}$ except for the
eigenvalue $1-\lambda_{1}=0$ which is mapped into 1. This fact enable us to
express Kemeny’s constant in terms of the eigenvalues $\lambda_{i}$ of $P$ as

$\kappa(P)=\sum_{i=2}^{n}\frac{1}{1-\lambda_{i}}.$ (7)

Recall that given two square stochastic matrices $A,B$ of the same size, the
products $AB$ and $BA$ are stochastic as well and have the same characteristic
polynomial, see [17, Theorem 1.3.20]. So, their eigenvalues coincide and have
the same algebraic multiplicities. Thus, it follows from (7) that

$\kappa(AB)=\kappa(BA).$ (8)

A slightly different situation is encountered for non-square matrices. Suppose
that $A$ and $B$ are stochastic matrices of size $m\times n$ and $n\times m$,
respectively, where $n>m$. Then $AB$ and $BA$ are square stochastic matrices
of size $m\times m$ and $n\times n$, respectively. Moreover, their nonzero
eigenvalues coincide, thus we deduce from (7) that

$\kappa(BA)=\kappa(AB)+n-m.$ (9)

###  2.3 Kemeny’s constant and graphs

A graph is a pair $G=(V,E)$, where $V$ is a finite set of _vertices_ of
cardinality $|V|=n$, and $E$ is a set of pairs $\\{i,j\\}$ where $i,j\in V$,
called _edges_. A graph is _weighted_ if it is equipped with a function
$a:E\to\mathbb{R}_{+}$ that associates a non-negative _weight_ to each edge.
For unweighted graphs, the weight of each edge is assumed to be equal to 1.
The _adjacency matrix_ of a graph on $n$ vertices is the $n\times n$ matrix
$A=[a_{ij}]$ such that $a_{ij}$ is equal to the weight of edge $\\{i,j\\}$ if
$\\{i,j\\}\in E$, and zero otherwise.

Given a graph $G$, we can associate the Markov chain having transition matrix
$P=D^{-1}A$, where $D=\textrm{diag}(\mathbf{d})$,
$\mathbf{d}=A\operatorname{\mathbf{1}}$, and $A$ is the adjacency matrix of
the graph. Here, we assume that $d_{i}\neq 0$ for any $i$, that is there are
no isolated vertices. The matrix $P$ describes a random walk on the graph $G$.
We can then define Kemeny’s constant of the graph as $\kappa(G):=\kappa(P)$.

A graph is said _bipartite_ if there exist disjoint sets of vertices $V_{1}$
and $V_{2}$, of cardinality $m_{1}$ and $m_{2}$, respectively, such that
$V=V_{1}\cup V_{2}$, and the vertices in $V_{1}$ as well the vertices in
$V_{2}$ are connected by no edges. The adjacency matrix $H$ of a bipartite
graph can be partitioned into 4 blocks $H_{i,j}$, $i,j=1,2$ where $H_{11}$ and
$H_{22}$ are the null square matrices of size $m_{1}$, and $m_{2}$,
respectively. If the graph is unweighted complete bipartite then the blocks
$H_{12}$ and $H_{21}$ of size $m_{1}\times m_{2}$, $m_{2}\times m_{1}$,
respectively, have all the entries equal to 1.

##  3 Kemeny’s constant and stochastic complement

In this section we relate Kemeny’s constant of the transition matrix $P$ in
(1) to Kemeny’s constants of the stochastic complements (2). We begin with
introducing the following preliminary result.

######  Lemma 3.1.

Let $P$ be the transition matrix given by (1), and let
$Z=(I-P+\operatorname{\mathbf{1}}\mathbf{h}^{T})^{-1}$ where
$\mathbf{h}\in\mathbb{R}^{n}$ with
$\mathbf{h}^{T}\operatorname{\mathbf{1}}=1$. Partition $\mathbf{h}$
conformally with $P$ as
$\mathbf{h}^{T}=\begin{bmatrix}\mathbf{h}_{1}^{T}&\mathbf{h}_{2}^{T}\end{bmatrix}$.
Then

$\operatorname{Tr}(Z)=\operatorname{Tr}(Z_{1})+\operatorname{Tr}(Z_{2})$ (10)

where

$\begin{split}Z_{1}=(I-P_{1}+\sigma_{2}^{-1}\mathbf{u}_{1}\mathbf{v}_{1}^{T})^{-1},\quad\quad
Z_{2}=(I-P_{2}+\sigma_{1}^{-1}\mathbf{u}_{2}\mathbf{v}_{2}^{T})^{-1},\end{split}$
(11)

and

$\begin{split}&\mathbf{u}_{1}=\operatorname{\mathbf{1}}_{m}+P_{12}(I-P_{22})^{-1}\operatorname{\mathbf{1}}_{n-m},~{}~{}\mathbf{u}_{2}=\operatorname{\mathbf{1}}_{n-m}+P_{21}(I-P_{11})^{-1}\operatorname{\mathbf{1}}_{m},\\\
&\mathbf{v}_{1}^{T}=\mathbf{h}_{1}^{T}+\mathbf{h}_{2}^{T}(I-P_{22})^{-1}P_{21},~{}~{}\mathbf{v}_{2}^{T}=\mathbf{h}_{2}^{T}+\mathbf{h}_{1}^{T}(I-P_{11})^{-1}P_{12},\\\
&\sigma_{i}=1+\mathbf{h}_{i}^{T}(I-P_{ii})^{-1}\operatorname{\mathbf{1}}\quad\text{for}\quad
i=1,2.\end{split}$

###### Proof.

By partitioning $Z$ conformally with $P$ in (1), we find that

$Z=\begin{bmatrix}I-P_{11}+\operatorname{\mathbf{1}}_{m}\mathbf{h}_{1}^{T}&-P_{12}+\operatorname{\mathbf{1}}_{m}\mathbf{h}_{2}^{T}\\\
-P_{21}+\operatorname{\mathbf{1}}_{n-m}\mathbf{h}_{1}^{T}&I-P_{22}+\operatorname{\mathbf{1}}_{n-m}\mathbf{h}_{2}^{T}\end{bmatrix}^{-1}.$

We shall omit the subscripts of all ones vectors. Then

$Z=\begin{bmatrix}Z_{1}&*\\\ *&Z_{2}\end{bmatrix},$

where

$\begin{split}&Z_{1}=\left(I-P_{11}+\operatorname{\mathbf{1}}\mathbf{h}_{1}^{T}-(P_{12}-\operatorname{\mathbf{1}}\mathbf{h}_{2}^{T})(I-P_{22}+\operatorname{\mathbf{1}}\mathbf{h}_{2}^{T})^{-1}(P_{21}-\operatorname{\mathbf{1}}\mathbf{h}_{1}^{T})\right)^{-1},\\\
&Z_{2}=\left(I-P_{22}+\operatorname{\mathbf{1}}\mathbf{h}_{2}^{T}-(P_{21}-\operatorname{\mathbf{1}}\mathbf{h}_{1}^{T})(I-P_{11}+\operatorname{\mathbf{1}}\mathbf{h}_{1}^{T})^{-1}(P_{12}-\operatorname{\mathbf{1}}\mathbf{h}_{1}^{T})\right)^{-1},\end{split}$
(12)

and $*$ denotes a generic entry. In particular,
$\operatorname{Tr}(Z)=\operatorname{Tr}(Z_{1})+\operatorname{Tr}(Z_{2})$. By
using the Sherman–Morrison formula, we find that

$(I-P_{22}+\operatorname{\mathbf{1}}\mathbf{h}_{2}^{T})^{-1}=(I-P_{22})^{-1}-\sigma_{2}^{-1}(I-P_{22})^{-1}\operatorname{\mathbf{1}}\mathbf{h}_{2}^{T}(I-P_{22})^{-1},$

where
$\sigma_{2}=1+\mathbf{h}_{2}^{T}(I-P_{22})^{-1}\operatorname{\mathbf{1}}$.
Replacing this expression in the first formula in (12) yields the formula for
$Z_{1}$ in (11). We proceed similarly for $Z_{2}$. ∎

Observe that the vectors $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$ in Lemma 3.1
depend on the vector $\mathbf{h}$ that can be chosen arbitrarily under the
condition $\mathbf{h}^{T}\operatorname{\mathbf{1}}=1$. In the next proposition
we show that $\mathbf{h}$ can be chosen with
$\mathbf{v}_{1}=\hat{}\boldsymbol{\pi}_{1}$ and
$\mathbf{v}_{2}=\hat{}\boldsymbol{\pi}_{2}$. This choice will allow us to
express Kemeny’s constant of $P$ in terms of Kemeny’s constant of the
stochastic complements $P_{1}$ and $P_{2}$, which will be shown in Theorem
3.3.

######  Proposition 3.2.

There exist vectors $\mathbf{h}_{1}\in\mathbb{R}^{m}$ and
$\mathbf{h}_{2}\in\mathbb{R}^{n-m}$ such that

$\begin{split}&\mathbf{h}_{1}^{T}+\mathbf{h}_{2}^{T}(I-P_{22})^{-1}P_{21}=\hat{}\boldsymbol{\pi}_{1}^{T},\\\
&\mathbf{h}_{2}^{T}+\mathbf{h}_{1}^{T}(I-P_{11})^{-1}P_{12}=\hat{}\boldsymbol{\pi}_{2}^{T},\end{split}$
(13)

and
$\mathbf{h}_{1}^{T}\operatorname{\mathbf{1}}_{m}+\mathbf{h}_{2}^{T}\operatorname{\mathbf{1}}_{n-m}=1$.

###### Proof.

Observe that (13) is equivalent to

$\begin{split}&\mathbf{h}_{1}^{T}\left(I-(I-P_{11})^{-1}P_{12}(I-P_{22})^{-1}P_{21}\right)=\hat{}\boldsymbol{\pi}_{1}^{T}-\hat{}\boldsymbol{\pi}_{2}^{T}(I-P_{22})^{-1}P_{21},\\\
&\mathbf{h}_{2}^{T}=\hat{}\boldsymbol{\pi}_{2}^{T}-\mathbf{h}_{1}^{T}(I-P_{11})^{-1}P_{12}.\end{split}$
(14)

Since the matrix $T=(I-P_{11})^{-1}P_{12}(I-P_{22})^{-1}P_{21}$ is stochastic,
the vector $\operatorname{\mathbf{1}}_{m}$ is orthogonal to the rows of $I-T$.
On the other hand, since $(I-P_{22})^{-1}P_{21}$ is stochastic and
$\hat{}\boldsymbol{\pi}^{T}_{1}\operatorname{\mathbf{1}}_{m}=\hat{}\boldsymbol{\pi}_{2}^{T}\operatorname{\mathbf{1}}_{n-m}=1$,
the vector $\operatorname{\mathbf{1}}_{m}$ is orthogonal to the right hand
side of the first equation in (14). Hence, the system has solution. If
$\mathbf{h}_{1}$ is a solution, then $\mathbf{h}_{2}$ can be recovered from
the second equation in (14). By multiplying to the right by
$\operatorname{\mathbf{1}}_{m}$ the first equation in (13), we find that
$\mathbf{h}_{1}^{T}\operatorname{\mathbf{1}}_{m}+\mathbf{h}_{2}^{T}\operatorname{\mathbf{1}}_{n-m}=1$.
∎

Proposition 3.2, together with Equation (10) and Lemma 2.1, implies the
following result.

######  Theorem 3.3.

Let $P$ be a stochastic irreducible matrix, partitioned as in (1), and denote
by $P_{1}$ and $P_{2}$ the stochastic complements (2). Then

$\kappa(P)=\kappa(P_{1})+\kappa(P_{2})+\gamma,$ (15)

where

$\gamma=\begin{bmatrix}\widehat{}\boldsymbol{\pi}_{1}^{T}&-\widehat{}\boldsymbol{\pi}_{2}^{T}\end{bmatrix}(I-P+\operatorname{\mathbf{1}}\boldsymbol{\pi}^{T})^{-1}\begin{bmatrix}\|\boldsymbol{\pi}_{2}\|\operatorname{\mathbf{1}}\\\
-\|\boldsymbol{\pi}_{1}\|\operatorname{\mathbf{1}}\end{bmatrix}.$

###### Proof.

From (6), $\kappa(P)=\operatorname{Tr}(Z)-1$, where
$Z=(I-P+\operatorname{\mathbf{1}}\mathbf{h}^{T})^{-1}$ and $\mathbf{h}$ is a
vector such that $\mathbf{h}^{T}\operatorname{\mathbf{1}}=1$. Partition
$\mathbf{h}$ as $\mathbf{h}^{T}=[\mathbf{h}_{1}^{T},\mathbf{h}_{2}^{T}]$,
where $\mathbf{h}_{1}$ and $\mathbf{h}_{2}$ satisfy (13). We shall maintain
the same notation in Lemma 3.1. Then
$\mathbf{v}_{1}^{T}=\hat{}\boldsymbol{\pi}_{1}^{T}$ and
$\mathbf{v}_{2}^{T}=\hat{}\boldsymbol{\pi}_{2}^{T}$. Since
$\hat{}\boldsymbol{\pi}_{1}^{T}(I-P_{1})=0$, the eigenvalues of
$I-P_{1}+\sigma_{2}^{-1}\mathbf{u}_{1}\hat{}\boldsymbol{\pi}_{1}^{T}$ are the
eigenvalues of $I-P_{1}$, except for the eigenvalue equal to zero, which is
replaced by $\sigma_{2}^{-1}\hat{}\boldsymbol{\pi}_{1}^{T}\mathbf{u}_{1}$.
Similarly, the eigenvalues of
$I-P_{2}+\sigma_{1}^{-1}\mathbf{u}_{2}\hat{}\boldsymbol{\pi}_{2}^{T}$ are the
eigenvalues of $I-P_{2}$, except for the eigenvalue equal to zero, which is
replaced by $\sigma_{1}^{-1}\hat{}\boldsymbol{\pi}_{2}^{T}\mathbf{u}_{2}$. In
view of (10),
$\operatorname{Tr}(Z)=\operatorname{Tr}(Z_{1})+\operatorname{Tr}(Z_{2})$. From
the expression of $\mathbf{u}_{1}$ in Lemma 3.1, we find that

$\hat{}\boldsymbol{\pi}_{1}^{T}\mathbf{u}_{1}=1+\hat{}\boldsymbol{\pi}_{1}^{T}P_{12}(I-P_{22})^{-1}\operatorname{\mathbf{1}}_{n-m}=1+\frac{\boldsymbol{\pi}_{2}^{T}\operatorname{\mathbf{1}}_{n-m}}{\|\boldsymbol{\pi}_{1}\|}=1+\frac{\|\boldsymbol{\pi}_{2}\|}{\|\boldsymbol{\pi}_{1}\|}=\frac{1}{\|\boldsymbol{\pi}_{1}\|}.$

Hence,
$\operatorname{Tr}(Z_{1})=\kappa(P_{1})+\sigma_{2}{\|\boldsymbol{\pi}_{1}\|}$.
Similarly,
$\hat{}\boldsymbol{\pi}_{2}^{T}\mathbf{u}_{2}=\frac{1}{\|\boldsymbol{\pi}_{2}\|}$
and so
$\operatorname{Tr}(Z_{2})=\kappa(P_{2})+\sigma_{1}{\|\boldsymbol{\pi}_{2}\|}$.
Since
$\operatorname{Tr}(Z)=\kappa(P_{1})+\kappa(P_{2})+{\sigma_{2}}{\|\boldsymbol{\pi}_{1}\|}+{\sigma_{1}}{\|\boldsymbol{\pi}_{2}\|}$,
we have $\kappa(P)=\kappa(P_{1})+\kappa(P_{2})+\gamma$ where
$\gamma={\sigma_{2}}{\|\boldsymbol{\pi}_{1}\|}+{\sigma_{1}}{\|\boldsymbol{\pi}_{2}\|}-1$.

Set $\mathbf{k}_{i}^{T}=\mathbf{h}_{i}^{T}(I-P_{ii})^{-1}$, $i=1,2$. Then

$\gamma=\|\boldsymbol{\pi}_{1}\|\mathbf{k}_{2}^{T}\operatorname{\mathbf{1}}+\|\boldsymbol{\pi}_{2}\|\mathbf{k}_{1}^{T}\operatorname{\mathbf{1}}.$

We find from (13) that $\mathbf{k}_{1}$ and $\mathbf{k}_{2}$ solve the linear
system

$\begin{bmatrix}\mathbf{k}_{1}^{T}&-\mathbf{k}_{2}^{T}\end{bmatrix}\begin{bmatrix}I-P_{11}&-P_{12}\\\
-P_{21}&I-P_{22}\end{bmatrix}=\begin{bmatrix}\widehat{}\boldsymbol{\pi}_{1}^{T}&-\widehat{}\boldsymbol{\pi}_{2}^{T}\end{bmatrix}.$

Since $(I-P)\operatorname{\mathbf{1}}=0$ and $\boldsymbol{\pi}^{T}(I-P)=0$,
this system is consistent and a solution is

$\begin{bmatrix}\mathbf{k}_{1}^{T}&-\mathbf{k}_{2}^{T}\end{bmatrix}=\begin{bmatrix}\widehat{}\boldsymbol{\pi}_{1}^{T}&-\widehat{}\boldsymbol{\pi}_{2}^{T}\end{bmatrix}(I-P+\operatorname{\mathbf{1}}\boldsymbol{\pi}^{T})^{-1}.$

Therefore,

$\gamma=\begin{bmatrix}\widehat{}\boldsymbol{\pi}_{1}^{T}&-\widehat{}\boldsymbol{\pi}_{2}^{T}\end{bmatrix}(I-P+\operatorname{\mathbf{1}}\boldsymbol{\pi}^{T})^{-1}\begin{bmatrix}\|\boldsymbol{\pi}_{2}\|\operatorname{\mathbf{1}}\\\
-\|\boldsymbol{\pi}_{1}\|\operatorname{\mathbf{1}}\end{bmatrix}.\qed$

Observe that

$\begin{split}(I-P+\operatorname{\mathbf{1}}\boldsymbol{\pi}^{T})^{-1}\begin{bmatrix}\|\boldsymbol{\pi}_{2}\|\operatorname{\mathbf{1}}\\\
-\|\boldsymbol{\pi}_{1}\|\operatorname{\mathbf{1}}\end{bmatrix}&=\left(I+\sum_{j=1}^{\infty}(P^{j}-\operatorname{\mathbf{1}}\boldsymbol{\pi}^{T})\right)\begin{bmatrix}\|\boldsymbol{\pi}_{2}\|\operatorname{\mathbf{1}}\\\
-\|\boldsymbol{\pi}_{1}\|\operatorname{\mathbf{1}}\end{bmatrix}\\\
&=\sum_{j=1}^{\infty}\left(P^{j}\begin{bmatrix}\|\boldsymbol{\pi}_{2}\|\operatorname{\mathbf{1}}\\\
-\|\boldsymbol{\pi}_{1}\|\operatorname{\mathbf{1}}\end{bmatrix}\right)\end{split}$

therefore

$\gamma=\sum_{j=1}^{\infty}\left(\begin{bmatrix}\widehat{}\boldsymbol{\pi}_{1}^{T}&-\widehat{}\boldsymbol{\pi}_{2}^{T}\end{bmatrix}P^{j}\begin{bmatrix}\|\boldsymbol{\pi}_{2}\|\operatorname{\mathbf{1}}\\\
-\|\boldsymbol{\pi}_{1}\|\operatorname{\mathbf{1}}\end{bmatrix}\right).$

By partitioning $P^{j}$ according to the partitioning of $P$, we have

$P^{j}=\begin{bmatrix}P_{11}^{(j)}&P_{12}^{(j)}\\\
P_{21}^{(j)}&P_{22}^{(j)}\end{bmatrix}.$

Since $P^{j}\operatorname{\mathbf{1}}=\operatorname{\mathbf{1}}$, we find that

$\begin{bmatrix}\widehat{}\boldsymbol{\pi}_{1}^{T}&-\widehat{}\boldsymbol{\pi}_{2}^{T}\end{bmatrix}P^{j}\begin{bmatrix}\|\boldsymbol{\pi}_{2}\|\operatorname{\mathbf{1}}\\\
-\|\boldsymbol{\pi}_{1}\|\operatorname{\mathbf{1}}\end{bmatrix}=\widehat{}\boldsymbol{\pi}_{1}^{T}P_{11}^{(j)}\operatorname{\mathbf{1}}+\widehat{}\boldsymbol{\pi}_{2}^{T}P_{22}^{(j)}\operatorname{\mathbf{1}}-1=-\lambda_{2}^{(j)},$

where $\lambda_{2}^{(j)}$ is the eigenvalue different from 1 of the $2\times
2$ aggregated matrix

$\begin{bmatrix}\widehat{}\boldsymbol{\pi}_{1}^{T}P_{11}^{(j)}\operatorname{\mathbf{1}}&\widehat{}\boldsymbol{\pi}_{1}^{T}P_{12}^{(j)}\operatorname{\mathbf{1}}\\\
\widehat{}\boldsymbol{\pi}_{2}^{T}P_{21}^{(j)}\operatorname{\mathbf{1}}&\widehat{}\boldsymbol{\pi}_{2}^{T}P_{22}^{(j)}\operatorname{\mathbf{1}}\end{bmatrix}.$

By means of formal manipulations, we can derive different expressions for the
constant $\gamma$. In fact, we have

$\gamma=\|\boldsymbol{\pi}_{1}\|\theta-\|\boldsymbol{\pi}_{2}\|,$ (16)

where $\theta>0$ has one of the following equivalent expressions

$\displaystyle\theta$
$\displaystyle=\hat{\boldsymbol{\pi}}_{2}^{T}\left(I+(I-P_{22})^{-1}P_{21}(I-P_{1}+\operatorname{\mathbf{1}}\hat{\boldsymbol{\pi}}_{1}^{T})^{-1}P_{12}\right)(I-P_{22})^{-1}\operatorname{\mathbf{1}}$
(17)
$\displaystyle=\begin{bmatrix}\mathbf{0}^{T}&\hat{\boldsymbol{\pi}}_{2}^{T}\end{bmatrix}\left(I-P+\begin{bmatrix}\operatorname{\mathbf{1}}\\\
\mathbf{0}\end{bmatrix}\begin{bmatrix}\hat{\boldsymbol{\pi}}_{1}^{T}&\mathbf{0}^{T}\end{bmatrix}\right)^{-1}\begin{bmatrix}\mathbf{0}\\\
\operatorname{\mathbf{1}}\end{bmatrix}$ (18)
$\displaystyle=\hat{\boldsymbol{\pi}}_{2}^{T}\left(I-P_{22}-P_{21}(I-P_{11}+\operatorname{\mathbf{1}}\hat{\boldsymbol{\pi}}_{1}^{T})^{-1}P_{12}\right)^{-1}\operatorname{\mathbf{1}}$
(19)
$\displaystyle=\hat{\boldsymbol{\pi}}_{2}^{T}\left(I-P_{2}+\frac{1}{1+\hat{\boldsymbol{\pi}}_{1}^{T}S\operatorname{\mathbf{1}}}P_{21}S\operatorname{\mathbf{1}}\hat{\boldsymbol{\pi}}_{1}^{T}SP_{12}\right)^{-1}\operatorname{\mathbf{1}},\quad
S=(I-P_{11})^{-1}.$ (20)

##  4 Application of Theorem 3.3

Here we apply the main result of the previous section to structured stochastic
matrices.

###  4.1 Periodic Markov chains

Let $P$ be the transition matrix of an ergodic Markov chain with $n$ states.
The period $p_{i}$ of state $i$ is the greatest common divisor of all natural
numbers $m$ such that $(P^{m})_{i,i}>0$. A Markov chain is called periodic if
$p_{i}\geq 2$ for all $1\leq i\leq n$, and it is called aperiodic otherwise.
The period of a periodic Markov chain is the greatest common divisor of
periods of all states. It is well-known that if the Markov chain is periodic,
then $P$ is permutationally similar to a block-cyclic matrix (see [8]).

Given a periodic Markov chain with $n$ states and period $d$, we may assume
that the transition matrix $P$ is given by

$P=\begin{bmatrix}P_{11}&P_{12}\\\
P_{21}&P_{22}\end{bmatrix}=\left[\begin{array}[]{c|cccc}0&0&\ldots&0&A_{d}\\\
\hline\cr A_{1}&0&\ddots&\ddots&0\\\ 0&A_{2}&\ddots&\ddots&\vdots\\\
\vdots&\ddots&\ddots&0&\vdots\\\ 0&\ldots&0&A_{d-1}&0\end{array}\right]$ (21)

where all block diagonal matrices are square and each of $A_{i}$’s is a
rectangular stochastic matrix. Let $n_{i}$ be the size of $i^{\text{th}}$
block diagonal matrix of $P$ for $i=1,\dots,d$. A set of the states
corresponding to a block diagonal matrix is called a cyclic class. The
cyclicity index is the number of cyclic classes. One can find that each of

$A_{d}A_{d-1}\cdots A_{1},\quad A_{1}A_{d}\cdots A_{2},\quad\dots\quad,\quad
A_{d-1}\cdots A_{1}A_{d}$ (22)

is a square stochastic matrix and its corresponding Markov chain is aperiodic.
It can be seen that

$(I-P_{22})^{-1}=\left[\begin{array}[]{c|c|c|c|c}I&0&&0&0\\\
A_{2}&I&&\vdots&\vdots\\\ A_{3}A_{2}&A_{3}&\cdots&0&0\\\ \vdots&\vdots&&I&0\\\
A_{d-1}\cdots A_{3}A_{2}&A_{d-1}\cdots
A_{4}A_{3}&&A_{d-1}&I\end{array}\right].$ (23)

It is immediate to find that the stochastic complements $P_{1}$ and $P_{2}$ of
$P_{11}$ and $P_{22}$ are, respectively,

$P_{1}=A_{d}A_{d-1}\cdots A_{1},\quad
P_{2}=\begin{bmatrix}0&\ldots&\ldots&0&A_{1}A_{d}\\\
A_{2}&0&\ddots&\ddots&0\\\ 0&A_{3}&\ddots&\ddots&\vdots\\\
\vdots&\ddots&\ddots&0&\vdots\\\ 0&\ldots&0&A_{d-1}&0\end{bmatrix}.$ (24)

We note that the censored Markov chain for $P_{2}$ is periodic.

######  Proposition 4.1.

Let $P$ be the transition matrix of a periodic Markov chain. Suppose that $P$
is of form (21). Then

$\kappa(P)=\kappa(P_{1})+\kappa(P_{2})+\frac{1}{2}$

where $P_{1}$ and $P_{2}$ are given in (24).

###### Proof.

It suffices to show that $\gamma$ in (15) is $\frac{1}{2}$. In order to find
$\gamma$, we shall find the value of $\theta$ introduced in Equation (16).
From (17), we may look at $\theta$ as the sum of two terms $\theta_{1}$ and
$\theta_{2}$, where

$\begin{split}&\theta_{1}=\hat{}\boldsymbol{\pi}_{2}^{T}(I-P_{22})^{-1}\operatorname{\mathbf{1}},\\\
&\theta_{2}=\hat{}\boldsymbol{\pi}_{2}^{T}(I-P_{22})^{-1}P_{21}(I-P_{1}+\operatorname{\mathbf{1}}\hat{}\boldsymbol{\pi}_{1}^{T})^{-1}P_{12}(I-P_{22})^{-1}\operatorname{\mathbf{1}}.\end{split}$

We first find expressions of $\hat{}\boldsymbol{\pi}_{2}^{T}$ and
$(I-P_{22})^{-1}\operatorname{\mathbf{1}}$. Let

$\boldsymbol{\pi}^{T}=[\operatorname{\mathbf{s}}_{1}^{T},\operatorname{\mathbf{s}}_{2}^{T},\ldots,\operatorname{\mathbf{s}}_{d}^{T}]$

be conformally partitioned with $P$ and be the stationary distribution vector.
From the condition $\boldsymbol{\pi}^{T}P=\boldsymbol{\pi}^{T}$, we have
$\operatorname{\mathbf{s}}_{i}^{T}=\operatorname{\mathbf{s}}_{i+1}^{T}A_{i}$,
$i=1,\ldots,d-1$, and
$\operatorname{\mathbf{s}}_{d}^{T}=\operatorname{\mathbf{s}}_{1}^{T}A_{d}$.
Then

$\operatorname{\mathbf{s}}_{2}^{T}=\operatorname{\mathbf{s}}_{1}^{T}A_{d}A_{d-1}\cdots
A_{2},\quad\operatorname{\mathbf{s}}_{3}^{T}=\operatorname{\mathbf{s}}_{1}^{T}A_{d}A_{d-1}\cdots
A_{3},\quad\dots\quad,\quad\operatorname{\mathbf{s}}_{d-1}^{T}=\operatorname{\mathbf{s}}_{1}^{T}A_{d}A_{d-1}.$

Moreover, since $A_{i}\operatorname{\mathbf{1}}=\operatorname{\mathbf{1}}$, it
follows that
$\operatorname{\mathbf{s}}_{i}^{T}\operatorname{\mathbf{1}}=\operatorname{\mathbf{s}}_{i+1}^{T}\operatorname{\mathbf{1}}$
and
$\operatorname{\mathbf{s}}_{d}^{T}\operatorname{\mathbf{1}}=\operatorname{\mathbf{s}}_{1}^{T}\operatorname{\mathbf{1}}$.
Hence, $\|\operatorname{\mathbf{s}}_{j}\|=\frac{1}{d}$ for $1\leq j\leq d$.
Note that $\boldsymbol{\pi}_{1}=\operatorname{\mathbf{s}}_{1}$. Therefore,

$\boldsymbol{\pi}^{T}=\frac{1}{d}\begin{bmatrix}\hat{}\boldsymbol{\pi}_{1}^{T}&\hat{}\boldsymbol{\pi}_{1}^{T}(A_{d}A_{d-1}\cdots
A_{2})&\hat{}\boldsymbol{\pi}_{1}^{T}(A_{d}A_{d-1}\cdots
A_{3})&\cdots&\hat{}\boldsymbol{\pi}_{1}^{T}A_{d}\end{bmatrix},$

whence

$\hat{}\boldsymbol{\pi}_{2}^{T}=\frac{1}{d-1}\begin{bmatrix}\hat{}\boldsymbol{\pi}_{1}^{T}(A_{d}A_{d-1}\cdots
A_{2})&\hat{}\boldsymbol{\pi}_{1}^{T}(A_{d}A_{d-1}\cdots
A_{3})&\cdots&\hat{}\boldsymbol{\pi}_{1}^{T}A_{d}\end{bmatrix}.$ (25)

Observing $(I-P_{22})^{-1}$ in (23), we have

$\displaystyle(I-P_{22})^{-1}\operatorname{\mathbf{1}}=\begin{bmatrix}\operatorname{\mathbf{1}}^{T}&2\operatorname{\mathbf{1}}^{T}&\cdots&(d-1)\operatorname{\mathbf{1}}^{T}\end{bmatrix}^{T},$
(26) $\displaystyle
P_{12}(I-P_{22})^{-1}\operatorname{\mathbf{1}}=(d-1)\operatorname{\mathbf{1}},$
$\displaystyle(I-P_{22})^{-1}P_{21}\operatorname{\mathbf{1}}=\operatorname{\mathbf{1}}.$

By using the expression of $\hat{}\boldsymbol{\pi}_{2}$ in (25) and the
expression of $(I-P_{22})^{-1}\operatorname{\mathbf{1}}$ in (26) together with
the property $A_{i}\operatorname{\mathbf{1}}=\operatorname{\mathbf{1}}$, we
find that

$\theta_{1}=\frac{1}{d-1}\sum_{i=1}^{d-1}i\hat{}\boldsymbol{\pi}_{1}^{T}\operatorname{\mathbf{1}}=\frac{1}{2}d.$

Note that since
$(I-P_{1}+\operatorname{\mathbf{1}}\hat{}\boldsymbol{\pi}_{1}^{T})\operatorname{\mathbf{1}}=\operatorname{\mathbf{1}}$,
we have
$(I-P_{1}+\operatorname{\mathbf{1}}\hat{}\boldsymbol{\pi}_{1}^{T})^{-1}\operatorname{\mathbf{1}}=\operatorname{\mathbf{1}}$.
Concerning $\theta_{2}$, we find from (26) that

$\displaystyle\theta_{2}$
$\displaystyle=\hat{}\boldsymbol{\pi}_{2}^{T}(I-P_{22})^{-1}P_{21}(I-P_{1}+\operatorname{\mathbf{1}}\hat{}\boldsymbol{\pi}_{1}^{T})^{-1}P_{12}(I-P_{22})^{-1}\operatorname{\mathbf{1}}=d-1.$

Thus, $\theta=\theta_{1}+\theta_{2}=\frac{3}{2}d-1$ and therefore
$\gamma=\|\boldsymbol{\pi}_{1}\|\theta-\|\boldsymbol{\pi}_{2}\|=\frac{1}{2}$.
∎

######  Theorem 4.2.

Let $P$ be the transition matrix of a periodic Markov chain. Suppose that $P$
is of form (21). Then

$\kappa(P)=d\kappa(A_{d}A_{d-1}\cdots A_{1})+n-dn_{1}+\frac{d-1}{2}.$ (27)

Moreover, if $n_{1}=\cdots=n_{d}$ then

$\kappa(P)=d\kappa(A_{d}A_{d-1}\cdots A_{1})+\frac{d-1}{2}.$ (28)

###### Proof.

Consider $P_{2}$ in (24) with $d\geq 3$. Note that the Markov chain
corresponding to $P_{2}$ is periodic. Applying Proposition 4.1, we find that

$\kappa(P_{2})=\kappa(A_{1}A_{d}A_{d-1}\cdots
A_{2})+\kappa(P_{3})+\frac{1}{2},$

where

$P_{3}=\left[\begin{array}[]{c|ccc}0&\ldots&0&A_{2}A_{1}A_{d}\\\ \hline\cr
A_{3}&\ddots&\ddots&0\\\ &\ddots&\ddots&\vdots\\\
0&&A_{d-1}&0\end{array}\right].$

Hence, $\kappa(P)=\kappa(A_{d}A_{d-1}\cdots
A_{1})+\kappa(A_{1}A_{d}A_{d-1}\cdots A_{2})+\kappa(P_{3})+1$. If $d\geq 4$,
one can apply the proposition to $P_{3}$ with the partition. In this manner,
recursively applying the proposition, we obtain

$\displaystyle\kappa(P)$ $\displaystyle=\kappa(A_{d}A_{d-1}\cdots
A_{1})+\kappa(A_{1}A_{d}A_{d-1}\cdots A_{2})+\cdots$
$\displaystyle+\kappa(A_{d-1}A_{d-2}\cdots A_{1}A_{d})+\frac{d-1}{2}.$

Recall that $A_{i}$ is of size $n_{i+1}\times n_{i}$ for $1\leq i\leq d-1$ and
$A_{d}$ is of size $n_{1}\times n_{d}$. Note that $A_{i-1}A_{k-2}\cdots
A_{1}A_{d}A_{d-1}\cdots A_{i}$ is of size $n_{i}\times n_{i}$ for $1\leq i\leq
d$. From (9), we obtain

$\displaystyle\kappa(A_{k-1}A_{k-2}\cdots A_{1}A_{d}A_{d-1}\cdots A_{k})$
$\displaystyle=\kappa(A_{k-2}A_{k-3}\cdots A_{1}A_{d}A_{d-1}\cdots A_{k-1})$
$\displaystyle+n_{k}-n_{k-1}$

for $2\leq k\leq d$. It follows that

$\displaystyle\kappa(A_{d}$ $\displaystyle A_{d-1}\cdots
A_{1})+\cdots+\kappa(A_{d-2}\cdots A_{1}A_{d}A_{d-1})+\kappa(A_{d-1}\cdots
A_{1}A_{d})$ $\displaystyle=\;d\kappa(A_{d}A_{d-1}\cdots
A_{1})+(d-1)(n_{2}-n_{1})+\cdots$
$\displaystyle+\;2(n_{d-1}-n_{d-2})+n_{d}-n_{d-1}=\;d\kappa(A_{d}A_{d-1}\cdots
A_{1})+n-dn_{1}.\qed$

######  Remark 4.3.

The characteristic polynomials of the matrices $P$ in (21) and
$P_{1}=A_{d}A_{d-1}\cdots A_{1}$ are related by the equation $\det(\lambda
I-P)=\lambda^{\ell}\det(\lambda^{d}I-P_{1})$, where $\ell=n-dn_{1}$.
Therefore, from the computational point of view, if the eigenvalues of $P_{1}$
are explicitly known, then also the eigenvalues of $P$ are explicitly known
and we can recover $\kappa(P)$ from (7). On the other hand, if the eigenvalues
of $P_{1}$ are not known, but the value of $\kappa(P_{1})$ is available, then
Kemeny’s constant of $P$ can be directly obtained from (27).

######  Example 4.4.

A collection of random walks on undirected graphs is one of the most
accessible families of Markov chains. If a random walk is periodic, then the
underlying graph is necessarily bipartite. Let $P$ be an irreducible
stochastic matrix with the structure

$P=\begin{bmatrix}0&P_{12}\\\ P_{21}&0\end{bmatrix}$

where $P_{12}$ and $P_{21}$ have size $n_{1}\times n_{2}$ and $n_{2}\times
n_{1}$, respectively. Then

$\kappa(P)=2\kappa(P_{1})+n_{1}-n_{2}+\frac{1}{2},$ (29)

where $P_{1}=P_{12}P_{21}$.

######  Example 4.5.

Consider the transition matrix $P$ in (21). Suppose that for $i=1,\dots,d$,
$A_{i}$ is the matrix with all entries equal to $\frac{1}{n_{i}}$. Then
$A_{d}A_{d-1}\cdots A_{1}$ is the $n_{1}\times n_{1}$ matrix with all entries
equal to $\frac{1}{n_{1}}$. Since $\kappa(A_{d}A_{d-1}\cdots A_{1})=n_{1}-1$,
we have

$\kappa(P)=n-\frac{d+1}{2}.$

Now we provide a lower bound on Kemeny’s constant of a periodic Markov chain.

######  Corollary 4.6.

Let $P$ be the transition matrix of a periodic Markov chain. Suppose that
$P=[p_{xy}]$ is of form (21) with $n_{1}\leq n_{i}$ for $1\leq i\leq d$. Then

$\kappa(P)\geq n-\frac{dn_{1}+1}{2},$

with equality when $p_{xy}$ is given as in (30).

###### Proof.

From (27), $\kappa(P)=d\kappa(A_{d}A_{d-1}\cdots
A_{1})+n-dn_{1}+\frac{d-1}{2}$. It is known in [26, Remark 2.14] that given an
$m\times m$ irreducible stochastic matrix $P$, $\kappa(P)\geq\frac{m-1}{2}$ if
and only if $P$ is the adjacency matrix of a directed $m$-cycle. Hence, it is
enough to find a periodic Markov chain such that $A_{d}A_{d-1}\cdots A_{1}$ is
the adjacency matrix of a directed $n_{1}$-cycle.

Note that $n_{1}\leq n_{i}$ for $1\leq i\leq d$. We may suppose that for
$j=1,\dots,d$, the cyclic class corresponding to $j^{\text{th}}$ block
diagonal matrix of $P=[p_{xy}]$ is partitioned into $n_{1}$ subsets
$C_{1}^{j},\dots,C_{n_{1}}^{j}$. Let $p_{xy}$ be given as follows:

$\displaystyle\begin{split}p_{xy}=\begin{cases}\frac{1}{|C_{\ell}^{d}|},&\text{if
$x\in C_{\ell}^{1}$ and $y\in C_{\ell}^{d}$ for $1\leq\ell\leq n_{1}$;}\\\
\frac{1}{|C_{\ell-1}^{j-1}|},&\text{if $x\in C_{\ell}^{j}$ and $y\in
C_{\ell}^{j-1}$ for $1\leq\ell\leq n_{1}$ and $3\leq j\leq d$;}\\\ 1,&\text{if
$x\in C_{\ell}^{2}$ and $y\in C_{\ell-1}^{1}$ for $2\leq\ell\leq n_{1}$, or
$x\in C_{1}^{2}$ and $y\in C_{n_{1}}^{1}$;}\\\ 0,&\text{otherwise.}\\\
\end{cases}\end{split}$ (30)

It can be seen that $A_{d}A_{d-1}\cdots A_{1}$ is the adjacency matrix of a
directed $n_{1}$-cycle, which completes the proof. ∎

###  4.2 Kronecker product of stochastic matrices

Given stochastic matrices $A$ and $B$, $A\otimes B$ is also stochastic where
$\otimes$ denotes the Kronecker product. We will provide an expression of
Kemeny’s constant of $A\otimes B$ after applying Theorem 3.3.

Partition $P=A\otimes B$ as follows:

$P=\begin{bmatrix}P_{11}&P_{12}\\\
P_{21}&P_{22}\end{bmatrix}=\left[\begin{array}[]{c|ccc}a_{11}B&a_{12}B&\dots&a_{1n}B\\\
\hline\cr a_{21}B&&&\\\ \vdots&&A_{1}\otimes B&\\\
a_{n1}B&&&\end{array}\right].$ (31)

We use $\mathbf{e}_{i}$ to denote the vector with a $1$ in the $i\text{th}$
entry and zeros elsewhere. We have the following

######  Proposition 4.7.

Let $A=[a_{ij}]$ and $B$ be stochastic matrices with stationary distribution
vectors $\mathbf{x}=(x_{i})$ and $\mathbf{y}$, respectively. Let $P=A\otimes
B$. Then

$\kappa(P)=\kappa(P_{1})+\kappa(P_{2})+\frac{1}{1-x_{1}}(\mathbf{e}_{1}^{T}(I-A+\operatorname{\mathbf{1}}\mathbf{x}^{T})^{-1}\mathbf{e}_{1}-x_{1}),$

where $P_{1}$ and $P_{2}$ are the stochastic complements in (31) of $P_{11}$
and $P_{22}$, respectively.

###### Proof.

Let $A=[a_{ij}]$ and $B$ be matrices of size $n\times n$ and $m\times m$,
respectively. It suffices to provide an explicit expression for $\gamma$ in
Theorem 3.3.

Let $\mathbf{x}=(x_{i})\in\mathbb{R}^{n}$ and $\mathbf{y}\in\mathbb{R}^{m}$ be
the stationary distribution vector for $A$ and $B$, respectively. Then
$\mathbf{x}^{T}A=\mathbf{x}^{T}$, $\mathbf{y}^{T}B=\mathbf{y}^{T}$ and
$\mathbf{x}^{T}\operatorname{\mathbf{1}}=\mathbf{y}^{T}\operatorname{\mathbf{1}}=1$.
Let
$\boldsymbol{\pi}^{T}=\begin{bmatrix}\boldsymbol{\pi}_{1}^{T}&\boldsymbol{\pi}_{2}^{T}\end{bmatrix}$
be conformally partitioned with $P$ and be the stationary distribution vector
for $P$. Since $\boldsymbol{\pi}=\mathbf{x}\otimes\mathbf{y}$, we have
$\boldsymbol{\pi}_{1}=x_{1}\mathbf{y}$ and
$\boldsymbol{\pi}_{2}=(x_{2},\dots,x_{n})\otimes\mathbf{y}$. So,
$\|\boldsymbol{\pi}_{1}\|=x_{1}$ and $\|\boldsymbol{\pi}_{2}\|=1-x_{1}$.

Set $q=mn-m$. Consider

$\begin{bmatrix}\|\boldsymbol{\pi}_{2}\|\operatorname{\mathbf{1}}_{m}\\\
-\|\boldsymbol{\pi}_{1}\|\operatorname{\mathbf{1}}_{q}\end{bmatrix}=-\|\boldsymbol{\pi}_{1}\|\begin{bmatrix}\operatorname{\mathbf{1}}_{m}\\\
\operatorname{\mathbf{1}}_{q}\end{bmatrix}+\begin{bmatrix}\operatorname{\mathbf{1}}_{m}\\\
\mathbf{0}_{q}\end{bmatrix}.$

Then

$(I-P+\operatorname{\mathbf{1}}\boldsymbol{\pi}^{T})^{-1}\begin{bmatrix}\|\boldsymbol{\pi}_{2}\|\operatorname{\mathbf{1}}_{m}\\\
-\|\boldsymbol{\pi}_{1}\|\operatorname{\mathbf{1}}_{q}\end{bmatrix}=-\|\boldsymbol{\pi}_{1}\|\begin{bmatrix}\operatorname{\mathbf{1}}_{m}\\\
\operatorname{\mathbf{1}}_{q}\end{bmatrix}+(I-P+\operatorname{\mathbf{1}}\boldsymbol{\pi}^{T})^{-1}(\mathbf{e}_{1}\otimes\operatorname{\mathbf{1}}_{m}).$
(32)

We claim that

$(I-P+\operatorname{\mathbf{1}}\boldsymbol{\pi}^{T})^{-1}(\mathbf{e}_{1}\otimes\operatorname{\mathbf{1}}_{m})=\mathbf{z}\otimes\operatorname{\mathbf{1}}_{m}$

where
$\mathbf{z}=(I-A+\operatorname{\mathbf{1}}_{n}\mathbf{x}^{T})^{-1}\mathbf{e}_{1}$.
Consider

$(I-A\otimes
B+\operatorname{\mathbf{1}}(\mathbf{x}\otimes\mathbf{y})^{T})(\mathbf{z}\otimes\operatorname{\mathbf{1}}_{m})=\mathbf{e}_{1}\otimes\operatorname{\mathbf{1}}_{m}.$

This system has the form

$(\mathbf{z}\otimes\operatorname{\mathbf{1}}_{m})-(A\otimes
B)(\mathbf{z}\otimes\operatorname{\mathbf{1}}_{m})+(\operatorname{\mathbf{1}}_{n}\otimes\operatorname{\mathbf{1}}_{m})(\mathbf{x}\otimes\mathbf{y})^{T}(\mathbf{z}\otimes\operatorname{\mathbf{1}}_{m})=(\mathbf{e}_{1}\otimes\operatorname{\mathbf{1}}_{m}).$

This implies that

$(\mathbf{z}-A\mathbf{z}+(\operatorname{\mathbf{1}}_{n}\mathbf{x}^{T})\mathbf{z}-\mathbf{e}_{1})\otimes\operatorname{\mathbf{1}}_{m}=\mathbf{0}.$

Hence, our desired claim is established.

Now, applying (32) together with the claim to the expression of $\gamma$ given
in Theorem 3.3, we obtain

$\begin{split}\gamma=&\begin{bmatrix}\widehat{}\boldsymbol{\pi}_{1}^{T}&-\widehat{}\boldsymbol{\pi}_{2}^{T}\end{bmatrix}(I-P+\operatorname{\mathbf{1}}\boldsymbol{\pi}^{T})^{-1}\begin{bmatrix}\|\boldsymbol{\pi}_{2}\|\operatorname{\mathbf{1}}_{m}\\\
-\|\boldsymbol{\pi}_{1}\|\operatorname{\mathbf{1}}_{q}\end{bmatrix}\\\
=&\begin{bmatrix}\hat{}\boldsymbol{\pi}_{1}^{T}&-\hat{}\boldsymbol{\pi}_{2}^{T}\end{bmatrix}(-\|\boldsymbol{\pi}_{1}\|\operatorname{\mathbf{1}}+\mathbf{z}\otimes\operatorname{\mathbf{1}}_{m})\\\
=&(\begin{bmatrix}1&-\frac{1}{1-\|\boldsymbol{\pi}_{1}\|}(x_{2},\ldots,x_{n})\end{bmatrix}\otimes\mathbf{y}^{T})(\mathbf{z}\otimes\operatorname{\mathbf{1}}_{m})\\\
=&\frac{1}{1-x_{1}}\begin{bmatrix}1-x_{1}&x_{2}&\dots&x_{n}\end{bmatrix}\mathbf{z}\\\
=&\frac{1}{1-x_{1}}(\mathbf{e}_{1}-\mathbf{x})^{T}(I-A+\operatorname{\mathbf{1}}_{n}\mathbf{x}^{T})^{-1}\mathbf{e}_{1}\\\
=&\frac{1}{1-x_{1}}(\mathbf{e}_{1}^{T}(I-A+\operatorname{\mathbf{1}}_{n}\mathbf{x}^{T})^{-1}\mathbf{e}_{1}-x_{1}).\end{split}$

Therefore, we obtain the desired result. ∎

###  4.3 Sub-stochastic matrices with constant row sums

Here is the result of this subsection.

######  Proposition 4.8.

Let

$P=\begin{bmatrix}P_{11}&P_{12}\\\ P_{21}&P_{22}\end{bmatrix},$

be an irreducible stochastic matrix, where
$P_{11}\operatorname{\mathbf{1}}=r_{1}\operatorname{\mathbf{1}}$ and
$P_{22}\operatorname{\mathbf{1}}=r_{2}\operatorname{\mathbf{1}}$ for some
$0\leq r_{1},r_{2}<1$. Denote by $P_{1}$ and $P_{2}$ the stochastic
complements of $P_{11}$ and $P_{22}$, respectively. Then

$\displaystyle\|\boldsymbol{\pi}_{1}\|=\frac{1-r_{2}}{2-r_{1}-r_{2}},\quad\|\boldsymbol{\pi}_{2}\|=\frac{1-r_{1}}{2-r_{1}-r_{2}},$
$\displaystyle\kappa(P)=\kappa(P_{1})+\kappa(P_{2})+\frac{1}{2-r_{1}-r_{2}}.$

###### Proof.

Note that $P_{12}\operatorname{\mathbf{1}}=(1-r_{1})\operatorname{\mathbf{1}}$
and $P_{21}\operatorname{\mathbf{1}}=(1-r_{2})\operatorname{\mathbf{1}}$. From
(17),

$\displaystyle\theta$
$\displaystyle=\hat{\boldsymbol{\pi}}_{2}^{T}(I-P_{22})^{-1}\operatorname{\mathbf{1}}+\hat{\boldsymbol{\pi}}_{2}^{T}(I-P_{22})^{-1}P_{21}(I-P_{1}+\operatorname{\mathbf{1}}\hat{\boldsymbol{\pi}}_{1}^{T})^{-1}P_{12}(I-P_{22})^{-1}\operatorname{\mathbf{1}}$
$\displaystyle=\frac{1}{1-r_{2}}+\frac{1-r_{1}}{1-r_{2}}=\frac{2-r_{1}}{1-r_{2}}.$

On the other hand, we find from (4) that the aggregated matrix $S$ is given by

$S=\begin{bmatrix}r_{1}&1-r_{1}\\\ 1-r_{2}&r_{2}\end{bmatrix},$

and so $\alpha_{1}=\|\boldsymbol{\pi}_{1}\|=(1-r_{2})/(2-r_{1}-r_{2})$,
$\alpha_{2}=\|\boldsymbol{\pi}_{2}\|=(1-r_{1})/(2-r_{1}-r_{2})$, and
$\|\boldsymbol{\pi}_{1}\|(1+\theta)=(3-r_{1}-r_{2})/(2-r_{1}-r_{2})$. Whence,
in view of (15) and (16), we deduce that

$\kappa(P)=\kappa(P_{1})+\kappa(P_{2})+\frac{1}{2-r_{1}-r_{2}}.\qed$

##  5 Some bounds

Assume we are given the values of $\kappa(P_{1})$ and $\kappa(P_{2})$.
According to Equation (15), we can provide lower and upper bounds to the value
of $\kappa(P)$ once we are given bounds to the constant $\gamma$. In order to
do this, we need to determine upper and lower bounds to the value of
$\|\boldsymbol{\pi}_{1}\|$.

From the equation $\boldsymbol{\pi}^{T}P=\boldsymbol{\pi}^{T}$ we find that

$\begin{split}\boldsymbol{\pi}_{1}^{T}=\boldsymbol{\pi}_{2}^{T}P_{21}(I-P_{11})^{-1},\quad\quad\boldsymbol{\pi}_{2}^{T}=\boldsymbol{\pi}_{1}^{T}P_{12}(I-P_{22})^{-1}.\end{split}$

Taking the infinity norms of both sides and using the identity
$\|\mathbf{x}^{T}\|_{\infty}=\|\mathbf{x}\|_{1}$, we obtain

$\|\boldsymbol{\pi}_{1}\|\leq\|\boldsymbol{\pi}_{2}\|\,\|P_{21}\|_{\infty}\|(I-P_{11})^{-1}\|_{\infty}\leq\|\boldsymbol{\pi}_{2}\|\frac{\|P_{21}\|_{\infty}}{1-\|P_{11}\|_{\infty}},$

where the latter inequality is valid if $\|P_{11}\|_{\infty}<1$. A similar
inequality can be obtained for $\|\boldsymbol{\pi}_{2}\|$. Combining both
inequalities, under the assumption
$\|P_{11}\|_{\infty},\|P_{22}\|_{\infty}<1$, we get

$\frac{1-\|P_{22}\|_{\infty}}{\|P_{12}\|_{\infty}}\leq\frac{\|\boldsymbol{\pi}_{1}\|}{\|\boldsymbol{\pi}_{2}\|}\leq\frac{\|P_{21}\|_{\infty}}{1-\|P_{11}\|_{\infty}}.$

Moreover, since $\|\boldsymbol{\pi}_{2}\|=1-\|\boldsymbol{\pi}_{1}\|$, we
obtain

$\frac{1-\|P_{22}\|_{\infty}}{1-\|P_{22}\|_{\infty}+\|P_{12}\|_{\infty}}\leq\|\boldsymbol{\pi}_{1}\|\leq\frac{\|P_{21}\|_{\infty}}{1-\|P_{11}\|_{\infty}+\|P_{21}\|_{\infty}}.$
(33)

Since the matrix $P$ is known, the above bounds can be actually computed at a
low computational cost.

Now a bound to the constant $\gamma$ can be obtained by relying on Equation
(16) coupled with (17) (one can use (18) or (19)). From (16),
$\gamma=(1+\theta)\|\boldsymbol{\pi}_{1}\|-1$. We see from (33) that

$(\theta+1)\frac{1-\|P_{22}\|_{\infty}}{1-\|P_{22}\|_{\infty}+\|P_{12}\|_{\infty}}-1\leq\gamma\leq(\theta+1)\frac{\|P_{21}\|_{\infty}}{1-\|P_{11}\|_{\infty}+\|P_{21}\|_{\infty}}-1.$

Note that $\|(I-P_{22})^{-1}P_{21}\|_{\infty}=1$. Concerning $\theta$, it
follows from (17) that

$\theta\leq\;\|(1-P_{22})^{-1}\|_{\infty}\left(1+\|P_{12}\|_{\infty}\|(I-P_{1}+\operatorname{\mathbf{1}}\widehat{}\boldsymbol{\pi}_{1}^{T})^{-1}\|_{\infty}\right).$
(34)

The upper bound to $\gamma$ is expressed in terms of $(I-P_{22})^{-1}$,
$(I-P_{1}+\operatorname{\mathbf{1}}\widehat{}\boldsymbol{\pi}_{1}^{T})^{-1}$
and the norms of the blocks $P_{ij}$.

For another bound on $\gamma$, consider the expression of $\gamma$ in Theorem
3.3, i.e.,

$\gamma=\begin{bmatrix}\widehat{}\boldsymbol{\pi}_{1}^{T}&-\widehat{}\boldsymbol{\pi}_{2}^{T}\end{bmatrix}(I-P+\operatorname{\mathbf{1}}\boldsymbol{\pi}^{T})^{-1}\begin{bmatrix}||\boldsymbol{\pi}_{2}\|\operatorname{\mathbf{1}}\\\
-\|\boldsymbol{\pi}_{1}\|\operatorname{\mathbf{1}}\end{bmatrix}.$

Taking the infinity norm of both sides in the above equation yields

$|\gamma|\leq
2\|(I-P+\operatorname{\mathbf{1}}\boldsymbol{\pi}^{T})^{-1}\|_{\infty}\max(\|\boldsymbol{\pi}_{1}\|,1-\|\boldsymbol{\pi}_{1}\|).$

###  5.1 A perturbation result

Let $P$ be an $n\times n$ stochastic and irreducible matrix, and $E$ be
$n\times n$ matrix such that $E\operatorname{\mathbf{1}}=0$. Let $\epsilon>0$.
Then $(P+\epsilon E)\operatorname{\mathbf{1}}=\operatorname{\mathbf{1}}$ .
Assume that $P(\epsilon):=P+\epsilon E$ is stochastic and $\|E\|_{\infty}\leq
1$.

Here, our goal is to relate $\kappa(P(\epsilon))$ and $\kappa(P)$. We have

$\begin{split}\kappa(P(\epsilon))=\hbox{Tr}((I-P(\epsilon)+\operatorname{\mathbf{1}}\mathbf{h}^{T})^{-1}),\quad\quad\kappa(P)=\hbox{Tr}((I-P+\operatorname{\mathbf{1}}\mathbf{h}^{T})^{-1}),\end{split}$

where $\mathbf{h}$ is any vector such that
$\mathbf{h}^{T}\operatorname{\mathbf{1}}=1$. By subtracting both sides of the
above equations we get

$\kappa(P(\epsilon))-\kappa(P)=\hbox{Tr}((I-P(\epsilon)+\operatorname{\mathbf{1}}\mathbf{h}^{T})^{-1}(P(\epsilon)-P)(I-P+\operatorname{\mathbf{1}}\mathbf{h}^{T})^{-1}).$

Neglecting $O(\epsilon^{2})$ terms, we obtain

$\begin{split}\kappa(P(\epsilon))-\kappa(P)&=\epsilon\,\hbox{Tr}((I-P+\operatorname{\mathbf{1}}\mathbf{h}^{T})^{-1}E(I-P+\operatorname{\mathbf{1}}\mathbf{h}^{T})^{-1})\\\
&=\epsilon\,\hbox{Tr}((I-P+\operatorname{\mathbf{1}}\mathbf{h}^{T})^{-2}E)+O(\epsilon^{2}).\end{split}$

This estimate can be also deduced as a specific case of [7, Lemma 3.2].

From the Cauchy–Schwarz inequality, we have
$\hbox{Tr}(AB)\leq\|A\|_{F}\|B\|_{F}$, where $\|\cdot\|_{F}$ is the Frobenius
norm. It follows that

$|\kappa(P(\epsilon))-\kappa(P)|\leq\epsilon\|(I-P+\operatorname{\mathbf{1}}\mathbf{h}^{T})^{-2}\|_{F}\|E\|_{F}+O(\epsilon^{2}).$
(35)

Combining the above bound with Example 4.4 yields the following result that
concerns Kemeny’s constant of a matrix associated with an almost bipartite
graph.

######  Corollary 5.1.

Let $P$ be the stochastic matrix defined in Example 4.4 and $E$ be a matrix
such that $E\operatorname{\mathbf{1}}=0$ and $\|E\|_{\infty}\leq 1$. Set
$P(\epsilon)=P+\epsilon E$ and assume that $P(\epsilon)$ is stochastic in a
neighborhood of 0. Then, up to within $O(\epsilon^{2})$ terms we have

$\left|\kappa(P(\epsilon))-\kappa(P_{1})-\kappa(P_{2})-\frac{1}{2}\right|\leq\epsilon\|(I-P+\operatorname{\mathbf{1}}\mathbf{h}^{T})^{-2}\|_{F}\|E\|_{F}+O(\epsilon^{2}),$

where $\mathbf{h}$ is any vector such that
$\mathbf{h}^{T}\operatorname{\mathbf{1}}=1$

Similar bounds can be obtained for the “stochastic” perturbation of a cyclic
matrix having arbitrary cyclicity index and for the perturbation of the
Kronecker product of two stochastic matrices.

Analogously, we may provide a “perturbed version” of Proposition 4.8 obtained
by combining Proposition 4.8 itself with Equation (35).

##  6 A divide-and-conquer algorithm

We can use the analysis done in Section 3 to construct a divide-and-conquer
algorithm for computing Kemeny’s constant of a stochastic sparse matrix $P$.
In fact, it is sufficient to identify a partitioning of the form (1) to start
a recursion procedure. This procedure employs Theorem 3.3 to express
$\kappa(P)$ in terms of $\kappa(P_{1})$ and $\kappa(P_{2})$ of the censored
chains and the $\gamma$, where $P_{1}$ and $P_{2}$ are the stochastic
complements. Algorithm 1 compactly describes the entire procedure in a
recursive formulation.

1 _Kemeny(_$P$ ,$\boldsymbol{\pi}$_)_ /* Implementation as a recursive
function */

Input: Stochastic matrix $P$, stationary distribution vector
$\boldsymbol{\pi}$, an integer $n_{0}>0$.

Output: Kemeny’s constant $\kappa(P)$ of $P$.

2 /* In the following the $\leftarrow$ means variable assignment while the $=$
are used to define shorthand. */

3 $n\leftarrow\operatorname{size}(P)$;

4 if  _$n <n_{0}$_ then

5
$\kappa(P)\leftarrow\operatorname{Tr}\left((I-P+\nicefrac{{\mathbf{1}\mathbf{1}^{T}}}{{n}})^{-1}\right)-1$;

6 return

7 else

8 $m\leftarrow\lfloor\nicefrac{{n}}{{2}}\rfloor$;

9 $P_{11}=P(1:m,1:m)$;

10 /* Block structure */

11 $P_{12}=P(1:m,m+1:n)$;

12 $P_{21}=P(m+1:n,1:m)$;

13 $P_{22}=P(m+1:n,m+1:n)$;

14 $\boldsymbol{\pi}_{1}=\boldsymbol{\pi}(1:m)$;

15 $\boldsymbol{\pi}_{2}=\boldsymbol{\pi}(m+1:n)$;

16
$\hat{}\boldsymbol{\pi}_{i}\leftarrow\frac{1}{\|\boldsymbol{\pi}_{i}\|}\boldsymbol{\pi}_{i},~{}~{}i=1,2$;

17 /* Remark: we can compute the $LU$ factorization of $I-P_{22}$ once and
reuse it. */

18 $P_{1}\leftarrow P_{11}+P_{12}(I-P_{22})^{-1}P_{21}$;

19 /* stochastic complements */

20 $P_{2}\leftarrow P_{22}+P_{21}(I-P_{11})^{-1}P_{12}$;

21 $\mathbf{x}\leftarrow(I-P_{22})^{-1}\operatorname{\mathbf{1}}$;

22 /* Computation of $\theta$ using (17) */

23
$\mathbf{y}\leftarrow(I-P_{1}+\operatorname{\mathbf{1}}\hat{\boldsymbol{\pi}}_{1}^{T})^{-1}(P_{12}\mathbf{x})$;

24 $\mathbf{y}\leftarrow(I-P_{22})^{-1}(P_{21}\mathbf{y})$;

25 $\theta\leftarrow\hat{\boldsymbol{\pi}}_{2}^{T}(\mathbf{x}+\mathbf{y})$;

26 $\gamma\leftarrow\|\boldsymbol{\pi}_{1}\|\theta-\|\boldsymbol{\pi}_{2}\|$;

27 /* Compute correction as in (16) */

28 $\kappa(P_{1})\leftarrow\;$Kemeny(_$P_{1}$ ,$\hat{}\boldsymbol{\pi}_{1}$_);

29 /* Recursion */

30 $\kappa(P_{2})\leftarrow\;$Kemeny(_$P_{2}$ ,$\hat{}\boldsymbol{\pi}_{2}$_);

31 /* Recursion */

32 $\kappa(P)\leftarrow\kappa(P_{1})+\kappa(P_{2})+\gamma$;

33 return

34

Algorithm 1 Divide-and-conquer algorithm for the computation of Kemeny’s
constant.

To effectively apply this strategy, it is crucial to efficiently address
several computational subproblems. Let us set aside the computation of the
stationary distribution vector $\boldsymbol{\pi}$ for the starting chain $P$,
which is the essential component required to initiate the entire procedure (to
this regard, for large-scale problems, the algorithms proposed in [13] might
be used). The most significant part of the computation lies in the solution of
the linear systems to lines 1–1 of Algorithm 1.

Let us focus on the solution of the following systems:

$(I-P_{ii})\mathbf{x}=\mathbf{b},\qquad i=1,2,$

where we further assume that the block matrices $\\{P_{ii}\\}_{i=1}^{2}$ can
benefit from sparse storage. By construction, the blocks $P_{ii}$, $i=1,2$ are
sub-stochastic, i.e., $(P_{ii})_{p,q}\geq 0$ for all $p,q$, and
$P_{ii}\operatorname{\mathbf{1}}\leq\operatorname{\mathbf{1}}$,
$P_{ii}\operatorname{\mathbf{1}}\neq\operatorname{\mathbf{1}}$. Since $P$ is
irreducible, the matrices $I-P_{ii}$, $i=1,2$, are non-singular M-matrices.
This property allows us to resort to different efficient iterative strategies
for the solution of the systems involved. Since in general matrices will be
non-symmetric, one can consider using GMRES [33, Section 6.5] or BiCGstab [33,
Section 7.4.2] for solving the different linear systems. (In any case it will
be necessary to have a preconditioner available to accelerate the convergence
of the Krylov method in question.) Since we are working with M-matrices, a
natural choice is to use incomplete factorizations. Specifically, we can use
Incomplete LU factorizations (ILU), that is, we can approximate the matrices
as

$I-P_{ii}=\widetilde{L}_{i}\widetilde{U}_{i}+R_{i}\approx\widetilde{L}_{i}\widetilde{U}_{i}\qquad
i=1,2,$

with the residual matrix $R_{i}$ satisfying certain constraints, such as
having zero entries in some prescribed locations–either static, determined on
the base of the natural occurring fill-in during the computation, or via
thresholding on their entries. For M-matrices, the existence of such objects
is guaranteed by the fact that Gaussian elimination and non-diagonal dropping
of the entries preserves the property of being a non singular M-matrix–see
[15] for the original proof or [33, Theorem 10.1] for a modern explanation.
Similarly, to precondition the system on the line 1 of the algorithm that
contains the matrix $P_{1}$, we can consider an incomplete factorization of an
approximation of the matrix $I-P_{1}$, namely,

$\begin{split}I-P_{1}=&\;I-P_{11}-P_{12}(I-P_{22})^{-1}P_{21}\\\
\approx&\;I-P_{11}-P_{12}(\operatorname{diag}(I-P_{22}))^{-1}P_{21}\approx\widetilde{L}_{3}\widetilde{U}_{3},\end{split}$

where $\operatorname{diag}(I-P_{22})$ is the diagonal matrix formed by the
diagonal entries of $I-P_{22}$. It is known that the computation of LU
factorization, incomplete or not, benefits from the reordering of the entries
of the matrix itself, see, e.g., [2]. Furthermore, as seen in Proposition 4.8,
the closer a matrix is to block diagonal after appropriate permutation, the
easier the global calculation will be. Since the off-diagonal matrices in the
block decomposition (1) will be made up of a few nonzero elements, the natural
choice for the permutation algorithm is to use Nested Dissection permutation
algorithm [20]; an example of the application of this algorithm is given in
Figure 1, from which we observe that the resulting matrix has the desired
“quasi block-diagonal” structure.

Figure 1: Application of the Nested Dissection algorithm from [20] to the
Gaertner/nopoly matrix from the SuiteSparse Matrix Collection (formerly the
University of Florida Sparse Matrix Collection) [12].

We want to underline that this initial permutation step would also be
advisable if one wants to obtain Kemeny’s constant directly from the
expression (6), since the related effect of this choice is also that of a
permutation which is fill-reducing for the computation of $LU$ factorization
of sparse matrices, a step needed for the efficient computation of the matrix
inverse in (6).

###  6.1 Low-precision randomized approximation

We note that since we want to solve intermediate linear systems using an
iterative method, what we are actually calculating is an approximation of
Kemeny’s constant. For this reason, it makes sense to consider randomized
algorithms for the direct approximation of (6) through the trace of a matrix.
Such approaches have been also used for instance in [30], for reversible
Markov chains, and in [35], for Markov chains modeling a random walk on an
undirected graphs. We consider here the application of the Hutch++ algorithm
[32], which is a straightforward improvement on the Hutchinson estimator [18]
used in [35]. In our case this means having an _oracle_ that computes

$\mathbf{y}=(I-P+\operatorname{\mathbf{1}}\mathbf{h}^{T})^{-1}\mathbf{x}\equiv
A\mathbf{x},\quad\mathbf{h}=\mathbf{1}/n.$ (36)

Then the following results gives us the number of oracle calls, i.e., linear
system solutions, we have to approximate
$\operatorname{Tr}\left((I-P+\operatorname{\mathbf{1}}\mathbf{h}^{T})^{-1}\right)$
within a given tolerance.

######  Theorem 6.1 ([32, Theorem 1]).

If Hutch++ is implemented with
$\ell=O(\nicefrac{{\sqrt{\log(\nicefrac{{1}}{{\delta}})}}}{{\epsilon}}+\log(\nicefrac{{1}}{{\delta}}))$
matrix-vector multiplication queries, then for any positive semidefinite
matrix $A$, with probability $\geq 1-\delta$, the output of
$\text{Hutch\raisebox{1.50696pt}{++}}(A)$ satisfies

$(1-\varepsilon)\operatorname{Tr}(A)\leq\text{Hutch\raisebox{1.50696pt}{++}}(A)\leq(1+\varepsilon)\operatorname{Tr}(A).$

Similarly to what was done for the solution of linear systems in the divide-
and-conquer algorithm, we can use the PCG with an Incomplete Cholesky
preconditioner calculated on $I-P$, i.e.,

$A\approx\widetilde{L}^{-T}\left(I-\frac{1}{1+(\mathbf{h}^{T}\widetilde{L}^{-T})(\widetilde{L}^{-1}\operatorname{\mathbf{1}})}(\widetilde{L}^{-1}\operatorname{\mathbf{1}})(\mathbf{h}^{T}\widetilde{L}^{-T})\right)\widetilde{L}^{-1},$

as an oracle for the calculation of the products $A\mathbf{x}$ necessary for
Hutch++. We point out that the Hutch++ algorithm works under the assumption
that the matrix of which we approximate the trace is symmetric positive
semidefinite, and this assumption is verified for Markov chains modeling a
random walk on an undirected graph.

###  6.2 Numerical examples

This section contains some numerical examples in which the performance of the
algorithms obtained starting from the theoretical analysis is analyzed. All
experiments are reproducible starting from the code contained in the
repository github.com/Cirdans-Home/Kemeny-and-Conquer. All experiments are
performed on a vertex of the Toeplitz cluster at the University of Pisa
equipped with an Intel® Xeon® CPU E5-2650 v4 at 2.20GHz and 250 Gb of RAM,
using MATLAB 9.10.0.1602886 (R2021a). The Markov chains used in the examples
are built employing matrices from the SuiteSparse Matrix Collection (formerly
the University of Florida Sparse Matrix Collection) [12]. Specifically, we
build the probability transition matrix $P$ from irreducible adjacency
matrices $A$ as

$P=\operatorname{diag}(\hat{A}\operatorname{\mathbf{1}})^{-1}\hat{A},$

for $\hat{A}$ the matrix obtained from $A$ with all weights set to $1$. For
the low-precision randomized case we use instead

$P=\operatorname{diag}(\hat{A}\operatorname{\mathbf{1}})^{-\nicefrac{{1}}{{2}}}\hat{A}\operatorname{diag}(\hat{A}\operatorname{\mathbf{1}})^{-\nicefrac{{1}}{{2}}},$

for cases with $A=A^{T}$, i.e., the graph is undirected. All the relative
errors in the following experiments are computed with respect to Kemeny’s
constant obtained directly, i.e., applying (6) by computing the whole matrix
inverse.

####  6.2.1 Low-precision randomized approximation

We first focus on using randomized estimators from Section 6.1 for the trace
of a matrix. Table 1 contains the estimates obtained using Hutch++ with the
parameters $(\delta,\epsilon)$ of Theorem 6.1 chosen as
$\delta=\nicefrac{{1}}{{4}}$, and $\epsilon=10^{-1}$, i.e., we use a set of
$l=13$ random sample vectors.

|  | Time (s) |
---|---|---|---
Matrix | $n$ | Direct | Hutch++ | Rel. Error
Pajek/USpowerGrid | 4941 | 2.11 | 0.43 | 2.69e-04
Gaertner/nopoly | 10774 | 21.97 | 2.64 | 7.58e-03
Gleich/minnesota | 2640 | 0.33 | 0.15 | 1.94e-03
Table 1: Randomized approximation of Kemeny’s constant with Hutch++algorithm.
The relative error and the timings (s) for the Hutch++is the average over 100
repetitions. The inner iterative solver is PCG preconditioned by the
ICHOL(0)-based preconditioner from 6.1 and with a tolerance of $10^{-2}$. The
number of samples is $l=13$, that corresponds to $\delta=1/4$ and
$\epsilon=10^{-1}$.

As an internal solver for the oracle calculation we use the PCG preconditioned
with an ICHOL(0), i.e., such as to preserve the sparsity pattern of the
starting matrix in $\widetilde{L}$. Since the external precision we expect to
achieve is of the order of $10^{-2}$, the tolerance for the iterative method
is chosen to be $10^{-3}$. In all cases, we always start from the matrix
reordered with the Nested Dissection permutation algorithm [20].

####  6.2.2 Divide-and-conquer algorithm

Here we test the application of the divide-and-conquer strategy on the
transition matrices built from the matrices in SuiteSparse Matrix Collection
[12] in the previous section. We use two variants of Algorithm 1: one called
Recursive in Table 2 exploits an incomplete ILU(0) factorization and the GMRES
to solve systems with left term $I-P_{22}$, and the other called direct-
recursive exploits the LU factorization of the matrix, already needed for the
calculation of the stochastic complement, also for the solution of the two
auxiliary systems with the same matrix for the calculation of $\theta$.

|  |  | Time (s) | Rel. Error
---|---|---|---|---
Matrix | $n$ | $\kappa(P)$ | Direct | Recursive | Dir-Rec | Recursive | Dir-Rec
Gaertner/big | 13209 | 58134.53 | 17 | 5.81 | 4.15 | 1.57e-08 | 4.09e-09
vanHeukelum/cage10 | 11397 | 15378.12 | 11.06 | 44.83 | 56.76 | 8.28e-13 | 1.25e-11
vanHeukelum/cage11 | 39082 | 51177.08 | 315.27 | 1259.11 | 1591.04 | 2.59e-10 | 1.53e-10
HB/gre_1107 | 1107 | 1483.57 | 0.04 | 0.10 | 0.11 | 1.97e-09 | 1.09e-09
Gaertner/nopoly | 10774 | 171656.87 | 9.56 | 9.14 | 6.71 | 3.25e-07 | 3.09e-07
Gaertner/pesa | 11738 | 131250.78 | 11.45 | 6.36 | 3.30 | 2.15e-07 | 2.28e-07
Gleich/usroads-48 | 126146 | 1818057.53 | 8243.57 | 743.37 | 542.88 | 8.65e-07 | 1.60e-07
Barabasi/NotreDame_www$\dagger$ | 34643 | 1173610.94 | 172.36 | 94.28 | 93.20 | 2.85e-04 | 3.58e-05
Pajek/USpowerGrid | 4941 | 30166.55 | 1.32 | 1.48 | 0.58 | 3.01e-08 | 2.80e-08
Gleich/minnesota$\dagger$ | 2640 | 18243.53 | 0.30 | 0.34 | 0.22 | 8.37e-08 | 3.53e-08
Table 2: Performance of the recursive implementation of the divide-and-conquer
algorithm for computing Kemeny’s constant on some test matrices. If the matrix
name has a $\dagger$ symbol, then the experiment has been run on the largest
connected component of the graph, i.e, on the largest irreducible sub-chain.

From the results in Table 2 we observe that in most cases the divide-and-
conquer algorithm manages to reduce the computational time compared to the
direct computation of Kemeny’s constant. In some cases we observe the absence
of an improvement, investigating in detail what we observe is that the
decomposition into blocks done by halving is far from being optimal. This
causes both the creation of denser stochastic complements and higher solution
times for auxiliary linear systems. The implementation of nested dissection in
MATLAB does not have in output the limitation of the _clusters_ obtained,
being able to use that should increase the advantage of the recursive version
compared to the one in which the recursion is done by simple halving.

##  7 Conclusions

Kemeny’s constant $\kappa(P)$ of a stochastic matrix $P$ has been expressed in
terms of Kemeny’s constants of the stochastic complements obtained from a
block partitioning, and the constant $\gamma$. Explicit expressions of
$\kappa(P)$ have been provided for the transition matrix of a periodic Markov
chain, the Kronecker product of stochastic matrices, and sub-stochastic
matrices with constant row sums. The main result, Theorem 3.3, has been used
to design a divide-and-conquer algorithm for recursively computing
$\kappa(P)$. Numerical experiments applied to real world graphs show the
effectiveness of this approach especially in the case of nearly completely
decomposable matrices.

As Kemeny’s constant measures the expected time of a Markov chain $X$ to
travel between two randomly chosen states, a natural question arises: “What
interpretation can be ascribed to Kemeny’s constant of a censored Markov chain
$X_{1}$ associated with the original Markov chain $X$?”. Since $X_{1}$ is
induced from $X$, it would be interesting to provide the insights into what
features of $X$ are captured in Kemeny’s constant of $X_{1}$. Partitioning the
state space of $X$ into two subsets yields two censored Markov chains. We can
see from Theorem 3.3 that Kemeny’s constants of these censored Markov chains
are interdependent. Understanding how partitioning the state space influences
Kemeny’s constants of censored Markov chains would be beneficial for gaining
insights into the question.

## References

  * [1] D. Altafini, D. A. Bini, V. Cutini, B. Meini, and F. Poloni.  An edge centrality measure based on the Kemeny constant.  SIAM J. Matrix Anal. Appl., 44(2):648–669, 2023.
  * [2] M. Benzi, D. B. Szyld, and A. van Duin.  Orderings for incomplete factorization preconditioning of nonsymmetric problems.  SIAM J. Sci. Comput., 20(5):1652–1670, 1999.
  * [3] D. Bini, J. J. Hunter, G. Latouche, B. Meini, and P. Taylor.  Why is Kemeny’s constant a constant?  Journal of Applied Probability, 55(4):1025–1036, 2018.
  * [4] A. Brauer.  Limits for the characteristic roots of a matrix. VII.  Duke Math. J., 25:583–590, 1958.
  * [5] J. Breen, S. Butler, N. Day, C. DeArmond, K. Lorenzen, H. Qian, and J. Riesen.  Computing Kemeny’s constant for barbell-type graphs.  Electron. J. Linear Algebra, 35:583–598, 2019.
  * [6] J. Breen, N. Faught, C. Glover, M. Kempton, A. Knudson, and A. Oveson.  Kemeny’s constant for nonbacktracking random walks.  Random Structures Algorithms, 63(2):343–363, 2023.
  * [7] J. Breen and S. Kirkland.  A structured condition number for Kemeny’s constant.  SIAM J. Matrix Anal. Appl., 40(4):1555–1578, 2019.
  * [8] P. Brémaud.  Markov chains: Gibbs fields, Monte Carlo simulation, and queues, volume 31.  Springer Science & Business Media, 2001.
  * [9] L. Ciardo.  The Braess’ paradox for pendent twins.  Linear Algebra and its Applications, 590:304–316, 2020.
  * [10] L. Ciardo, G. Dahl, and S. Kirkland.  On Kemeny’s constant for trees with fixed order and diameter.  Linear and Multilinear Algebra, 70(12):2331–2353, 2022.
  * [11] E. Crisostomi, S. Kirkland, and R. Shorten.  A google-like model of road network dynamics and its application to regulation and control.  International Journal of Control, 84(3):633–651, 2011.
  * [12] T. A. Davis and Y. Hu.  The University of Florida sparse matrix collection.  ACM Trans. Math. Software, 38(1):Art. 1, 25, 2011.
  * [13] H. De Sterck, K. Miller, E. Treister, and I. Yavneh.  Fast multilevel methods for Markov chains.  Numer. Linear Algebra Appl., 18(6):961–980, 2011.
  * [14] P. G. Doyle.  The Kemeny constant of a Markov chain.  arXiv preprint arXiv:0909.2636, 2009.
  * [15] K. Fan.  Note on $M$-matrices.  Quart. J. Math. Oxford Ser. (2), 11:43–49, 1960.
  * [16] N. Faught, M. Kempton, and A. Knudson.  A 1-separation formula for the graph Kemeny constant and Braess edges.  Journal of Mathematical Chemistry, 60(1):49–69, 2022.
  * [17] R. A. Horn and C. R. Johnson.  Matrix analysis.  Cambridge University Press, Cambridge, second edition, 2013.
  * [18] M. F. Hutchinson.  A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines.  Comm. Statist. Simulation Comput., 19(2):433–450, 1990.
  * [19] J. Jang, S. Kim, and M. Song.  Kemeny’s constant and Wiener index on trees.  Linear Algebra Appl., 674:230–243, 2023.
  * [20] G. Karypis and V. Kumar.  A fast and high quality multilevel scheme for partitioning irregular graphs.  SIAM J. Sci. Comput., 20(1):359–392, 1998.
  * [21] J. G. Kemeny.  Generalization of a fundamental matrix.  Linear Algebra Appl., 38:193–206, 1981.
  * [22] J. G. Kemeny and J. L. Snell.  Finite Markov chains.  Undergraduate Texts in Mathematics. Springer-Verlag, New York-Heidelberg, 1976.  Reprinting of the 1960 original.
  * [23] S. Kim.  Families of graphs with twin pendent paths and the Braess edge.  The Electronic Journal of Linear Algebra, pages 9–31, 2022.
  * [24] S. Kim, J. Breen, E. Dudkina, F. Poloni, and E. Crisostomi.  On the effectiveness of random walks for modeling epidemics on networks.  Plos one, 18(1):e0280277, 2023.
  * [25] S. Kim, N. Madras, A. Chan, M. Kempton, S. Kirkland, and A. Knudson.  Bounds on Kemeny’s constant of a graph and the Nordhaus-Gaddum problem.  arXiv preprint arXiv:2309.05171, 2023.
  * [26] S. Kirkland.  Fastest expected time to mixing for a Markov chain on a directed graph.  Linear Algebra and its Applications, 433(11-12):1988–1996, 2010\.
  * [27] S. Kirkland.  Directed forests and the constancy of Kemeny’s constant.  Journal of Algebraic Combinatorics, 53(1):81–84, 2021.
  * [28] S. Kirkland, Y. Li, J. McAlister, and X. Zhang.  Edge addition and the change in Kemeny’s constant.  arXiv preprint arXiv:2306.04005, 2023.
  * [29] S. Kirkland and Z. Zeng.  Kemeny’s constant and an analogue of Braess’ paradox for trees.  Electronic Journal of Linear Algebra, 31:444–464, 2016.
  * [30] S. Li, X. Huang, and C.-H. Lee.  An Efficient and Scalable Algorithm for Estimating Kemeny’s Constant of a Markov Chain on Large Graphs.  In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, KDD ’21, page 964–974, New York, NY, USA, 2021. Association for Computing Machinery.
  * [31] C. D. Meyer.  Stochastic complementation, uncoupling Markov chains, and the theory of nearly reducible systems.  SIAM Review, 31(2):240–272, 1989.
  * [32] R. A. Meyer, C. Musco, C. Musco, and D. P. Woodruff.  Hutch++: optimal stochastic trace estimation.  In Symposium on Simplicity in Algorithms (SOSA), pages 142–155. [Society for Industrial and Applied Mathematics (SIAM)], Philadelphia, PA, 2021.
  * [33] Y. Saad.  Iterative methods for sparse linear systems.  Society for Industrial and Applied Mathematics, Philadelphia, PA, second edition, 2003.
  * [34] X. Wang, J. L. Dubbeldam, and P. Van Mieghem.  Kemeny’s constant and the effective graph resistance.  Linear Algebra and its Applications, 535:231–244, 2017.
  * [35] W. Xu, Y. Sheng, Z. Zhang, H. Kan, and Z. Zhang.  Power-law graphs have minimal scaling of Kemeny constant for random walks.  page 46 – 56, 2020.
  * [36] S. Yilmaz, E. Dudkina, M. Bin, E. Crisostomi, P. Ferraro, R. Murray-Smith, T. Parisini, L. Stone, and R. Shorten.  Kemeny-based testing for COVID-19.  Plos One, 15(11):e0242401, 2020.