# Learning the Influence Graph of a High-Dimensional Markov Process with
Memory

Smita Bagewadi and Avhishek Chatterjee Authors are with the Department of
Electrical Engineering, Indian Institute of Technology Madras, Chennai,
India.AC thanks the Department of Science and Technology, India, for support
through the grant INSPIRE/04/2016/001171.

###### Abstract

Motivated by multiple applications in social networks, nervous systems, and
financial risk analysis, we consider the problem of learning the underlying
(directed) influence graph or causal graph of a high-dimensional multivariate
discrete-time Markov process with memory. At any discrete time instant, each
observed variable of the multivariate process is a binary string of random
length, which is parameterized by an unobservable or hidden $[0,1]$-valued
scalar. The hidden scalars corresponding to the variables evolve according to
discrete-time linear stochastic dynamics dictated by the underlying influence
graph whose nodes are the variables. We extend an existing algorithm for
learning i.i.d. graphical models to this Markovian setting with memory and
prove that it can learn the influence graph based on the binary observations
using logarithmic (in number of variables or nodes) samples when the degree of
the influence graph is bounded. The crucial analytical contribution of this
work is the derivation of the sample complexity result by upper and lower
bounding the rate of convergence of the observed Markov process with memory to
its stationary distribution in terms of the parameters of the influence graph.

###### Index Terms:

Causal graph; Markov process; Sample complexity

##  I Introduction

Learning the underlying influence graph, a.k.a. causal graph, of a high
dimensional multivariate Markov process is an intriguing mathematical problem
with many applications [24, 2, 8, 30]. It is central to many disciplines,
including but not limited to social network analysis, neural and biological
signal processing, and financial risk analysis. In this paper, we study the
problem of learning the influence graph of the following high dimensional
$d$-Markov process, i.e., a Markov process with $d$ memory, from its samples.

There is an underlying directed graph $G$ with nodes $V$ and directed edges
$E$. At a given time $t$, any node $v\in V$ has an internal parameter
$X_{v}(t)\in[0,1]$ and generates multiple (possibly random number of)
independent instances of Bernoulli$(X_{v}(t))$ random variables. These
Bernoulli variables at $v$ over the last $d$ time steps, in turn, impact the
future $X_{u}(t+1)$ of node $u$ through weighted averaging if there is a
directed edge $(v,u)\in E$ from $v$ to $u$. Our problem is to learn the graph
$G$ by observing only the instances of the Bernoulli variables across all
nodes.

Below, we discuss some application areas where the above problem is often
encountered.

Social networks: In many online social platforms, agents express their
opinions in binary form as likes or dislikes. The expressed opinions on the
platform influence the internal or private beliefs of other agents. In
particular, an agent’s internal belief is a weighted average of the recently
expressed opinions of the agents it regularly follows [1],[2],[3],[4],[5]. In
turn, the internal or private belief of an agent dictates the nature of its
expressed opinions. Thus, the process of expressed opinions of all agents in a
social network is indeed a high-dimensional $d$-Markov process of the above
kind. It is known that only a few agents on the friend list of an agent on a
social platform have a real influence on it, and often, an agent not on the
friend list has an influence either through the physical world or through
other platforms. Hence, one cannot directly use the friend list as the
influence graph; rather, one has to learn the true influencers using the high-
dimensional process of expressed binary opinions.

Brain and biological networks: In neuroscience, there has been significant
interest in finding the hidden logical connections between different centers
of the brain and nervous systems [7], [8], [9]. The only observable data are
the EEG or other measurements of different brain and nerve centers. Neuronal
triggers or firings have inherent casuality and Markovianity [29]. The firing
rate of a nerve center, i.e., the number of (binary) spikes in a time interval
of appropriate duration, reflects the internal excitement of that center.
These spikes then contribute to the excitement of other nervous centers
causally and cumulatively [29], [28]. Thus, the overall process of neuronal
triggers across many centers behaves like a Markov process with binary
observables. Similar dynamics are also observed in other biological networks,
for example, gene regulatory networks [6].

Financial risk analysis:  Major financial changes in a company are often
influenced by other companies and worldwide financial and administrative
bodies, some of which may not even have any direct financial link with the
company. Thus, for accurate financial risk analysis, learning such hidden
connections is extremely important [19],[20]. In practice, such connections
have to be learned based on the raw data of the performance of different
companies and the actions of different bodies in the past. This, in essence,
is a problem of learning the causal graph of a high-dimensional Markov
process. Though the above-mentioned binary model is not a perfect
representation of this scenario, it is a tractable first-order approximation.

###  I-A Main Contribution

Under the condition that the in-degree of any node is $O(1)$, we present a
conditional directed entropy-based polynomial-time greedy algorithm for
learning the causal influence graph based on the binary observations. For
provable recovery of the true causal graph with high probability, we provide
an upper bound on sample complexity that scales logarithmically with the
number of nodes. The greedy algorithm presented here is an extension of a
greedy algorithm for learning i.i.d. graphical models [13] to our Markovian
dynamics.

Our main mathematical contribution lies in the derivation of a sample
complexity bound. An important component of our proof technique builds on a
relatively recent and tighter result regarding a Hoeffding-style Markov
concentration [27] using the second eigenvalue. The main novelty of our proof
technique lies in establishing an upper bound on the second eigenvalue in
terms of the parameters of the Markov dynamics. We achieve this by upper
bounding the total variation distance from the stationary distribution in
terms of the parameters of the dynamics and comparing that with a lower bound
obtained in terms of the second eigenvalue.

###  I-B Organization

The rest of the paper is organized as follows. In Sec. II, we discuss existing
literature that shares some connection with our work. We present the
$d$-Markov model in detail in Sec. III. In Sec. IV, we discuss some useful
mathematical structures of the model. Next, we present the algorithm and its
performance guarantees in Sec. V followed by the detailed proofs in Sec. VI.
In Sec. VII, we present some simulation results that corroborate and
complement our analytical bounds. Finally, we conclude in Sec. VIII.

##  II Related Literature

To the best of our knowledge, no analytical work has been reported in the
literature on the particular mathematical problem considered here. However,
our problem has partial but interesting overlaps with multiple areas of
theoretical machine learning and applied mathematics. From the point of view
of applications, as discussed in Sec. I, this work is related to learning
influence in social networks and learning logical connections in nervous
systems. Below, we discuss some of the related literature, starting with
application-oriented works.

###  II-A Data Driven Applications

In [1], the authors learn the strength of the influence in a social network
graph on each edge, assuming continuous time linear opinion evolution and the
knowledge of the structure of the influence graph, given only the expressed
opinions. An agent’s opinion at a time is measured as a real number by using
natural language processing based sentiment analysis of the posted texts. A
discrete-time continuous-valued opinion evolution model was considered in [3],
[2], and the underlying influence strengths were learned by posing well-
behaved optimization problems. In [3], the authors exploit the linear
structure to learn the influences using only opinions at time $t=0$ and
$\infty$ across multiple sample paths. On the other hand, influences are
learned in [2] using one sample path by exploiting the underlying sparsity of
the influence graph.

Over the last decade or so, there has been interest in understanding causal
influence patterns between different parts of the brain or nervous system
while performing a specific cerebral or motor task [7, 8, 9, 29]. The problem
of understanding brain communities using temporal information flow in the
brain was studied in [7]. Zhang et al. worked on estimating the functional
connectivity in a brain using a sparse hidden Markov model [9]. Quinn et al.
used estimates of directed information to understand causal relationships
between different neural spike trains from different brain centers [8].

The similarity between genes having the same phenotype was learned using
random walks on a network of interacting genes [6]. On the other hand,
learning the underlying causal influence graph for financial networks has been
undertaken in [19, 20].

###  II-B Analytical Guarantees

Learning the underlying graphical model of a Markov Random Field (MRF) has
been a flourishing area of research in the last decade. Extensive research has
been done on learning graphical models using i.i.d. samples from the MRF, for
example, [10, 11, 12, 14, 15, 16]. One of the major challenges in learning
graphical models was the locally diamond-like structure [12] and the
assumption on correlation decay. These were settled in [13], which proved
efficient learning in the presence of a diamond-like structure and without
invoking correlation decay. However, the non-degeneracy assumption for
conditional entropy was necessary. Further, in [32], algorithms were presented
for learning antiferromagnetic Ising models, where there is no correlation
decay.

Learning causality in a multi-variate process has a rich history [21, 23, 22,
24, 25]. Directed mutual information was used to learn causality between two
stationary and ergodic processes in [8]. Learning the underlying parameters of
two causally related Markov processes have been studied in [18].

In [17], authors provide sample complexity guarantees for learning the
influence graph of a special Markov process, the Bernoulli autoregressive
(BAR) process. The binary random variables associated with the nodes are
Bernoulli with probability derived from the past states of their neighbors.
Though [17] has quite a few similarities with our work, there are some major
differences as well: (i) in our problem, each node generates multiple
(possibly random number) of Bernoulli variables at each time, (ii) we do not
assume the knowledge of the degree of each node and a lower bound on the
stationary probabilities and (iii) we use directed conditional entropy based
greedy algorithms.

The work closest to the current paper is [24], where learning the causal graph
of a general discrete-time Markov chain is considered. However, in [24], many
assumptions had to be imposed on the Markov process to obtain sample
complexity guarantees. These assumptions include a lower bound on the smallest
value of the stationary distribution, the absence of a diamond-like directed
structure, and correlation decay-like assumption for Markov processes. By
using ideas for learning i.i.d. graphical model from [13], some recent results
in concentration inequalities [27] and an approach developed by us for
bounding the second eigenvalue of the particular Markov chain in question, we
could relax those assumptions. Also, unlike the bound in [24], which is in
terms of the mixing time, our bound is directly in terms of the edge weights
and the graph structure.

##  III Model

The pattern of influence is represented by a directed graph $G=(V,E)$. The
relationship between the nodes $u,v\in V$ is represented by the edge $(u,v)\in
E$. The edge $(u,v)\in E$ implies that node $u$ influences node $v$. We
consider directed edges for capturing asymmetry in the influence, i.e., for
$u\neq v$, $(u,v)\neq(v,u)$. The edge $(v,v)$ captures the self-influence of
$v$ and is assumed to be present for each node. The neighborhood
$\mathcal{N}_{v}$ of node $v$ is defined as the set of nodes that influence
$v$, i.e., $\mathcal{N}_{v}:=\\{u\neq v:(u,v)\in E\\}$.

We consider a setting where a binary observed variable is evolving with time
over the above network. In particular, we consider a discrete-time system,
where at time $t$, node $v$ generates observation ‘$1$’ with probability
$X_{v}(t)$, where $X_{v}(t)$ is an internal parameter associated with the
node.

Based on the internal parameter value at time $t$, each node generates
$M_{v}(t)$ instances of the observed variable. We model $M_{v}(t)$ as a random
variable with distribution
$\min(\text{Poisson}(\mu_{v}(X_{v}(t))),\bar{M})+1$, where $\bar{M}$ is a
constant and for each $v$, $\mu_{v}(\cdot):[0,1]\to\mathbb{R}$ is a
$L$-Lipschitz function. Among these $M_{v}(t)$ binary random variables,
$N_{v}(t)$ are positive (i.e ‘1’), where
$N_{v}(t)\sim\text{Binomial}(M_{v}(t),X_{v}(t))$.

Corresponding to each edge $(u,v)$ there are $d$ non-negative weights
$a^{(r)}_{uv}$, $r=0,1,\ldots d-1$. These weights capture the strength of the
influence of the observed value of $u$ at $t-r$ on the internal parameter of
$v$ at $t$. For each $v$, the weights satisfy the condition
$\sum_{r=0}^{d-1}\sum_{u}a^{(r)}_{uv}=1$.

The evolution of the internal parameter of any node $v$ is affected by the
observed variables of the nodes in $\mathcal{N}_{v}$, its own inner bias
$l_{v}$, and the random fluctuations $Z_{v}(t)\in[0,1]$ with mean
$0<\bar{z}_{v}<1$. The internal parameters evolve according to the following
update equation. For any $v\in V$,

$\displaystyle X_{v}(t+1)$
$\displaystyle=(1-\alpha_{v})((1-\beta)Z_{v}(t+1)+\beta l_{v})$
$\displaystyle+\alpha_{v}\sum_{u\in\mathcal{N}_{v}\cup\\{v\\}}\sum_{r=0}^{(d-1)}a_{uv}^{(r)}\left(\frac{N_{u}(t)}{M_{u}(t)}\right)$
(1)

where $\alpha_{v},\beta\in(0,1)$.

Intuitively, the parameter $X_{v}(t)$ represents some internal process that is
affected by the generated instances of the observed variable. It evolves as a
weighted average of the past instances of the observed variable of its
neighbors, its own inner bias, and some unpredictable fluctuations. The
parameter $a_{uv}^{(r)}$ represents the weight given by node $v$ at time $t$
to the observed variables of node $u$ at time $t-r$. The parameter
$\alpha_{v}$ captures the openness of node $v$, i.e., the weight it gives to
the status of other nodes against its own inner bias and random fluctuations.

The learning problem: The central problem in this paper is to learn the
underlying influence graph based on the observed variables over a finite time
window. Mathematically, given $\\{(N_{v}(t),{M}_{v}(t)):v\in V\\}_{t=0}^{T-1}$
learn the underlying influence graph $G=(V,E)$, i.e., learn the neighborhood
$\mathcal{N}_{v}$ for all nodes $v\in V$.

##  IV Useful Mathematical Structures

Given the observation $\\{(N_{v}(t),M_{v}(t)):v\in V\\}_{t=0}^{T-1}$, we
construct the processes $Y_{v}(t):=\frac{N_{v}(t)}{M_{v}(t)}$,
$\mathbf{Y}^{(d)}_{v}(t):=\left(Y_{v}(t),Y_{v}(t-1),\ldots,Y_{v}(t-d+1)\right)$,
$\mathbf{Y}^{(d)}_{Q}(t):=\\{\mathbf{Y}^{(d)}_{v}(t):v\in Q\subseteq V\\}$,
$\mathbf{Y}^{(d)}_{u,Q}(t):=\\{\mathbf{Y}^{(d)}_{v}(t):v\in
Q\cup\\{u\\},Q\subseteq V\\}$ and
$\mathbf{Y}^{(d)}(t):=\\{\mathbf{Y}^{(d)}_{v}(t):v\in V\\}$. From the dynamics
in (III), it can be seen that
$\mathbf{Y}^{(d)}(t)=\\{\mathbf{Y}^{(d)}_{v}(t):v\in V\\}$ satisfies the
following lemma.

######  Lemma 1.

The process $\mathbf{Y}^{(d)}(t)$ is a Markov chain. Moreover, for any node
$v$ and at any time $t$, given
$\mathbf{Y}^{(d)}_{\mathcal{N}_{v}\cup\\{v\\}}(t)$,
$\mathbf{Y}^{(d)}_{v}(t+1)$ is independent of all other past process values
$\\{\mathbf{Y}^{(d)}_{V\setminus\\{\mathcal{N}_{v}\cup\\{v\\}\\}}(\tau):\tau\leq
t\\}$ and $\\{\mathbf{Y}^{(d)}_{\mathcal{N}_{v}\cup\\{v\\}}(\tau):\tau<t\\}$.

When $\mathbf{Y}^{(d)}(t)$ is stationary, for nodes $u,v\in V$,
$\mathcal{N}_{v}\subseteq Q$, and $u\notin Q$, the probability
$\mathbf{P}(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q})=\mathbf{P}(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q,u})$
for all $y_{v_{+}},\mathbf{y}^{(d)}_{v,Q,u}$, where $\mathbf{y}^{(d)}_{v,Q}$
represents the value of the process $\mathbf{Y}^{(d)}_{v,Q}(t)$ at the current
time instant and $y_{v_{+}}$ represents the value of the process $Y_{v}(t)$ at
the next time instant.

On the other hand, if $\mathcal{N}_{v}\not\subseteq Q$, then
$\mathbf{P}(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q})\neq\mathbf{P}(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q,u})$
for some $y_{v_{+}},\mathbf{y}^{(d)}_{v,Q,u}$. In fact, the following lower
bound on the total variation distance between these two conditional
distributions holds when $M_{v}(t)=1$ for all $v,t$

######  Lemma 2.

If $d=1$, $M_{v}(t)=1$ for all $t$ and $v$, and
$\mathcal{N}_{v}\setminus\\{u\\}\subset Q$ for some $Q\subset V$, then there
exists $y_{u}\in\\{0,1\\}$ such that for all values of
$\mathbf{y}^{(d)}_{v,Q}$
$||\mathbf{P}(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q},y_{u})-\mathbf{P}(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q})||_{1}$
is strictly bounded away from zero by a constant that does not depend on the
number of nodes.

###### Proof of Lemma 2.

We can write $\mathbf{P}(Y_{v_{+}}=1|\mathbf{Y}^{(d)}_{v,Q,u})$ as
$\mathbb{E}[X_{v_{+}}|\mathbf{Y}^{(d)}_{v,Q,u}].$ By using eq III, for
$M_{v}(t)=1,d=1,$ the quantity
$|\mathbf{P}(Y_{v_{+}}=1|\mathbf{Y}^{(d)}_{v,Q,u})-\mathbf{P}(Y_{v_{+}}=1|\mathbf{Y}^{(d)}_{v,Q})|$
equals

$\displaystyle\Big{\lvert}\alpha_{v}[a_{uv}(Y_{u}-\mathbb{E}[Y_{u}|Y_{v,Q}])+$
$\displaystyle\sum\limits_{w\in\mathcal{N}_{v}\setminus\\{v,Q,u\\}}\left(\mathbb{E}[Y_{w}|Y_{v,Q,u}]-\mathbb{E}[Y_{w}|Y_{v,Q}]\right)]\Big{\rvert}.$
(2)

When $\mathcal{N}_{v}\subseteq Q\cup\\{u\\}$, i.e., when $u$ is the only
neighbor not contained in $Q$, the set $\mathcal{N}_{v}\setminus\\{v,Q,u\\}$
is empty. Thus, in this case,

$\displaystyle|\mathbf{P}(Y_{v_{+}}=1|\mathbf{Y}_{v,Q,u})-\mathbf{P}(Y_{v_{+}}=1|Y_{v,Q})|$
$\displaystyle=\alpha_{v}a_{uv}|(Y_{u}-\mathbb{E}[Y_{u}|Y_{v,Q}])|$ (3)

$Y_{u}\in\\{0,1\\}$ and for any
$y_{v,Q}\in\\{0,1\\}^{|Q|+1},\mathbb{E}[Y_{u}|Y_{v,Q}=y_{v,Q}])\in(0,1).$ If
we consider the event $Y_{u}=0,$ we have

$\displaystyle|\mathbf{P}(Y_{v_{+}}=1|\mathbf{Y}_{v,Q},Y_{u}=0)-\mathbf{P}(Y_{v_{+}}=1|Y_{v,Q})|$
$\displaystyle=\alpha_{v}a_{uv}|\mathbb{E}[Y_{u}|Y_{v,Q}]|$ (4)

$Y_{u}$ is a Bernoulli with probability $X_{u}$,
$\mathbb{E}[Y_{u}|Y_{v,Q}=y_{v,Q}])=\mathbf{P}(Y_{u}=1|Y_{v,Q}=y_{v,Q})=\mathbb{E}[X_{u}|Y_{v,Q}=y_{v,Q}]$.
The quantity $\mathbb{E}[X_{u}|Y_{v,Q}=y_{v,Q}]$ is lower bounded by
$(1-\alpha_{v})((1-\beta)Z_{v}(t)+\beta l_{v}).$ Which means that when $u$ is
the only neighbor not included in $Q$, for all $\mathbf{y}_{v,Q},\exists\text{
a }y_{u}$ such that
$|\mathbf{P}(Y_{v_{+}}=1|\mathbf{y}_{v,Q})-\mathbf{P}(Y_{v_{+}}=1|\mathbf{y}_{v,Q},y_{u}=0)|>c_{1}$
where $c_{1}$ is a constant strictly greater than $0$

∎

The above result also seems to be true for $M_{v}(t)>1$ and $d>1$. However,
the proof remains elusive for the general case due to certain technicalities.
We make the following assumption for the general case.

######  Assumption 1.

Consider any $v\in V$, $u\in\mathcal{N}_{v}$, $Q\subset V$ and $u\not\in Q$.
Then for any choice of $y_{v}$ and $y_{Q}$, there exists $y_{u}$ such that
$\sum_{y_{v_{+}}}|\mathbf{P}(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q})-\mathbf{P}(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q,u})|>\sqrt{\epsilon^{\prime}},$
where $\epsilon^{\prime}>0$.

The algorithm presented here utilizes the conditional independence property
presented in Lemma 1 to determine the true neighborhoods of all the nodes. To
check these independence conditions, the algorithm uses the directed
conditional entropy of each node or variable of the process
$\mathbf{Y}^{(d)}(t)$.

The directed conditional entropy of node $v$ given the set $Q$ is given by:
$-\sum_{\mathbf{Y}^{(d)}_{Q,v}(t)}\sum_{\mathbf{Y}^{(d)}_{v}(t+1)}\mathbf{P}(\mathbf{Y}^{(d)}_{v}(t+1),\mathbf{Y}^{(d)}_{Q,v}(t))\log\mathbf{P}(\mathbf{Y}^{(d)}_{v}(t+1)|\mathbf{Y}^{(d)}_{Q,v}(t))$,
for some $Q\subset V$ and $v\not\in Q$. When the system is stationary, this
quantity would not depend on $t$, and we denote it as $H(v_{+}|v,Q)$. This can
be computed by computing the joint entropy $H(v_{+},v,Q)$ and $H(v,Q)$, where
$H(v_{+},v,Q)=-\sum_{\mathbf{Y}^{(d)}_{Q,v}(t),\mathbf{Y}^{(d)}_{v}(t+1)}\mathbf{P}(\mathbf{Y}^{(d)}_{v}(t+1),\mathbf{Y}^{(d)}_{Q,v}(t))\log\mathbf{P}(\mathbf{Y}^{(d)}_{v}(t+1),\mathbf{Y}^{(d)}_{Q,v}(t))$
and
$H(v,Q)=-\sum_{\mathbf{Y}^{(d)}_{Q,v}(t)}\mathbf{P}(\mathbf{Y}^{(d)}_{Q,v}(t))\log\mathbf{P}(\mathbf{Y}^{(d)}_{Q,v}(t))$.

Based on Lemma 1 and Assumption 1 (extension of Lemma 2), we obtain the
following useful fact.

######  Lemma 3.

For any $v$ and $u\not\in\mathcal{N}_{v}$, if $Q\subset V$ such that
$\mathcal{N}_{v}\subseteq Q$, then $H(v_{+}|v,Q,u)=H(v_{+}|v,Q)$. Further,
under Assumption 1, for any $Q\subset V$ and $u\in\mathcal{N}_{v}$ and
$u\not\in Q$, $H(v_{+}|v,Q)-H(v_{+}|v,Q,u)>\epsilon$, where $\epsilon>0$.

Proof of this lemma is presented in Sec. VI. Intuitively, the lemma says that
a node is a neighbor of $v$ only if adding that node to the conditioning set
$Q$ decreases the directed conditional entropy $H(v_{+}|v,Q)$ of node $v$. The
algorithms presented here build on this insight to learn the neighborhood
using directed conditional entropy. Note that when $Q$ does not contain
$\mathcal{N}_{v}$, Lemma 3 does not guarantee that every node $u$ which
exceeds the $\epsilon$ threshold is a true neighbor. The algorithm has to work
around this fact so that no node is falsely added to the neighborhood.

In this problem, however, $H(v_{+}|v,Q)$ cannot be directly computed since the
distributions are not known. Hence, they have to be estimated based on
samples, which we denote by $\hat{H}(v_{+}|v,Q)$. We first obtain plug-in
estimates the joint entropy $\hat{H}(v_{+},v,Q)$ and $\hat{H}(v,Q)$ using the
estimated joint distributions
$\hat{p}(\mathbf{y}^{(d)}_{v},\mathbf{y}^{(d)}_{Q})$ and
$\hat{p}(y_{v_{+}},\mathbf{y}^{(d)}_{v},\mathbf{y}^{(d)}_{Q})$ respectively.
Then, we subtract them to obtain the conditional entropy estimates
$\hat{H}(v_{+}|v,Q)$, which are used by the algorithms.

Building on [13], we observe that in our case, if the error in the estimation
of directed conditional entropy can be limited below a certain value, it is
possible to select a threshold that can distinguish between neighbors and non-
neighbors. This is stated mathematically in Lemma 4:

######  Lemma 4.

Under Assumption 1, if
$|\hat{H}(v_{+}|v,P)-H(v_{+}|v,P)|\leq\frac{\epsilon}{4}\text{ for all
}v\text{ and any }P\subseteq V$ , then for all $u\in\mathcal{N}_{v},v\in
V,\text{ and }Q\subset V\setminus\\{u,v\\}$,

$\hat{H}(v_{+}|v,Q)-\hat{H}(v_{+}|v,Q,u)>\frac{\epsilon}{2}$ (5)

and for all $w\notin\mathcal{N}_{v},Q\subset V\setminus\\{u,v,w\\}\text{ such
that }\mathcal{N}_{v}\subseteq Q$,

$\hat{H}(v_{+}|v,Q)-\hat{H}(v_{+}|v,Q,w)<\frac{\epsilon}{2}$ (6)

Next, we define a few quantities that would be useful later. For
$r=0,1,\ldots,d-1$, let $\tilde{A}(r)$ be a matrix whose $(v,u)^{th}$ element
$[\tilde{A}(r)]_{vu}=\alpha_{v}a^{(r)}_{uv}$. Let $\rho(\tilde{A}(r))$ be the
largest eigenvalue of $\tilde{A}(r)$ and let us define
$\rho(\tilde{A}):=\max\limits_{r}\rho(\tilde{A}(r))$. Also, let us define
$\bar{\mu}:=\max\limits_{v}\sup\limits_{x\in[0,1]}\mu_{v}(x)$.

##  V Algorithm

In this section, we present a greedy algorithm for recovering the underlying
directed graph $G$ from observations $\\{(N_{v}(t),M_{v}(t)):v\in
V\\}_{t=0}^{T-1}$, which is an adaptation of the RecGreedy algorithm
introduced for i.i.d. graphical model learning [13]. Our main technical
contribution is in providing sample complexity guarantee for our Markovian
dynamics. We prove that the algorithm can recover the true graph with high
probability using observations over $O(\log|V|)$ time steps.

The algorithm uses estimates of directed conditional entropy,
$\hat{H}(v_{+}|v,Q)$ for $Q\subset V$. For each node $v$, the algorithm starts
with empty neighborhood estimates $\hat{U}(v)$ and $\hat{T}(v)$; at every
step, it chooses the node that produces the maximum reduction in the directed
conditional entropy, i.e., the node for which
$\hat{H}(v_{+}|v,\hat{U}(v))-\hat{H}(v_{+}|v,\hat{U}(v),u)$ is maximum. If
this reduction is greater than a threshold, the corresponding node is added to
the neighborhood estimate ($\hat{U}(v)$). Nodes are added until no new node
results in a reduction greater than the threshold.

When the threshold is chosen as $\epsilon/2$, by Lemma 4, the last node added
to $\hat{U}(v)$ is always a true neighbor of $v$. Only the last node is added
to the final neighborhood estimate $\hat{T}(v)$. The whole process is
repeated, after initializing $\hat{U}(v)=\hat{T}(v)$. When $\hat{T}(v)$
becomes equal to the true neighborhood, no new node can be added to
$\hat{T}(v)$, and the algorithm stops at this point.

The sample complexity of this algorithm mainly depends on the dimension of the
probability distribution that needs to be estimated from the data. As the
algorithm starts with the smallest subsets and stops when there is no further
drop in entropy, the dimensionality of the estimated distributions is bounded.

To derive the algorithm’s sample complexity, we introduce a notation
$|P|_{max}$ to denote the maximum size of the subset $P$ from Lemma 4. We will
prove the upper bound on the value of $|P|_{max}$ in section VI. This means
that the error limit of $\epsilon/4$ is only required for estimating entropies
conditioned on subsets of size up to $|P|_{max}$. The support of random
variable $Y_{v}(t)$ is denoted by $\chi$ and
$|\xi|:=|\chi|^{(1+d(|P|_{max}+1))}$ is the maximum size of the joint
distribution that needs to be estimated. This, along with bounds on the
evolution of the Markov dynamics, results in logarithmic sample complexity,
which is formally stated in the Theorem 1.

######  Theorem 1.

The RecGreedy($\epsilon$) algorithm recovers the true influence graph with
probability at least $1-\gamma$ using data for

$\displaystyle T\geq
d+\frac{\left(\log(|V|)\cdot(|P|_{max}+1)+\log(2|\xi|/\gamma)\right)\cdot|\xi|^{2}}{(1-2(\bar{\mu}+L)|\rho(\tilde{A})|^{\frac{1}{d}})\delta^{2}}$

time steps where $|\xi|:=|\chi|^{(1+d(|P|_{max}+1))}$,
$|\chi|\leq(\frac{\bar{M}(\bar{M}+1)}{2}+2),|P|_{max}\leq\frac{2\log(\frac{\bar{M}(\bar{M}+1)}{2}+2)}{\epsilon}+1$,
$\delta=\frac{\epsilon^{2}}{8|\xi|}$, if
$2(\bar{\mu}+L)|\rho(\tilde{A})|^{\frac{1}{d}}<1$ and Assumption 1 holds. The
computational complexity is $\mathcal{O}(|V|)$.

Algorithm 1 RecGreedy($\epsilon$) Algorithm for $d$-Markov

Input: $\\{N_{v}(t):v\in V\\}_{t=0}^{T-1},\\{M_{v}(t):v\in V\\}_{t=0}^{T-1},$
Output: ${\hat{T}(v):v\in V}$

1:for $v\in V$ do

2: $\hat{T}(v),\leftarrow\phi$

3: repeat

4: $\hat{U}(v)\leftarrow\hat{T}(v)$

5: $\text{LastNode}\leftarrow\phi$

6: repeat

7: $\displaystyle u$ $\displaystyle=\arg\max\limits_{k\in
V\setminus\hat{U}(v)}\hat{H}(v_{+}|v,\hat{U}(v))$
$\displaystyle-\hat{H}(v_{+}|v,\hat{U}(v),k)$

8: $\Delta_{u}=\hat{H}(v_{+}|v,\hat{U}(v))-\hat{H}(v_{+}|v,\hat{U}(v),u)$

9: if $\Delta_{u}>\epsilon/2$ then

10: $\hat{U}(v)\leftarrow\hat{U}(v)\cup\\{u\\}$

11: $\text{LastNode}\leftarrow u$

12: else

13: $\hat{T}(v)\leftarrow\hat{T}(v)\cup\text{LastNode}$

14: end if

15: until New node is not added to $\hat{U}(v)$

16: until New node is not added to $\hat{T}(v)$

17:end for

##  VI Proofs

Our main results are stated in terms of the quantity $|P|_{max}$, which
denotes the maximum size of the subset of nodes the algorithm considers in any
of its iterations. This is because the algorithm is based on conditional
entropy estimation from the data; the accuracy of entropy estimation depends
on the dimensionality of the distribution under consideration. Therefore, we
structure our proof as follows: In Lemma 5, we state the value of $|P|_{max}$
for Algorithm 1. Then, in Lemma 6, we state the number of samples that are
sufficient to achieve the entropy estimation accuracy stated in Lemma 4.

After stating these results, we provide our proof for Theorem 1. The detailed
proofs of Lemmas 3\- 6 are provided in section VI-A. A novel contribution of
this work is the proof technique used for Lemma 6, where we use an improved
Markov Hoeffding inequality from [27] and obtain a bound on the second largest
eigenvalue of the concerned Markov chain. The bound on the second eigenvalue
was obtained by upper and lower bounding total variation distance with respect
to stationary distribution and then by comparing these bounds. This proof may
be of independent mathematical interest. We also give an indirect proof of the
non-degeneracy assumption in [13] for our particular $d$-Markov chain under
certain assumptions (Lemma 3).

######  Lemma 5.

The maximum possible size of the set $P\subset V$ for which the Algorithm 1
evaluates directed conditional entropy has an upper bound given by

$|P|_{max}\leq\frac{2\log(\frac{\bar{M}(\bar{M}+1)}{2}+2)}{\epsilon}+1$

######  Lemma 6.

If Assumption 1 holds and the number of time steps, $T\geq
d+\frac{\left(\log(|V|)\cdot(|P|_{max}+1)+\log(2|\xi|/\gamma)\right)\cdot|\xi|^{2}}{(1-2(\bar{\mu}+L)|\rho(\tilde{A})|^{\frac{1}{d}})\delta^{2}}$
where $\delta=\frac{\epsilon^{2}}{8|\xi|}$, then
$|\hat{H}(v_{+}|v,P)-H(v_{+}|v,P)|\leq\frac{\epsilon}{4}$ with probability at
least $1-\gamma$ .

Proof is given in section VI-A

###### Proof of Theorem 1:.

To recover the graph correctly, the RecGreedy($\epsilon$) algorithm requires
the following conditions: a) If the greedy step encounters a neighbor
$u\in\mathcal{N}_{v}$ of a node $v\in V$, it will be added. b) The last node
added by the greedy step is always a neighbor. From Lemma 4, we see that these
conditions are satisfied if the error in entropy estimation is limited to
$\epsilon/4$. From Lemma 6, we get the number of samples sufficient to limit
the entropy estimation error to $\epsilon/4$ with probability at least
$1-\gamma.$ Using Lemma 5, we substitute the value of
$|P|_{max}\leq\frac{2\log(\frac{\bar{M}(\bar{M}+1)}{2}+2)}{\epsilon}+1$ in
Lemma 6 and obtain the sample complexity stated in Theorem 1.

Proof of computational complexity: For every node, the greedy step runs at
most $|P|_{max}$ times. Therefore, the computational complexity is
$\mathcal{O}(|V|)$.

∎

###  VI-A Detailed Proofs

###### Proof of Lemma 3.

$H(v_{+}|v,Q)-H(v_{+}|v,Q,u)$ can be written as

$\displaystyle H(v_{+}|v,Q)-H(v_{+}|v,Q,u)$
$\displaystyle=\mathbb{E}_{v,Q,u,v_{+}}\left[\log\left(\frac{p(Y_{v_{+}}|\mathbf{Y}^{(d)}_{v,Q,u})}{p(Y_{v_{+}}|\mathbf{Y}^{(d)}_{v,Q})}\right)\right]$

Thus we see that $H(v_{+}|v,Q)-H(v_{+}|v,Q,u)$ is actually the expected K-L
divergence between two distributions over $Y_{v_{+}}$, where expectation is
taken with respect to $\mathbf{Y}^{(d)}_{v}$, $\mathbf{Y}^{(d)}_{Q}$ and
$\mathbf{Y}^{(d)}_{u}$, i.e.,

$\displaystyle H(v_{+}|v,Q)-H(v_{+}|v,Q,u)$
$\displaystyle=\mathbb{E}_{u,v,Q}\left[D_{y_{v_{+}}}\left({p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q,u})}||{p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q})}\right)\right]$

Using Pinsker’s inequality [31], for all $\mathbf{y}^{(d)}_{v,Q,u}$,

$\displaystyle\quad
D_{y_{v_{+}}}\left({p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q,u})}||{p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q})}\right)$
$\displaystyle\geq
2\|p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q,u})-p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q})\|_{1}^{2}$

Therefore, we have

$\displaystyle H(v_{+}|v,Q)-H(v_{+}|v,Q,u)$
$\displaystyle=\mathbb{E}_{v,Q}\mathbb{E}_{u}[D_{y_{v_{+}}}\left({p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q,u})}||{p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q})}\right)|v,Q]$
$\displaystyle\geq\mathbb{E}_{v,Q}\mathbb{E}_{u}[\|p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q,u})-p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q})\|_{1}^{2})|v,Q].$
$\displaystyle\geq\mathbb{E}_{v,Q}\big{[}\epsilon^{\prime}\mathbf{P}_{u}(\|p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q,u})-p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q})\|_{1}^{2}>\epsilon^{\prime})$
$\displaystyle\quad\mid v,Q]\big{]}$ (8)

By Assumption 1, there exists a $y_{u}$ such that
$\|p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q,u})-p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q})\|_{1}^{2})>\epsilon^{\prime}.$
For $M_{u}(t)=M$, we prove that the probability of such a $y_{u}$ is lower
bounded i.e., we prove that
$\inf_{y_{u},y_{v},y_{Q}}\mathbf{P}(Y_{u}=y_{u}|y_{v,Q})$ is lower bounded,
which would prove a lower bound on $H(v_{+}|v,Q)-H(v_{+}|v,Q,u)$.

$\mathbf{P}(Y_{u}=y_{u}|y_{v,Q})$ is a binomial distribution with parameter
$M$ and $\mathbb{E}[X_{u}|y_{v,Q}]$. The minimum probability of any binomial
with parameter $M$ and $a$ is $\min\\{a^{M},(1-a)^{M}\\}$. Now,

$\displaystyle X_{u}(t)=(1-\alpha_{u})((1-\beta)Z_{u}(t)+\beta l_{v})$
$\displaystyle+\alpha_{v}\sum_{u\in\mathcal{N}_{v}\bigcup{v}}\sum_{r=0}^{(d-1)}a_{uv}^{(r)}\\{\frac{N_{v}(t-1)}{M_{v}(t-1)}\\}$

and hence,
$a=\inf_{y_{v,Q}}\mathbb{E}[X_{u}|y_{v,Q}]=(1-\alpha_{u})((1-\beta)z_{u}+\beta
l_{v})+\alpha_{u}\theta_{u}$, where $0\leq\theta_{u}\leq 1$. Note that since
$0<z_{u}<1$, $0<\alpha<1$ we have $0<a<1$ and hence,
$\inf_{y_{u},y_{v},y_{Q}}\mathbf{P}(Y_{u}=y_{u}|y_{v,Q})=\min\\{a^{M},(1-a)^{M}\\}$.

Thus under Assumption 1,
$H(v_{+}|v,Q)-H(v_{+}|v,Q,u)>\epsilon:=\epsilon^{\prime}\cdot\min\\{a^{M},(1-a)^{M}\\}$
when $u$ is a neighbor of $v$.

Also by Lemma 1, if $Q$ contains $\mathcal{N}_{v}$ and
$u\notin\mathcal{N}_{v}$,
$p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q,u})=p(y_{v_{+}}|\mathbf{y}^{(d)}_{v,Q})$
for all $(y_{v_{+}},\mathbf{y}^{(d)}_{v,Q,u})$. Thus
$H(v_{+}|v,Q)-H(v_{+}|v,Q,u)=0$ when $Q$ contains the neighborhood.

∎

###### Proof of Lemma 4.

The proof follows on similar lines as [13]. If
$|\hat{H}(v_{+}|v,Q)-H(v_{+}|v,Q)|<\epsilon^{\prime\prime}$ for $v\in V\text{
and }Q\subseteq V\text{ such that }|Q|<|P|_{max}$, then

$H(v_{+}|v,Q)-\epsilon^{\prime\prime}<\hat{H}(v_{+}|v,Q)<\epsilon^{\prime\prime}+H(v_{+}|v,Q)$

Similar condition holds for $\hat{H}(v_{+}|v,Q,u),u\in V.$ When we combine
these two conditions and apply Lemma 3, we conclude that
$\hat{H}(v_{+}|v,Q)-\hat{H}(v_{+}|v,Q,u)$ is greater than
$\epsilon-2\epsilon^{\prime\prime}$ if $u$ is a neighbor and is less than
$2\epsilon^{\prime\prime}$ if $\mathcal{N}_{v}\subseteq Q\text{ and
}u\notin\mathcal{N}_{v}$.

Substituting $\epsilon^{\prime\prime}=\epsilon/4$ proves Lemma 4. ∎

###### Proof of Lemma 5.

In the context of this work, the algorithm starts with an initial entropy of
$\hat{H}(v_{+}|v)$, which is reduced by $\epsilon/2$ every time a neighbor is
added. The maximum number of nodes that can be estimated as neighbors is
(similar to [13]) $\frac{2\hat{H}(v_{+}|v)}{\epsilon}$. Therefore, the maximum
number of steps $|P|_{max}$ is $\frac{2\hat{H}(v_{+}|v)}{\epsilon}+1$. We know
that $\hat{H}(v_{+}|v)\leq\log|\chi|$ where $|\chi|$ is the support of
$y_{v_{+}}$. Here, $|\chi|\leq(\frac{\bar{M}(\bar{M}+1)}{2}+2)$.Thus,
$|P|_{max}\leq\frac{2\log(\frac{\bar{M}(\bar{M}+1)}{2}+2)}{\epsilon}+1$. ∎

Next, in order to prove Lemma 6, we first state a few results in Lemmas 7-10.

######  Lemma 7.

$|\hat{H}(v_{+}|v,P)-H(v_{+}|v,P)|<\epsilon/4$ if
$\|\hat{p}(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})-p(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})\|_{1}<\delta$,
where $\hat{p}(\cdot)$ is the estimated probability distribution, $p(\cdot)$
is the true probability distribution and $\delta=\frac{\epsilon^{2}}{8|\xi|}$

###### Proof.

The entropy estimation error depends on the error in the estimation of
probability as [26]

$\begin{split}|\hat{H}(v_{+},v,P)-H(v_{+},v,P)|\leq-\|\hat{p}(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})-\\\
p(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})\|_{1}\times\log\frac{\|\hat{p}(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})-p(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})\|_{1}}{|\xi|}\end{split}$
(9)

If
$\|\hat{p}(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})-p(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})\|_{1}\leq\delta$,
then $|\hat{H}(v_{+},v,P)-H(v_{+},v,P)|\leq\sqrt{\delta|\xi|}$. Since the
conditional entropy is computed from the estimated joint entropy as
$\hat{H}(v_{+}|v,P)=\hat{H}(v_{+},v,P)-\hat{H}(v,P)$, we have
$|\hat{H}(v_{+}|v,P)-H(v_{+}|v,P)|\leq 2\sqrt{\delta|\xi|}=\epsilon/4$. ∎

Next, let $\lambda^{*}$ denote the second largest eigenvalue of the Markov
chain $\\{\mathbf{Y}^{(d)}(t)\\}$. Lemma 8 states the sample complexity in
terms of $\lambda^{*}$ using the Markov Hoeffding inequality from [27]. Later,
in Lemmas 9 and 10, we upper bound the absolute value of $\lambda^{*}$.

######  Lemma 8.

If $|\lambda|^{*}<1$ and the number of data points $T$ available for
probability estimation is greater than
$d+\frac{\left(\log(|V|)\cdot(|P|_{max}+1)+\log(2|\xi|/\gamma)\right)\cdot|\xi|^{2}}{(1-|\lambda^{*}|)\delta^{2}}$,
then
$\|\hat{p}(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})-p(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})\|_{1}<\delta$
with probability at least $1-\gamma$.

###### Proof.

Consider the Markov chain
$\mathbf{Z}(t)=(\mathbf{Y}^{(d)}(t),\mathbf{Y}^{(d)}(t-1))$. We define a
function
$f(\mathbf{Z}(t))=\frac{\sum_{t=d}^{T-1}\mathbb{I}\left((Y_{v}(t),\mathbf{Y}^{(d)}_{v,P}(t-1))=(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})\right)}{T-d},$
where $\mathbb{I}$ is the indicator function. Using the concentration
inequality stated in [27] for $\mathbf{Z}(t)$ and $f$, we get

$\displaystyle~{}~{}~{}\mathbf{P}(|f(\mathbf{Z}(t))-\mathbb{E}f(\mathbf{Z}(t))|>\delta_{1})$
$\displaystyle\leq
2\exp\left({-\frac{(1-|\lambda^{*}|)2(T-d)^{2}\delta_{1}^{2}}{(1+|\lambda^{*}|)(T-d)}}\right)$

By applying union bound on $(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})$, we have

$\displaystyle~{}~{}\mathbf{P}(\|\hat{p}(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})-p(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})\|_{1}>\delta)$
$\displaystyle\leq
2|\xi|\exp\left({-\frac{(1-|\lambda^{*}|)2(T-d)(\delta/|\xi|)^{2}}{(1+|\lambda^{*}|)}}\right),$

which implies

$\displaystyle\mathbf{P}(\|\hat{p}(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})-p(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})\|_{1}>\delta)\leq$
$\displaystyle
2|\xi|\exp\left({-\frac{2(1-|\lambda^{*}|)(T-d)\delta^{2}}{(1+|\lambda^{*}|).(|\xi|^{2})}}\right)$

Applying union bound on $v$ and $P$,

$\displaystyle\mathbf{P}(\|\hat{p}(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})-p(y_{v_{+}},\mathbf{y}^{(d)}_{v,P})\|_{1}>\delta)\leq$
$\displaystyle
2|V|{{|V|}\choose{|P|_{max}}}|\xi|\exp\left({-\frac{2(1-|\lambda^{*}|)(T-d)\delta^{2}}{(1+|\lambda^{*}|)(|\xi|^{2})}}\right)$

This probability is less than or equal to $\gamma$ when

$\displaystyle\gamma\geq
2|V|\cdot{{|V|}\choose{|P|_{max}}}|\xi|\exp\left({-\frac{2(1-|\lambda^{*}|)(T-d)\delta^{2}}{(1+|\lambda^{*}|)(|\xi|^{2})}}\right)$
$\displaystyle\implies$ $\displaystyle\gamma\geq 2|V|\cdot|V|^{|P|_{max}}$
$\displaystyle~{}~{}~{}\times|\xi|\exp\left({-\frac{2(1-|\lambda^{*}|)(T-d)\delta^{2}}{(1+|\lambda^{*}|)|\xi|^{2}}}\right)$

Therefore,

$\displaystyle T\geq
d+\frac{(|\xi|^{2})(\log(|V|)(|P|_{max}+1)+\log(\frac{2|\xi|}{\gamma}))}{(1-|\lambda^{*}|)\delta^{2}}$

∎

######  Lemma 9.

$2c|\lambda^{*}|^{t}\leq\frac{1}{2}\max\limits_{\mathbf{y}^{(d)}_{0}}\|\mathcal{P}^{t}(\mathbf{y}^{(d)}_{0},.)-\bm{\pi}\|_{1}$,
where $\mathcal{P}$ is the transition matrix,
$\mathcal{P}^{t}(\mathbf{y}^{(d)}_{0},.)=\mathbf{P}(\mathbf{Y}^{(d)}(t)|\mathbf{Y}^{(d)}(0)=\mathbf{y}^{(d)}_{0})$,
$\bm{\pi}$ is the stationary distribution and $c$ is a constant independent of
$t$ and $|V|$.

###### Proof.

Consider the total variation distance
$\mbox{dist}(t):=\frac{1}{2}\max\limits_{\mathbf{y}^{(d)}_{0}}\|\mathcal{P}^{t}(\mathbf{y}^{(d)}_{0},.)-\bm{\pi}\|_{1}$.

$\displaystyle\mathcal{P}^{t}(\mathbf{y}^{(d)}_{0},.)-\bm{\pi}$
$\displaystyle=\mathbf{e}_{0}\left(\mathbbm{1}^{T}\bm{\pi}+\sum\limits_{i=2}^{N}\lambda_{i}^{t}\mathbf{w}_{i}\mathbf{z}_{i}^{T}\right)-\bm{\pi}$
$\displaystyle=\mathbf{e}_{0}\left(\sum\limits_{i=2}^{N}\lambda_{i}^{t}\mathbf{w}_{i}\mathbf{z}_{i}^{T}\right),$

where $\mathbf{e}_{0}$ is the initial state vector given the initial state is
$\mathbf{y}^{(d)}_{0}$, $\mathbf{w}_{i}$ and $\mathbf{z}_{i}$ are the right
and left eigenvectors of $\mathcal{P}$ and $N$ is the number of states. Now,

$\displaystyle\|\mathcal{P}^{t}(\mathbf{y}^{(d)}_{0},.)-\bm{\pi}\|_{1}=\left\|\mathbf{e}_{0}\left(\sum\limits_{i=2}^{N}\lambda_{i}^{t}\mathbf{w}_{i}\mathbf{z}_{i}^{T}\right)\right\|_{1}$
$\displaystyle=|\lambda^{*}|^{t}\left\|\mathbf{e}_{0}\left(\mathbf{w}^{*}\mathbf{z}^{*T}+\sum_{\begin{subarray}{c}i=3\end{subarray}}^{N}\left(\frac{\lambda_{i}}{\lambda^{*}}\right)^{t}\mathbf{w}_{i}\mathbf{z}_{i}^{T}\right)\right\|_{1}$

For all sufficiently large $t$,
$\left\|\mathbf{e}_{0}\sum_{\begin{subarray}{c}i=1\\\ i\neq
i^{*}\end{subarray}}^{N}\left(\frac{\lambda_{i}}{\lambda^{*}}\right)^{t}\mathbf{w}_{i}\mathbf{z}_{i}^{T}\right\|_{1}$
is exponentially smaller than
$\left\|\mathbf{e}_{0}\mathbf{w}^{*}\mathbf{z}^{*T}\right\|$ since
$\lambda_{i}<\lambda^{*}$ for each $i\geq 3$. Since
$\left\|\mathbf{e}_{0}\left(\mathbf{w}^{*}\mathbf{z}^{*T}+\sum_{\begin{subarray}{c}i=3\end{subarray}}^{N}\left(\frac{\lambda_{i}}{\lambda^{*}}\right)^{t}\mathbf{w}_{i}\mathbf{z}_{i}^{T}\right)\right\|_{1}$
is lower bounded by
$\left\|\mathbf{e}_{0}\mathbf{w}^{*}\mathbf{z}^{*T}\right\|_{1}-\left\|\mathbf{e}_{0}\sum_{\begin{subarray}{c}i=1\\\
i\neq
i^{*}\end{subarray}}^{N}\left(\frac{\lambda_{i}}{\lambda^{*}}\right)^{t}\mathbf{w}_{i}\mathbf{z}_{i}^{T}\right\|_{1}$,
we get that for all sufficiently large $t$,
$\left\|\mathbf{e}_{0}\left(\mathbf{w}^{*}\mathbf{z}^{*T}+\sum_{\begin{subarray}{c}i=3\end{subarray}}^{N}\left(\frac{\lambda_{i}}{\lambda^{*}}\right)^{t}\mathbf{w}_{i}\mathbf{z}_{i}^{T}\right)\right\|_{1}\geq\frac{1}{2}\left\|\mathbf{e}_{0}\mathbf{w}^{*}\mathbf{z}^{*T}\right\|_{1}=c$.
Thus, for sufficiently large $t_{1}$ we get

$\mbox{dist}(t)\geq 2c|\lambda^{*}|^{t}\quad\text{for all }t>t_{1}$ (10)

∎

######  Lemma 10.

$\mbox{dist}(t)\leq
c^{\prime}\left(2(\bar{\mu}+L)|\rho(\tilde{A})|^{\frac{1}{d}}\right)^{t}$ for
some constant $c^{\prime}$.

###### Proof.

We know that

$\displaystyle\mbox{dist}(t)\leq\bar{\mbox{dist}(t)}:=\frac{1}{2}\max\limits_{\mathbf{y}^{(d)}_{0},\mathbf{y}^{\prime(d)}_{0}}\Big{\|}\mathbf{P}(\mathbf{Y}^{(d)}(t)|\mathbf{Y}^{(d)}(0)=\mathbf{y}^{(d)}_{0})$
$\displaystyle-\mathbf{P}(\mathbf{Y}^{(d)}(t)|\mathbf{Y}^{(d)}(0)=\mathbf{y}^{\prime(d)}_{0})\Big{\|}_{1}$

Now we consider a coupling of the Markov chain, the process
$\\{(\mathbf{Y}^{(d)}(t),\mathbf{Y}^{\prime(d)}(t))\\}$ with the initial
conditions
$\mathbf{Y}^{(d)}(0)=\mathbf{y}^{(d)}(0),\mathbf{Y}^{\prime(d)}(0)=\mathbf{y}^{\prime(d)}(0)$.
We generate these two processes as follows:

At time $t$, we only need to generate $\mathbf{Y}(t)$ and
$\mathbf{Y}^{\prime}(t)$; the rest of $d-1$ terms are taken from
$\mathbf{Y}^{(d)}(t-1)$ and $\mathbf{Y}^{\prime(d)}(t-1)$. To generate
$\mathbf{Y}(t)$ and $\mathbf{Y}^{\prime}(t)$, we first generate
$\mathbf{X}(t)$ and $\mathbf{X}^{\prime}(t)$ from $\mathbf{Y}^{(d)}(t-1)$ and
$\mathbf{Y}^{\prime(d)}(t-1)$ according to eq(III).

Now, we denote introduce a random variable $R_{v}(t)$ to denote
$\text{Poisson}(\min\\{\mu_{v}(X_{v}(t)),\mu_{v}(X^{\prime}_{v}(t))\\})$ and
$S_{v}(t)$ to denote
$\text{Poisson}(\max\\{\mu_{v}(X_{v}(t)),\mu_{v}(X^{\prime}_{v}(t))\\}-\min\\{\mu_{v}(X_{v}(t)),\mu_{v}(X^{\prime}_{v}(t))\\}).$
If the minimum among $\mu_{v}(X_{v}(t))$ and $\mu_{v}(X^{\prime}_{v}(t))\\}$
is $\mu_{v}(X_{v}(t))$, then generate $M_{v}(t)=\min\\{R_{v}(t),\bar{M}\\}+1$
and $M^{\prime}_{v}(t)=\min\\{R_{v}(t)+S_{v}(t),\bar{M}\\}+1$.

Else, if If the minimum among $\mu_{v}(X_{v}(t))$ and
$\mu_{v}(X^{\prime}_{v}(t))\\}$ is $\mu_{v}(X^{\prime}_{v}(t))$, then generate
$M^{\prime}_{v}(t)=\min\\{R_{v}(t),\bar{M}\\}+1$ and
$M_{v}(t)=\min\\{R_{v}(t)+S_{v}(t),\bar{M}\\}+1$.

To generate $\mathbf{Y}^{(d)},\mathbf{Y}^{\prime(d)}$, we generate
$\min\\{M_{v}(t),M^{\prime}_{v}(t)\\}$ number of $\text{Unif}[0,1]$ random
variables (the $i^{th}$ uniform random variable is denoted by $U_{i}$). Then
compare each of these uniform random variables with $X_{v}(t)$ to generate
observed binary variables for the process $\mathbf{Y}^{(d)}(t)$ (if $X_{v}(t)$
is less than or equal to the uniform variable, it’s a ‘1’). Then, compare the
same uniform random variables with $X^{\prime}_{v}(t)$ to generate the
observed binary variables for $\mathbf{Y}^{\prime(d)}(t).$ Here,
$\min\\{M_{v}(t),M^{\prime}_{v}(t)\\}$ uniform random variables are common to
both processes. Then, generate the remaining binary variables using a new set
of $\max\\{M_{v}(t),M^{\prime}_{v}(t)\\}-\min\\{M_{v}(t),M^{\prime}_{v}(t)\\}$
uniform random variables.

Now we will consider the distance $\bar{d}(t)$ between these two processes.

$\bar{\mbox{dist}(t)}\leq\mathbf{P}(\tau_{couple}>t|\mathbf{Y}^{(d)}(0)=\mathbf{y}^{(d)}_{0},\mathbf{Y}^{\prime(d)}(0)=\mathbf{y}^{\prime(d)}_{0}),$

For ease of notation, from here on, we denote the initial condition event
$\\{\mathbf{Y}^{(d)}(0)=\mathbf{y}^{(d)}_{0},\mathbf{Y}^{\prime(d)}(0)=\mathbf{y}^{\prime(d)}_{0}\\}$
by $\mathbb{\tau_{0}}$ where
$\tau_{couple}:=\min\\{t:\mathbf{Y}^{(d)}(t)=\mathbf{Y}^{\prime(d)}(t)\\}$.

$\displaystyle\mathbf{P}(\tau_{couple}>t\mid\mathbb{\tau_{0}})=\mathbf{P}(\mathbf{Y}^{(d)}(t)\neq\mathbf{Y}^{\prime(d)}(t)\mid\mathbb{\tau_{0}})$

$\mathbf{Y}^{(d)}(t)$ and $\mathbf{Y}^{\prime(d)}(t)$ are not equal when
$\mathbf{Y}_{v}(k)\neq\mathbf{Y}^{\prime}_{v}(k)$ for at least one $(k,v)$
such that $k\in\\{t,t-1,\ldots,t-d+1\\},\text{ and }v\in\\{1,2,\ldots|V|\\}$.
Therefore, we apply the union bound on $k$ and $v$ to obtain eq(11)

$\displaystyle\mathbf{P}(\tau_{couple}>t\mid\mathbb{\tau_{0}})$
$\displaystyle=\mathbf{P}\Big{(}\bigcup\limits_{r=0}^{d-1}\\{\mathbf{Y}(t-r)\neq\mathbf{Y}^{\prime}(t-r)\\}\mid\mathbb{\tau_{0}}\Big{)}$
$\displaystyle\leq\sum\limits_{r=0}^{d-1}\mathbf{P}(\mathbf{Y}(t-r)\neq\mathbf{Y}^{\prime}(t-r)\mid\mathbb{\tau_{0}})$
$\displaystyle=\sum\limits_{k=t-d+1}^{t}\mathbf{P}(\mathbf{Y}(k)\neq\mathbf{Y}^{\prime}(k)\mid\mathbb{\tau_{0}})$
$\displaystyle\leq\sum\limits_{v=1}^{|V|}\sum\limits_{k=t-d+1}^{t}\mathbf{P}(Y_{v}(k)\neq
Y^{\prime}_{v}(k)\mid\mathbb{\tau_{0}})$ (11)

The event $Y_{v}(k)\neq Y^{\prime}_{v}(k)$ can be partitioned into events
based on whether the total number of binary random variables is equal for the
two processes. Therefore, $\mathbf{P}(Y_{v}(k)\neq
Y^{\prime}_{v}(k)\mid\mathbb{\tau_{0}})$ can be written as

$\displaystyle\mathbf{P}(Y_{v}(k)\neq Y^{\prime}_{v}(k)\mid\mathbb{\tau_{0}})$
$\displaystyle=\mathbf{P}(M_{v}(k)=M^{\prime}_{v}(k),Y_{v}(k)\neq
Y^{\prime}_{v}(k)|\mathbb{\tau_{0}})$ $\displaystyle+\mathbf{P}(M_{v}(k)\neq
M^{\prime}_{v}(k),Y_{v}(k)\neq Y^{\prime}_{v}(k)|\mathbb{\tau_{0}})$ (12)

Now, we consider the first term of eq (VI-A). The joint probability
$\mathbf{P}(M_{v}(k)=M^{\prime}_{v}(k),Y_{v}(k)\neq Y^{\prime}_{v}(k)|\cdot)$
is the same as $\mathbf{P}(M_{v}(k)=M^{\prime}_{v}(k),N_{v}(k)\neq
N^{\prime}_{v}(k)|\cdot)$, since $Y_{v}(k)=N_{v}(k)/M_{v}(k)$. Then, we apply
the law of total probability with respect to the smallest disjoint events
generated by $\\{\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1)\\}$ to
obtain
$\mathbb{E}_{\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1)}[\mathbf{P}(N_{v}(k)\neq
N^{\prime}_{v}(k),M_{v}(k)=M^{\prime}_{v}(k)\mid\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1);\mathbb{\tau_{0}})].$
This is an expectation over the distribution of
$\\{\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1)\\}$ given the initial
condition ${}^{\prime}\mathbb{\tau_{0}}^{\prime}$.

Also, any joint probability is less than or equal to conditional probability.
Therefore, the first term of eq (VI-A) turns out to be less than or equal to
$\mathbb{E}_{\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1)}[\mathbf{P}(N_{v}(k)\neq
N^{\prime}_{v}(k)\mid
M_{v}(k)=M^{\prime}_{v}(k);\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1);\mathbb{\tau_{0}})].$

$\displaystyle\mathbf{P}\big{(}M_{v}(k)=M^{\prime}_{v}(k),Y_{v}(k)\neq
Y^{\prime}_{v}(k)|\mathbb{\tau_{0}}\big{)}$
$\displaystyle\leq\mathbb{E}_{\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1)}[\mathbf{P}(N_{v}(k)\neq
N^{\prime}_{v}(k)$ $\displaystyle\quad\mid
M_{v}(k)=M^{\prime}_{v}(k);\mathbb{\tau_{0}};\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1))]$
(13)

Now, we analyze the term $\mathbf{P}(N_{v}(k)\neq N^{\prime}_{v}(k)\mid
M_{v}(k)=M^{\prime}_{v}(k);\mathbb{\tau_{0}};\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1)).$
Since $M_{v}(k)=M^{\prime}_{v}(k)$, all the uniform random variables are
common to both processes. $N_{v}(k)$ and $N^{\prime}_{v}(k)$ will not be equal
when $X_{v}(k)\leq U_{i}<X^{\prime}_{v}(k)$ or $X^{\prime}_{v}(k)\leq
U_{i}<X_{v}(k)$ for at least one out of the $M_{v}(k)$ uniform random
variables $\\{U_{i}\\}_{i=1}^{M_{v}(k)}$. Therefore, we can apply union bound
on the uniform random variables and use the law of total probability over
$M_{v}(k).$ Further, the probability that a uniform random variable lies in
$[X_{v}(k),X^{\prime}_{v}(k))\cup[X^{\prime}_{v}(k),X_{v}(k))$ is given by
$|X_{v}(k)-X^{\prime}_{v}(k)|$. Finally, using the model eq (III), we can
write this expression in terms of
$\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1).$

$\displaystyle\quad\mathbf{P}(N_{v}(k)\neq N^{\prime}_{v}(k)\mid$
$\displaystyle\quad
M_{v}(k)=M^{\prime}_{v}(k);\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1);\mathbb{\tau_{0}})$
$\displaystyle\leq\mathbb{E}_{M_{v}(k)}\Big{[}M_{v}(k)\mathbf{P}_{U}(\\{U\in[X_{v}(k),X^{\prime}_{v}(k))\cup[X^{\prime}_{v}(k),X_{v}(k))$
$\displaystyle\quad\mid
M_{v}(k)=M^{\prime}_{v}(k);\mathbb{\tau_{0}};\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1))\Big{]}$
$\displaystyle\leq\max\limits_{X_{v}(k)}\mathbb{E}_{M_{v}(k)}[M_{v}(k)]\lvert
X_{v}(k)-X^{\prime}_{v}(k)\rvert$
$\displaystyle=\bar{\mu}\left[\alpha_{v}\left|\sum_{u\in\mathcal{N}_{v}\cup\\{v\\}}\mathbf{a}^{(d)}_{uv}\mathbf{Y}^{(d)}(k-1)-\mathbf{a}^{(d)}_{uv}\mathbf{Y}^{\prime(d)}(k-1)\right|\right]$

Thus, we have

$\displaystyle\mathbf{P}\big{(}M_{v}(k)=M^{\prime}_{v}(k),Y_{v}(k)\neq
Y^{\prime}_{v}(k)|\mathbb{\tau_{0}}\big{)}$
$\displaystyle\leq\bar{\mu}\mathbb{E}_{\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1)}\Bigg{[}\alpha_{v}\Big{|}\sum_{u\in\mathcal{N}_{v}\cup\\{v\\}}\mathbf{a}^{(d)}_{uv}\mathbf{Y}^{(d)}(k-1)$
$\displaystyle-\mathbf{a}^{(d)}_{uv}\mathbf{Y}^{\prime(d)}(k-1)\Big{|}\mid\mathbb{\tau_{0}}\Bigg{]}$
(14)

Now we derive an expression for the second term of eq (VI-A), i.e.,
$\mathbf{P}(M_{v}(k)\neq M^{\prime}_{v}(k),Y_{v}(k)\neq
Y^{\prime}_{v}(k)|\mathbb{\tau_{0}})$ in terms of $\mathbf{Y}^{(d)}(k-1)$ and
$\mathbf{Y}^{\prime(d)}(k-1)$.

$\displaystyle\mathbf{P}(M_{v}(k)\neq M^{\prime}_{v}(k),Y_{v}(k)\neq
Y^{\prime}_{v}(k)|\mathbb{\tau_{0}})$
$\displaystyle\leq\mathbf{P}(M_{v}(k)\neq M^{\prime}_{v}(k),\mid Y_{v}(k)\neq
Y^{\prime}_{v}(k);$
$\displaystyle\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1);\mathbb{\tau_{0}})$
$\displaystyle=\mathbb{E}_{\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1)}[\mathbf{P}(M_{v}(k)\neq
M^{\prime}_{v}(k)\mid$ $\displaystyle Y_{v}(k)\neq
Y^{\prime}_{v}(k);\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1);\mathbb{\tau_{0}})]$

Since we have used two Poisson random variables
$R_{v}(t)=\text{Poisson}(\max\\{\mu_{v}(X_{v}(k),\mu_{v}(X^{\prime}_{v}(k)\\})$
and
$S_{v}(t)=\text{Poisson}(\max\\{\mu_{v}(X_{v}(k),\mu_{v}(X^{\prime}_{v}(k)\\}-\min\\{\mu_{v}(X_{v}(k),\mu_{v}(X^{\prime}_{v}(k)\\})$,
we have

$\displaystyle~{}~{}~{}\mathbf{P}(M_{v}(k)\neq M^{\prime}_{v}(k)\mid
Y_{v}(k)\neq Y^{\prime}_{v}(k);$
$\displaystyle~{}~{}~{}~{}~{}~{}\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1);\mathbb{\tau_{0}})$
$\displaystyle=\mathbf{P}(S_{v}(t)\neq 0\mid Y_{v}(k)\neq Y^{\prime}_{v}(k);$
$\displaystyle~{}~{}~{}\quad\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1);\mathbb{\tau_{0}})$
$\displaystyle\leq\mathbb{E}[1-\exp{(-|\mu_{v}(X_{v}(k))-\mu_{v}(X^{\prime}_{v}(k))|)}|$
$\displaystyle~{}~{}~{}~{}\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1),\mathbb{\tau_{0}}]$
$\displaystyle\leq\mathbb{E}[|\mu_{v}(X_{v}(k))-\mu_{v}(X^{\prime}_{v}(k))|]|$
$\displaystyle~{}~{}~{}~{}\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1),\mathbb{\tau_{0}}]$
$\displaystyle\leq
L\mathbb{E}[|X_{v}(k)-X^{\prime}_{v}(k)|\mid\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1)\text{
and }\mathbb{\tau_{0}}]$
$\displaystyle=L\left|\alpha_{v}\left(\sum_{u\in\mathcal{N}_{v}\cup\\{v\\}}\mathbf{a}^{(d)}_{uv}\mathbf{Y}^{(d)}(k-1)-\mathbf{a}^{(d)}_{uv}\mathbf{Y}^{\prime(d)}(k-1)\right)\right|$

Therefore,

$\displaystyle\mathbf{P}(Y_{v}(k)\neq Y^{\prime}_{v}(k)\mid\mathbb{\tau_{0}})$
$\displaystyle\leq(\bar{\mu}+L)\mathbb{E}_{\mathbf{Y}^{(d)}(k-1),\mathbf{Y}^{\prime(d)}(k-1)}\Bigg{[}\alpha_{v}\Big{|}\sum_{u\in\mathcal{N}_{v}\cup\\{v\\}}\mathbf{a}^{(d)}_{uv}\mathbf{Y}^{(d)}(k-1)$
$\displaystyle\quad-\mathbf{a}^{(d)}_{uv}\mathbf{Y}^{\prime(d)}(k-1)\Big{\rvert}\mid\mathbb{\tau_{0}}\Bigg{]}$

Since $0\leq Y_{v}(k),Y^{\prime}_{v}(k)\leq
1,\mathbb{E}[|Y_{v}(k)-Y^{\prime}_{v}(k)|]\leq\mathbf{P}(Y_{v}(k)\neq
Y^{\prime}_{v}(k))$. Therefore,

$\displaystyle\mathbf{P}(Y_{v}(k)\neq Y^{\prime}_{v}(k)\mid\mathbb{\tau_{0}})$
$\displaystyle=(\bar{\mu}+L)\sum\limits_{u}\sum\limits_{r=0}^{(d-1)}a_{uv}^{(r)}\mathbf{P}(Y_{u}(k-1-r)\neq
Y^{\prime}_{u}(k-1-r)$ $\displaystyle\mid\mathbb{\tau_{0}})$ (15)

Now, let $\eta_{v}(k):=\mathbf{P}(Y_{v}(k)\neq
Y^{\prime}_{v}(k)\mid\mathbb{\tau_{0}})$. Also, define a vector
$\eta(k):=[\eta_{1}(k),\eta_{2}(k),\ldots,\eta_{|V|}(k)]^{T}$

From eq (VI-A) we have the following element wise inequality for $\eta(k)$:

$\displaystyle\eta(k)$
$\displaystyle\leq(\bar{\mu}+L)\sum\limits_{r=0}^{d-1}\tilde{A}(r)\eta(k-1-r)$
(16)

We first choose an $r$ for which $\rho(\tilde{A}(r))$ is the maximum and
denote it by $r^{*}$. Next, we expand eq (16) until we reach the initial
condition, i.e. $\eta(0)=\eta(-1)=\ldots\eta(-d+1)=\mathbbm{1}.$ This
expansion gives an inequality with a sum of powers of $\tilde{A}$ multiplied
with the initial probability vector $\mathbbm{1}.$ Let the number of such
terms at time $k$ be $\kappa(k).$ We have
$\kappa(k)=\kappa(k-1)+\kappa(k-2)+\ldots+\kappa(k-d)\text{ for }k>d$ and
$\kappa(k)=\sum\limits_{l=0}^{k-2}2^{l}(d-(k-2-l))+(d-r+1),\text{ for
}k\in\\{1,2,\ldots d\\}.$ It follows after some calculations that any
$\kappa(k)$ satisfying the above recursion also satisfies $\kappa(k)\leq
2^{k+d}$. Also, the lowest power of $\tilde{A}(r)$ for any $r$ that appears is
$\tilde{A}(r)^{k/d}$ Therefore, we have

$\displaystyle\eta(k)$
$\displaystyle\leq(\bar{\mu}+L)^{k}\cdot{2}^{k+d}\tilde{A}(r^{*})^{k/d}\mathbbm{1}$

From the union bound in eq (11) we have

$\displaystyle\mbox{dist}(t)$ $\displaystyle\leq
d(\bar{\mu}+L)^{t}\cdot{2}^{t+d}\|\tilde{A}(r^{*})^{t/d}\mathbbm{1}\|_{1}$
(17)

Let us write $\tilde{A}(r^{*})$ in terms of its eigenvalues and eigenvectors
as
$\tilde{A}(r^{*})^{t}=\sum_{i}\lambda_{A_{i}}^{t}\mathbf{w}_{A_{i}}\mathbf{z}_{A_{i}}^{T}=\rho(\tilde{A})^{t}\left(\mathbf{w}_{A}\mathbf{z}_{A}^{T}+\sum_{i=2}^{|V|}(\frac{\lambda_{A_{i}}}{\rho(\tilde{A}})^{t}\mathbf{w}_{A_{i}}\mathbf{z}_{A_{i}}^{T}\right)$.
As $\lambda_{A_{i}}<\rho(\tilde{A})$,
$\left\|\sum_{i=2}^{|V|}(\frac{\lambda_{A_{i}}}{\rho(\tilde{A}})^{t}\mathbf{w}_{A_{i}}\mathbf{z}_{A_{i}}^{T})\mathbbm{1}\right\|$
is exponentially smaller than
$\|\mathbf{w}_{A_{i}}\mathbf{z}_{A_{i}}\mathbbm{1}\|$ Therefore,

$\displaystyle\mbox{dist}(t)\leq
d(\bar{\mu}+L)^{t}\cdot{2}^{t+d}\|\tilde{A}^{t/d}\mathbbm{1}\|_{1}$
$\displaystyle\leq d(\bar{\mu}+L)^{t}\cdot{2}^{t+d}|\rho(\tilde{A})|^{t/d}$
$\displaystyle\quad\times\left\|\left(\mathbf{w}_{A}\mathbf{z}_{A}^{T}+\sum_{i=2}^{|V|}(\frac{\lambda_{A_{i}}}{\rho(\tilde{A})})^{t/d}\mathbf{w}_{A_{i}}\mathbf{z}_{A_{i}}^{T}\right)\mathbbm{1}\right\|_{1}$
$\displaystyle\leq d(\bar{\mu}+L)^{t}\cdot{2}^{t+d}|\rho(\tilde{A})|^{t/d}$
$\displaystyle\quad\times\left(\|\mathbf{w}_{A}\mathbf{z}_{A}^{T}\mathbbm{1}\|_{1}+\left\|{}_{1}\sum_{i=2}^{|V|}(\frac{\lambda_{A_{i}}}{\rho(\tilde{A})})^{t/d}\mathbf{w}_{A_{i}}\mathbf{z}_{A_{i}}^{T}\mathbbm{1}\right\|_{1}\right)$
$\displaystyle\leq
d(\bar{\mu}+L)^{t}\cdot{2}^{t+d}|\rho(\tilde{A})|^{t/d}\cdot
2\|\mathbf{w}_{A}\mathbf{z}_{A}^{T}\mathbbm{1}\|_{1}$ for all $t>t_{2}$ where
$t_{2}$ is sufficiently large (18)

We get the final bound by choosing
$c^{\prime}=2^{d+1}d\|\mathbf{w}_{A}\mathbf{z}_{A}^{T}\mathbbm{1}\|_{1}$.

∎

###### Proof of Lemma 6.

Combining eq(10) and (VI-A), we have

$\displaystyle|\lambda|^{*}\leq(2c/c^{\prime})^{(1/t)}\left(2(\bar{\mu}+L)|\rho(\tilde{A})|^{\frac{1}{d}}\right)$
for all $t>\max\\{t_{1},t_{2}\\}$ $\displaystyle\therefore|\lambda|^{*}\leq
2(\bar{\mu}+L)|\rho(\tilde{A})|^{\frac{1}{d}}$

Next we substitute this upper bound on $|\lambda|^{*}$ in Lemma 8, along with
the appropriate value of $\delta$ (Lemma 7). This proves Lemma 6. ∎

##  VII Simulations

Though our main contribution is analytical in nature, we perform multiple
simulations to corroborate our analytical results. Interestingly, as is often
the case, the simulations also lead to certain insights, which are
analytically intractable.

(a)

(b)

(c)

Figure 1: (a) Probability of perfect recovery for $|V|=7$,$d=1$, $\bar{M}=0$,
(b)Probability of perfect recovery for $|V|=7,d=1$, $\bar{M}=2$ and (c)
Probability of perfect recovery for $|V|=7,d=2,\bar{M}=2$ Figure 2:
Probability of perfect recovery for $|V|=7,\bar{M}=2,d=2$ when
$2(\bar{\mu}+L)|\rho(\tilde{A})|^{\frac{1}{d}}>1$

We generate data using Eq. III for various values of $d$ and $\bar{M}$ and run
Algorithm 1 on that data. We chose the parameters of the model as follows: for
all nodes $v$, $\beta=0.75$, $l_{v}=0.167$, $\alpha_{v}=0.4$,
$Z_{v}(t)\sim\text{Ber}(0.5)$. We simulate multiple dynamics for different
memory ($d$) and $\bar{M}$, and three different causal or influence graph
structures ($G$): line, tree, and ring.

We run Algorithm 1 on data from all those different instants of the dynamics,
plot the dependence of the probability of graph recovery on the number of
samples, and observe the influence of $d$ and $\bar{M}$ on those. The
threshold $\epsilon/2$ in Algorithm 1 has to be chosen appropriately in each
case.

In Fig. 1a, we plot the probability of graph recovery against the number of
samples used for $d=1$, $\bar{M}=0$ (i.e $M_{v}(t)=1$ for all $v,t$). We plot
the same for $d=1$, $\bar{M}=2$ (i.e., $M_{v}(t)\leq 3$ for all $v,t$) in Fig.
1b and for $d=2$, $\bar{M}=2$ in Fig. 1c.

Note that from Fig. 1a to Fig. 1b, $d$ was fixed ($d=1$) and $M$ increases
from $0$ to $2$, where as from Fig. 1b to Fig. 1c, $\bar{M}$ remains the same
($\bar{M}=2$) and $d$ increases from $1$ to $2$. Interestingly, the number of
samples required for almost sure recovery changes from around $1000$ to $9000$
from Fig. 1a to Fig. 1b. On the other hand, the number of required samples
changes from $9000$ to $840000$ from Fig. 1b to Fig. 1c. In the second case,
the increase in the number of required samples is substantially more.

In these simulations, since the degree of all the nodes in all the graphs is
$1$, $|P|_{max}$ is actually a constant independent of $\bar{M}$. Hence, in
these cases, the analytical bound in Theorem 1 depends on $\bar{M}$ as
$\bar{M}^{2d}\log(\bar{M})$ through $|\xi|$. On the other hand, the analytical
bound has a term that is exponential in $d$. Thus, the dependence of the
constant terms (which are not dependent on $|V|$) of the bound on the
parameters $d$ and $\bar{M}$ corresponds well with the observations from
simulations.

Note that the bound in Theorem 1 is applicable when
$2(\bar{\mu}+L)|\rho(\tilde{A})|^{\frac{1}{d}}<1$. All the previous
simulations are done for a setting where that assumption is true. In Fig 2, we
plot the probability of success versus number of samples for $d=2$,
$\bar{M}=2$ when $2(\bar{\mu}+L)|\rho(\tilde{A})|^{\frac{1}{d}}>1$. We observe
that Algorithm 1 can still recover the true graphs perfectly but needs around
three times more samples than Fig. 1c.

For the above scenario, i.e.,
$2(\bar{\mu}+L)|\rho(\tilde{A})|^{\frac{1}{d}}>1$, we also observe the
dependence of the number of required samples on $d$ and $\bar{M}$. These are
plotted in Tables I and II. Note that here also we see that the number of
samples required for perfect recovery grows much faster with $d$ than with
$\bar{M}$ as in the analytical bound and as observed in Figures 1a–1c. This
seems to indicate that the $O(\log|V|)$ bound on sample complexity is likely
to be true even when $2(\bar{\mu}+L)|\rho(\tilde{A})|^{\frac{1}{d}}>1$. We
leave this for future analytical investigations.

TABLE I: Relationship of $d$ to the number of samples required for perfect recovery |  $\bar{M}$=0,d=1 |  $\bar{M}$=0, d=2 |  $\bar{M}$=0, d=3 |  $\bar{M}$=0, d=4
---|---|---|---|---
No. of samples | 3000 | 17000 | 48000 | 70000
TABLE II: Relationship of $\bar{M}$ to the number of samples required for perfect recovery | d=1,$\bar{M}$=0 | d=1,$\bar{M}$=1 | d=1,$\bar{M}$=2 | d=1,$\bar{M}$=3
---|---|---|---|---
No. of samples | 3000 | 4000 | 12000 | 51000

##  VIII Conclusion

We considered the problem of learning the influence or causal graph of a high-
dimensional Markov process with memory. This problem is vital in many real-
world applications, including social networks and nervous systems. We adapted
the RecGreedy algorithm for learning i.i.d graphical models from [13] by
introducing directed conditional entropy as the metric for learning the
neighborhood greedily. We proved a logarithmic sample complexity result by
deriving a bound on the second eigenvalue of the Markov process with memory,
which involved obtaining and comparing lower and upper bounds on the distance
from the stationary distribution as a function of time. Unlike RecGreedy [13],
where the non-degeneracy of conditional entropy for the i.i.d graphical model
is an unavoidable assumption, we presented a proof of non-degeneracy of the
conditional directed entropy in our Markov model. Simulations show that the
adapted RecGreedy algorithm has a high probability of success even in cases
where the analytical bound is not applicable.

## References

  * [1] A. De, I. Valera,N. Ganguly,S. Bhattacharya, Gomez Rodriguez “Learning and forecasting opinion dynamics in social networks.” Advances in neural information processing systems 29 (2016).
  * [2] C. Ravazzi, S. Hojjatinia, C. M. Lagoa and F. Dabbene, “Ergodic Opinion Dynamics Over Networks: Learning Influences From Partial Observations,” in IEEE Transactions on Automatic Control, vol. 66, no. 6, pp. 2709-2723, June 2021, doi: 10.1109/TAC.2021.3056362.
  * [3] C. Ravazzi, R. Tempo and F. Dabbene, “Learning Influence Structure in Sparse Social Networks,” in IEEE Transactions on Control of Network Systems, vol. 5, no. 4, pp. 1976-1986, Dec. 2018, doi: 10.1109/TCNS.2017.2781367.
  * [4] M. H. DeGroot, “Reaching a Consensus.” Journal of the American Statistical Association, vol. 69, no. 345, 1974, pp. 118–21. JSTOR, https://doi.org/10.2307/2285509. Accessed 9 Feb. 2024.
  * [5] A. Das, S. Gollapudi, and K. Munagala. “Modeling opinion dynamics in social networks.” Proceedings of the 7th ACM international conference on Web search and data mining. pp. 403-412. 2014.
  * [6] L. Cowen, T. Ideker,B.J. Raphael,R. Sharan, “Network propagation: a universal amplifier of genetic associations.” Nature Reviews Genetics 18, no. 9 (2017): 551-562.. https://doi.org/10.1038/nrg.2017.38
  * [7] Lazaro M. Sanchez-Rodriguez, Yasser Iturria-Medina, Pauline Mouches, Roberto C. Sotero, “Detecting brain network communities: Considering the role of information flow and its different temporal scales, NeuroImage, Volume 225, 2021, 117431, ISSN 1053-8119, https://doi.org/10.1016/j.neuroimage.2020.117431.
  * [8] Quinn CJ, Coleman TP, Kiyavash N, Hatsopoulos NG. Estimating the directed information to infer causal relationships in ensemble neural spike train recordings. J Comput Neurosci. 2011 Feb;30(1):17-44. doi: 10.1007/s10827-010-0247-2. Epub 2010 Jun 26. PMID: 20582566; PMCID: PMC3171872.
  * [9] G. Zhang et al., ”Estimating Dynamic Functional Brain Connectivity With a Sparse Hidden Markov Model,” in IEEE Transactions on Medical Imaging, vol. 39, no. 2, pp. 488-498, Feb. 2020, doi: 10.1109/TMI.2019.2929959.
  * [10] P. Ravikumar, M. J. Wainwright, and J. D. Lafferty. ”High-dimensional Ising model selection using l 1-regularized logistic regression.” The Annals of Statistics 38, no. 3 (2010): 1287-1319.
  * [11] Anandkumar, Animashree, and Ragupathyraj Valluvan. ”Learning loopy graphical models with latent variables: Efficient methods and guarantees.” The Annals of Statistics (2013): 401-435.
  * [12] P. Netrapalli, S. Banerjee, S. Sanghavi and S. Shakkottai, ”Greedy learning of Markov network structure,” 2010 48th Annual Allerton Conference on Communication, Control, and Computing (Allerton), Monticello, IL, USA, 2010, pp. 1295-1302, doi: 10.1109/ALLERTON.2010.5707063.
  * [13] A. Ray, S. Sanghavi and S. Shakkottai, ”Improved Greedy Algorithms for Learning Graphical Models,” in IEEE Transactions on Information Theory, vol. 61, no. 6, pp. 3457-3468, June 2015, doi: 10.1109/TIT.2015.2427354.
  * [14] J. Scarlett and V. Cevher, ”On the Difficulty of Selecting Ising Models With Approximate Recovery,” in IEEE Transactions on Signal and Information Processing over Networks, vol. 2, no. 4, pp. 625-638, Dec. 2016, doi: 10.1109/TSIPN.2016.2596439.
  * [15] G. Dasarathy, A. Singh, M.F. Balcan, and J. H. Park. ”Active learning algorithms for graphical model selection.” In Artificial Intelligence and Statistics, pp. 1356-1364. PMLR, 2016.
  * [16] L. Hamilton, F. Koehler, and A. Moitra, ”Information theoretic properties of Markov random fields, and their algorithmic applications.” Advances in Neural Information Processing Systems 30 (2017).
  * [17] D. Katselis, C. L. Beck and R. Srikant, ”Mixing Times and Structural Inference for Bernoulli Autoregressive Processes,” in IEEE Transactions on Network Science and Engineering, vol. 6, no. 3, pp. 364-378, 1 July-Sept. 2019, doi: 10.1109/TNSE.2018.2829520.
  * [18] W. Pieczynski, ”Pairwise Markov chains,” in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 25, no. 5, pp. 634-639, May 2003, doi: 10.1109/TPAMI.2003.1195998.
  * [19] B. Ribeiro, N. Chen, A. Kovacec, ”Shaping graph pattern mining for financial risk”, Neurocomputing, Volumes 326–327, 2019, Pages 123-131, ISSN 0925-2312,
  * [20] Y. Deng, Y. Liang, S.M. Yiu, “Towards interpretable stock trend prediction through causal inference”, Expert Systems with Applications, Volume 238, Part B, 2024, 121654, ISSN 0957-4174, https://doi.org/10.1016/j.eswa.2023.121654.
  * [21] R. Dahlhaus, Graphical interaction models for multivariate time series . Metrika 51, 157–172 (2000). https://doi.org/10.1007/s001840000055
  * [22] A. Jalali, and S. Sanghavi. ”Learning the dependence graph of time series with latent factors.” In Proceedings of the 29th International Coference on International Conference on Machine Learning, pp. 619-626. 2012.
  * [23] F. R. Bach and M. I. Jordan, ”Learning graphical models for stationary time series,” in IEEE Transactions on Signal Processing, vol. 52, no. 8, pp. 2189-2199, Aug. 2004, doi: 10.1109/TSP.2004.831032
  * [24] A. Chatterjee, A. S. Rawat, S. Vishwanath and S. Sanghavi, ”Learning the causal graph of Markov time series,” 2013 51st Annual Allerton Conference on Communication, Control, and Computing (Allerton), Monticello, IL, USA, 2013, pp. 107-114, doi: 10.1109/Allerton.2013.6736512.
  * [25] H. Hosseini, S. Kannan, B. Zhang and R. Poovendran, ”Learning Temporal Dependence from Time-Series Data with Latent Variables,” 2016 IEEE International Conference on Data Science and Advanced Analytics (DSAA), Montreal, QC, Canada, 2016, pp. 253-262, doi: 10.1109/DSAA.2016.34
  * [26] Thomas M. Cover, Joy A. Thomas, “Elements of information theory”
  * [27] J. Fan, B. Jiang, and Q. Sun. “Hoeffding’s Inequality for General Markov Chains and Its Applications to Statistical Learning.” Journal of Machine Learning Research 22.139 (2021): 1-35.
  * [28] N. Fourcaud, N. Brunel, “Dynamics of the firing probability of noisy integrate-and-fire neurons”, Neural Comput. 2002 Sep;14(9):2057-110. doi: 10.1162/089976602320264015. PMID: 12184844.
  * [29] N. Soltani and A. J. Goldsmith, ”Directed Information Between Connected Leaky Integrate-and-Fire Neurons,” in IEEE Transactions on Information Theory, vol. 63, no. 9, pp. 5954-5967, Sept. 2017, doi: 10.1109/TIT.2017.2700464.
  * [30] E. C. Hall, G. Raskutti and R. M. Willett, ”Learning High-Dimensional Generalized Linear Autoregressive Models,” in IEEE Transactions on Information Theory, vol. 65, no. 4, pp. 2401-2422, April 2019, doi: 10.1109/TIT.2018.2884673.
  * [31] Boucheron, Stéphane, Gábor Lugosi, and Pascal Massart, ’Basic Information Inequalities’, Concentration Inequalities: A Nonasymptotic Theory of Independence (Oxford, 2013; online edn, Oxford Academic, 23 May 2013), https://doi.org/10.1093/acprof:oso/9780199535255.003.0004,
  * [32] G. Bresler, D. Gamarnik and D. Shah. “Structure learning of antiferromagnetic Ising models.” Neural Information Processing Systems (2014).