1

Network Identification with Latent Nodes
via Auto-Regressive Models

arXiv:1601.04179v3 [] 7 Sep 2017

Erfan Nozari

Yingbo Zhao

Abstract—We consider linear time-invariant networks with
unknown topology where only a manifest subset of the nodes
can be directly actuated and measured while the state of the
remaining latent nodes and their number are unknown. Our
goal is to identify the transfer function of the manifest
subnetwork and determine whether interactions between
manifest nodes are direct or mediated by latent nodes. We
show that, if there are no inputs to the latent nodes, the
manifest transfer function can be approximated arbitrarily well
in the H∞ -norm sense by the transfer function of an
auto-regressive model and present a least-squares estimation
method to construct the auto-regressive model from measured
data. We show that the least-squares auto-regressive method
guarantees an arbitrarily small H∞ -norm error in the
approximation of the manifest transfer function, exponentially
decaying once the model order exceeds a certain threshold.
Finally, we show that when the latent subnetwork is acyclic, the
proposed method achieves perfect identification of the manifest
transfer function above a specific model order as the length of
the data increases. Various examples illustrate our results.

I. I NTRODUCTION
Network reconstruction problems are widespread in many
areas of science and engineering. In systems biology, for
instance, genetic network identification uses data from RNA
micro-array experiments to identify the interaction pattern
between genes in a regulatory network [2], [3]. In
neuroscience, researchers seek to understand how different
regions of the brain cooperate with each other by having
subjects perform certain goal-directed tasks while measuring
their brain activity via multi-channel recordings such as
electroencephalograms (EEG) [4]–[8]. Similar examples exist
in other areas including finance, social networks, and
physics. Roughly speaking, the objective in network
identification is to determine causal relationships among the
nodes in the network that model the direction and strength of
the interactions between them. While network control and
coordination has made much progress on problems where the
interaction topology is either given or the design objective
itself, not so much attention has been devoted to develop
techniques to address the identification of unknown
topologies from measured data. The need for the latter is
especially apparent in the context of complex, large-scale
networks, where it is often not possible to measure or actuate
all nodes, or even know their number. In this paper, we seek
to contribute to this body of work by studying the effect that
A preliminary version appeared as [1] at the American Control Conference.
E. Nozari and J. Cortés are with the Department of Mechanical and
Aerospace Engineering, University of California at San Diego, La Jolla, CA
92093, {enozari,cortes}@ucsd.edu. Y. Zhao is with Cymer Corporation, San Diego, CA 92127, yingbo.zhao@asml.com

Jorge Cortés

the presence of unmeasured nodes has on the identification
of networked linear systems with arbitrary topology.
Literature review: An increasing number of works study
topology identification problems to better understand the
interactions in large-scale networks and their role in
determining the network behavior. A complex network is
commonly represented as a directed graph, and the
interactions among neighboring nodes are represented by
directed edges whose weights reflect the interaction strength.
In this sense, topology identification aims at identifying the
adjacency matrix of the network graph [9] or its Boolean
structure [10]. The work [11] studies the complete
characterization of the interaction topology of consensus-type
networks using a series of node-knockout experiments,
where nodes are sequentially forced to broadcast a zero state
without being removed from the network. The work [12]
also uses node-knockout experiments to identify the topology
of directed linear time-invariant networks relying on the
cross-power spectral densities of the network response to
wide-sense stationary noise. The work [13] presents a
method to infer the topology of a network of coupled phase
oscillators from its stable response dynamics, assuming that
one can manipulate every individual node and perform large
number of experiments. In general, without such assumption,
it is difficult or impossible, depending on the additional
structural information available, to accurately identify the
topology of a general network. As a result, a main focus has
been on particular network realizations that explain the
measured data, such as the sparsest realization, sometimes
with a design parameter to manage the trade-off between
model accuracy and sparsity, see e.g., [3], [14]. Along these
lines, the work [15] considers the identification of networked
linear systems with tree topologies. The above-referenced
works rely on knowledge of the number of nodes in the
network. However, it is often impossible to sample the state
of all nodes, or even know the existence of some of them.
The work [16] studies the problem of learning latent tree
graphical models where samples are available only from a
subset of the nodes, and proposes computationally efficient
algorithms for learning trees without any redundant hidden
nodes. The work [17] proposes a method to identify the
latent nodes and consistently reconstruct the topology under
the assumptions that the network is a polytree and the degree
of each latent node is at least three, with out-degree of at
least two. Unlike the topology identification algorithms
proposed in [15], [17], our approach here allows for the
possibility of cycles in the network topology. Using the
notion of the dynamical structure function of a network with

2

latent nodes [18], the work [19] proposes a convex
optimization-based approach to find the best Boolean
structure among manifest nodes which consists of computing
and comparing the distance between an estimated transfer
matrix or data to all possible Boolean structures. The
problem of minimal state-space realization of a given
dynamical structure function was further studied in [20]. In
the present work, however, we use a least-square
autoregressive identification approach to identify not only
whether a pair of manifest nodes are dynamically connected,
but also whether this connection is direct or indirect
(latent-mediated) and, in the latter case, the length of the
shortest path between the two. Recent work has employed
sparse plus low-rank (S+L) decomposition to identify general
graphical models (with the possibility of cycles) with latent
variables for static [21] and dynamic [22] models. The
present paper has two main differences with respect to this
work. First, the S+L decomposition assumes that the
subnetwork among manifest nodes is sparse and the number
of latent nodes is (considerably) smaller than the number of
manifest ones, while our method is applicable to arbitrary
networks. Second, although the identification procedure
of [22] also leads to an auto-regressive (AR) model, it is
based on the so-called maximum-entropy covariance
extension. This method, with origins in seismic vibrations
and human voice analysis, seeks to maximize the prediction
error [23] (while our approach seeks to minimize it), leading
to very different models. Finally, our work is inspired by the
wide use in neuroscience of AR models to analyze brain
data via Granger causality and its variants and the study of
effective connectivity among different areas of the brain, see
e.g., [5], [6], [24]. The Granger causality measure is a
mainly descriptive tool that captures influence and
interconnection among time series. A popular variant of
Granger causality, direct directed transfer function
(dDTF) [8], [25] distinguishes between direct and indirect
interconnections between two nodes by multiplying the
directed transfer function (DTF, the normalized transfer
function between the two nodes) by the partial coherence
between them in the frequency domain. We are motivated
here by understanding to what extent the reconstruction
results obtained via methods that build on Granger causality
are sensitive to the presence of latent nodes. Furthermore, we
propose a method using (multivariate) AR models for
networks with latent nodes that distinguishes between direct
and indirect (i.e., latent-mediated) interconnections between
two nodes in the time domain based on the order of the
interconnection between them.
Statement of contributions: We consider a scenario where
one can only directly actuate and measure a subset of the
nodes, termed manifest, of a large linear time-invariant
network whose total number of nodes and interaction
topology are unknown. The objective is to identify the
manifest transfer function, which is the submatrix
corresponding to the manifest nodes of the transfer function
matrix of the entire network. To achieve this, we study the
transfer functions provided by linear AR models. Our
discussion shows how AR models can be used to effectively

distinguish direct interactions between manifest nodes from
indirect interactions mediated by latent nodes. Our first
contribution shows that, if no inputs act on the latent nodes,
then there exists a class of AR models whose transfer
functions converge exponentially in the H∞ norm to the
manifest transfer function as the model order increases. We
also show that, if the latent subnetwork is acyclic, then this
approximation is exact above a specific model order. Our
second contribution characterizes the properties of using
least-squares auto-regressive estimation to construct the AR
model from measured data. We establish that the
least-squares matrix estimate converges in probability to the
optimal matrix sequence identified in our first contribution,
enabling us to determine whether two manifest nodes
interact directly or indirectly through latent nodes. We also
show that the least-squares auto-regressive method
guarantees an arbitrarily small H∞ -norm error as the length
of data and the model order grow. In fact, once the order of
the AR model candidates exceeds a certain threshold, the
H∞ -norm error decays exponentially. Finally, we show that,
when the latent subnetwork is acyclic, the method achieves
perfect identification of the manifest transfer function.
Throughout a series of remarks in the paper, we also discuss
how our results can be extended to the identification of
linear network models of arbitrary order. Simulations on a
directed ring network, Erdős–Rényi random graphs, and real
EEG data illustrate our results.
Notation: For a vector x ∈ Rn , we use xi to denote its
n
i-th element. Given a sequence {x(k)}∞
and
k=0 ⊂ R
j2
j1 ≤ j2 ∈ Z≥0 , we use {x}j1 to denote the finite sequence
{x(j1 ), x(j1 + 1), . . . , x(j2 )}. We omit j1 if j1 = 0. We
 21
Pj2
T
denote k{x}jj21 k ,
. A sequence of
k=j1 x (k)x(k)
random variables {x} converges in probability to a random
variable X, denoted plimk→∞ x(k)
=
X, if
limk→∞ Pr(|x(k) − X| ≥ ε) = 0 for all ε > 0. Accordingly,
a sequence of random matrices {A} converges to a random
matrix A∞ in probability if plimk→∞ Aij (k) = A∞,ij for
all i, j. For a real matrix M ∈ Rm×n , we denote its singular
values in decreasing order as σ 1 (M ) ≥ σ 2 (M ) ≥ · · ·
≥ σ min(m,n) (M ) ≥ 0 and its spectral norm by
kM k
=
σ 1 (M ). The max norm of M is
kM kmax = maxi,j |Mij |. We denote by ρ(M ) the spectral
radius of a square matrix M . The matrix M is Schur stable
if and only if ρ(M ) < 1. We let 0m×n denote the m × n
matrix with all zero elements and by In the identity matrix
of dimension n × n. The H∞ -norm of a discrete transfer
function T is kT k∞ , sup−π≤ω≤π kT (ω)k.
II. P ROBLEM FORMULATION
We consider a discrete-time, linear time-invariant (LTI)
network dynamics with state-space representation
x(k + 1) = Ax(k) + u(k),
y(k) = Cx(k),

(1)
n

where k ∈ Z≥0 is the time index, x(k) ∈ R is the network
state (with xi (k) representing the state of node
i ∈ {1, . . . , n}), u(k) ∈ Rn is the control input (with ui (k)
acting on node i), and y(k) ∈ Rm is the network output.

3

Here, A ∈ Rn×n is the adjacency matrix of the network,
characterizing the interactions among neighboring nodes, and
C ∈ Rm×n is the output matrix. Since natural system are
usually driven by noise, the input, state, and output
sequences are in general stochastic processes over the
sample space of noise realizations. For simplicity, the
dynamical description (1) assumes that all nodes are of order
1, that is, x(k + 1) depends directly only on x(k) and is
conditionally independent of {x}k−1 given x(k).
Nevertheless, as we discuss later (see e.g., Remark 3.4), all
of the subsequent results are generalizable to systems whose
dynamics (in the original “physical” variables) are described
by difference equations of order higher than 1.
Even though there is a control input at every node in the
network dynamics (1), we do not assume that all the control
inputs are user-specified. In fact, in a large-scale network, it is
common that one can actuate only a small subset of the nodes
due to computational constraints, physical limitations, or cost.
A similar observation can be made regarding the number of
nodes whose state can be directly measured. For these reasons,
here we assume that the nodes of the network are divided
into nm ≤ n manifest nodes, which can be directly actuated
and measured by the user, and n − nm latent nodes, which
can neither be directly actuated nor measured by the user.
With this distinction, and using a permutation of the indices
in (1, 2, . . . , n) if necessary, we can decompose the network
and input state as x = [xm , xl ] and u = [um , ul ], respectively,
where the subindex ‘m’ corresponds to manifest nodes and the
subindex ‘l’ corresponds to latent nodes. With this convention,
the output matrix takes the form C = [Inm ×nm , 0nm ×(n−nm ) ].
With the decomposition of the nodes into manifest and latent,
the state-space representation (1) becomes

 

 

xm (k + 1)
A11 A12 xm (k)
um (k)
=
+
,
xl (k + 1)
A21 A22 xl (k)
ul (k)
y(k) = xm (k).

(2)

In the remainder of this paper, we consider the network in the
relabeled form (2). Fig. 1 illustrates this relabeling procedure
(corresponding to a linear transformation) in a ring.
u1

um,1

1

4

1

x m,1

x1

u4

u2

relabeling
2

4

ul ,2

ul ,1

3

x m,2

x3
3

u3

2

um,2

Fig. 1. Node relabeling in a directed ring with 4 nodes. Nodes 1 and 3 are
manifest, nodes 2 and 4 are latent. The permutation (1, 2, 3, 4) → (1, 3, 2, 4)
makes manifest and latent nodes have consecutive indices, as in (2).

Since the focus of this work is on network identification and
not stabilization, we make the following standard assumption.
Assumption 2.1: The adjacency matrix of the complete
network as well as the latent subnetwork are Schur stable,
i.e., ρ(A) < 1 and ρ(A22 ) < 1.
Remark 2.2: (Direct versus latent interactions). The
interaction graph of the manifest subnetwork is characterized

by A11 . In particular, the state of node p affects the state of
node q directly if and only if the entry on the q-th row and
the p-th column, denoted by A11 (q, p), is nonzero. However,
even if A11 (q, p) = 0, it is still possible that node p affects
node q indirectly through some latent nodes. The distinction
between direct and indirect connections is an important point
to which we come back later in our discussion.

We refer to a latent node as passive if its corresponding
input is zero. Throughout the paper, we only deal with
passive latent nodes, so that {ul } ≡ 0. We make the
following assumption on the input to the manifest nodes.
Assumption 2.3: The input {um } to the manifest
subnetwork is a zero-mean stochastic process with
independent and identically distributed (i.i.d.) absolutely
continuous1 random vectors um (k), with covariance Inm .
Assumption 2.3 guarantees that {um } is persistently
exciting of arbitrary order and its power spectral density
does not vanish at any frequency. Similar assumptions are
common in system identification, see e.g., [12], [26]. The
zero-mean assumption can be relaxed by assuming a nonzero
but known E[um (k)] corresponding to the scenario where the
designer injects a deterministic stimulating signal into every
manifest node, which itself is subject to the disturbance of a
zero-mean white noise. Without loss of generality and for
simplicity, we assume E[um (k)] ≡ 0nm .
Given the setup above, our objective is to identify the
transfer function Txm um (ω) of the manifest subnetwork, that
is, the transfer function from um to xm , absent any
knowledge of the latent nodes.
Problem 2.4: (Identification of the manifest transfer
function). Given the measured data {y}N
1 , find a linear
auto-regressive model of order τ , with N  τ , of the form
x̃m (k + 1) =

τ
−1
X

Ãi x̃m (k − i) + um (k),

(3)

i=0

such that the associated transfer function Tx̃m um from um to
x̃m and the transfer function Txm um from um to xm in (1)
are close in the H∞ -norm, i.e., kTx̃m um − Txm um k∞ is small.
There are alternative methods to identify the transfer
function matrix Txm um besides the AR method in (3). Our
adoption here of AR model candidates is motivated by their
widespread use in neuroscience to determine causality and
interconnections in human brain connectivity models, see
e.g., [5]–[7]2 . Equipped with time series data obtained during
the performance of a cognitive task, the conventional
procedure consists of first estimating an AR model, then
computing its associated transfer function matrix, and finally
evaluating the Granger causality connectivity measure, or
generalizations of it, in the frequency domain. We are
particularly motivated by the prospect of understanding the
sensitivity of these approaches to the presence of latent
1 Recall that an absolutely continuous random variable/vector is one that
has a probability density function (e.g., Gaussian).
2 In general, the main advantage of AR models over more general models
such as ARMA or BJ is their simplicity, only capturing the internal dynamics
and assuming negligible input noise correlation (though putting no restriction
on input signal correlation, which is significant in brain dynamics). As a
result, prediction error minimization has a closed-form solution for an AR
model while it is non-convex in the ARMA or BJ cases.

4

nodes corresponding to brain regions that are active during
the cognitive task but are not directly measured.
III. A SYMPTOTICALLY EXACT IDENTIFICATION OF THE
MANIFEST TRANSFER FUNCTION

In this section we establish that, given an arbitrary
precision, there exists an AR model solving Problem 2.4.
More precisely, we show that there exists a sequence of AR
models of the form (3) with increasing order whose transfer
functions converge to Txm um exponentially in the H∞ sense.
We later show that, if the latent subnetwork is acyclic, then
this approximation can be made exact.
We start our discussion with a useful auxiliary result.
Lemma 3.1: (Upper bound on kAi22 k). For any Schur stable
A22 ∈ Rnl ×nl and any ρ̄ ∈ (ρ(A22 ), 1), there exists κ ∈ R>0
such that kAi22 k ≤ κ · ρ̄i , for all i ∈ Z≥0 .
Proof. The result is an immediate consequence of the spectral
radius formula limi→∞ kAi22 k1/i = ρ(A22 ).
We are now ready to state the main result of this section.
Theorem 3.2: (AR model whose transfer function
converges to the manifest transfer function). Consider the
LTI network described by (2) where all the latent nodes are
passive. For any ρ̄ ∈ (ρ(A22 ), 1), there exists γ̄ ∈ R>0 such
that for all τ ∈ Z≥0 , the AR model (3) with
Ã∗0

= A11 ,

Ã∗i

=

A12 Ai−1
22 A21 ,

i ∈ {1, . . . , τ − 1}, (4)

guarantees
kTx̃m um (·, τ ) − Txm um k∞ ≤ γ̄ · ρ̄τ .

(5)

Proof. We obtain from (2) that
Txm um (ω) = (zInm − A11 − A12 (zInl − A22 )−1 A21 )−1
∞
X
(a)
−1
= (zInm − A11 −
z −i A12 Ai−1
, (6)
22 A21 )
i=1

where z = P
ejω and (a) follows by using the relation (zInl
∞
−1
A22 ) = i=1 z −i Ai−1
22 . Similarly, from (3) we obtain
Tx̃m um (ω, τ ) = (zInm −

τ
−1
X

z −i Ã∗i )−1 .

−

(7)

i=0

Here we write the transfer function as Tx̃m um (ω, τ ) to
emphasize its dependence on τ . It then follows directly that
kTx̃m um (·, τ ) − Txm um k∞
=

kTxm um (Tx−1
m um

−

(a)

≤ kTxm um k∞ kTx̃m um (·, τ )k∞ kTx−1
− Tx̃−1
(·, τ )k∞
m um
m um
∞
X
(b)
≤ kTxm um k∞ kTx̃m um (·, τ )k∞
kz −i A12 Ai−1
22 A21 k∞
i=τ

≤ kTxm um k∞ kTx̃m um (·, τ )k∞ kA12 kkA21 k

∞
X

kAi−1
22 k

i=τ
(d)

≤ γ(τ ) · ρ̄τ ,

where
γ(τ ) ,

κkTxm um k∞ kA12 kkA21 k
kTx̃m um (·, τ )k∞ .
ρ̄ − ρ̄2

−π≤ω≤π

(a)

=



inf

−π≤ω≤π

−1
σ min Tx̃−1
(ω, τ )
,
m um

−1
where (a) holds due to the fact that σ max (M ) = σ −1
)
min (M
for any invertible matrix M . To complete the proof, we only
need to show that

ϑ , inf
inf σ min Tx̃−1
(ω, τ ) > 0.
(10)
m um
τ ∈Z≥0 −π≤ω≤π

We show this in two steps.
(i) It follows from (6) and (7) that
lim Tx̃−1
(ω, τ ) = Tx−1
(ω),
m um
m um

τ →∞

∀ω ∈ [−π, π].

It is straightforward to show, using the exponential
decay of Aτ22 and definition of uniform convergence,
(·, τ ) converges uniformly to
that each entry of Tx̃−1
m um
the corresponding entry of Tx−1
. Hence, given the
m um
uniform continuity of matrix eigenvalues as a function
−1
of matrix entries [27, Thm 7.8c], σ min
 Tx̃m um (·, τ )
−1
converges uniformly to σ min
 Txm um . Thus, since
(ω)
= kTxm um k∞ > 0
inf −π≤ω≤π σ min Tx−1
m um
(which itself holds by Assumption 2.1), there exists
τ 0 ∈ Z≥0 such that

inf σ min Tx̃−1
(ω, τ ) > 0.
inf
m um
τ ≥τ 0 −π≤ω≤π

(ii) For any finite τ , we show that Tx̃m um (·, τ ) is BIBO
stable and thus has no poles on the unit circle (which
in

turn guarantees inf −π≤ω≤π σ min Tx̃−1
(ω,
τ
)
>
0).
m um
For any bounded input um , let the corresponding outputs
of Tx̃m um (·, τ ) and Txm um be denoted by x̃m and xm ,
resp. (with initial states set to zero). We then have
xm (k) − x̃m (k) = A12 Aτ22−1 xl (k − τ ),

Tx̃−1
(·, τ ))Tx̃m um (·, τ )k∞
m um

(c)

Here, (a) follows from the sub-multiplicativity of induced
norms, (b) follows by the sub-additivity of norms, (c)
follows by the definition of the H∞ -norm and also the
sub-multiplicativity of induced norms, and (d) follows from
Lemma 3.1. The remainder of the proof is devoted to
showing the existence of a uniform upper bound γ̄ for γ(τ ).
By the definition of the H∞ -norm,

kTx̃m um (·, τ )k∞ = sup σ max Tx̃m um (ω, τ )
(9)

(8)

where xl is the (internal) state of the latent nodes in
Txm um . By Assumption 2.1, both xm (k) and
A12 Aτ22−1 xl (k − τ ) are bounded, proving the BIBO
stability of Tx̃m um (·, τ ).
Hence, (10) follows by combining (i) and (ii) and the fact
that the decomposition Z≥0 = {0} ∪ {1} ∪ · · · ∪ {τ 0 − 1} ∪
{τ 0 , τ 0 + 1, . . . } is finite. Equivalently, there exists U > 0
such that kTx̃m um (·, τ )k∞ < U for all τ ∈ Z≥0 , so (5) holds
with γ̄ = κU kTxm um k∞ kA12 kkA21 k/(ρ̄ − ρ̄2 ).
Theorem 3.2 shows that the presence of latent nodes in the
network, as long as they do not receive any external input,
does not affect the achievable accuracy of the identification
via AR modeling of the manifest transfer function.

5

Remark 3.3: (Direct versus latent interactions – cont’d). It
follows from the network dynamics (2) that
xm (k + 1) =

k
X

Ã∗i xm (k − i) + A12 Ak22 xl (0) + um (k). (11)

i=0

By virtue of (11), we can distinguish whether two manifest
nodes interact directly or indirectly through latent nodes by
looking at the matrix sequence {Ã∗i }. First, the state of
manifest node p affects the state of manifest node q directly
if and only if Ã∗0 (q, p) = A11 (q, p) 6= 0. Similarly, the state
of manifest node p affects the state of manifest node q
indirectly through latent nodes if and only if Ã∗i (q, p) 6= 0
for some i ≥ 1. In particular, from the relation
Ã∗i = −A12 Ai−1
22 A21 , one can see that the state of p first
affects some latent nodes (that correspond to the nonzero
entries in the p-th column of A21 ) through A21 , then
propagates through the latent subnetwork, reflected by Ai−1
22 ,
and finally affects q through A12 . Furthermore, if the latent
subnetwork is acyclic, then Ã∗i (q, p) 6= 0 implies that there
are exactly i latent nodes in a path connecting p to q.

Remark 3.4: (Systems described by higher-order
difference equations). Unlike the system description in (1),
the dynamic behavior of many real-world complex systems
such as the brain cortical networks is described by difference
equations of orders significantly greater than 1, i.e.,
x(k + 1) = A(0) x(k) + A(1) x(k − 1) + · · ·
+A

(ν−1)

x(k − ν + 1) + u(k),

(12)
ν1

where x1 , . . . , xnm still denote the manifest (sensed and
actuated) nodes and xnm +1 , . . . , xn are the latent ones. In
this description, the vector x corresponds to some relevant
physical
variables.
Defining
the
state
vector
ξ(k) = [x(k)T x(k − 1)T · · · x(k − ν + 1)T ]T , one can
rewrite (12) in order-1 form as

 

 

ξ m (k + 1)
A11 A12
ξ m (k)
um (k)
=
+
,
ξ l (k + 1)
A21 A22
ξ l (k)
0
(13)
where ξ m (k) = xm (k), ξ l (k) = [xl (k)T xm (k − 1)T xl (k −
(0)
1)T · · · xm (k − ν + 1)T xl (k − ν + 1)T ]T , A11 = A11 , and
i
h
(1)
(τ −1)
(τ −1)
(1)
,
A12 = A(0)
A
A
A
·
·
·
A
11
12
11
12
12
h
iT
T
A21 = (A(0)
,
I nm 0 · · · 0 0
21 )
 (0) (1) (1)

(τ −2) (τ −2) (τ −1) (τ −1)
A21
A22
A22 A21 A22 · · · A21
A22
 0
0
0 ···
0
0
0
0 


 In 0
0
·
·
·
0
0
0
0 
 l

 .
..
.. . .
..
..
..
.. 
 ..
.
.
.
.
.
.
. 
A22 = 
.
 0
0
0 ···
0
0
0
0 


 0
0
0 ···
0
0
0
0 


 0
0
0 · · · Inm
0
0
0 
0
0
0 ···
0
I nl
0
0
In this description, we view ξ m as the actual “manifest
state” of the system while the whole vector ξ l is the “latent
state”. The reason for this interpretation is that, at any time
k, only xm (k) is directly sensed/actuated while

x(k − 1), . . . , x(k − ν + 1) are quantities stored in the
system. Interestingly, for the order-1 description (1), this
observation brings up the possibility of some of the latent
variables xl simply being a relayed version of manifest
variables. Note that, under this interpretation, the matrices
(1)
(ν−1)
A11 , . . . , A11
represent manifest-latent (rather than
manifest-manifest) interactions. From (13), it is clear that all
the treatment for (1) is readily applicable. Nevertheless, as ν
increases, larger τ is necessary in order for (3) to represent
the system accurately. This is both intuitive and clear
from (5) and (8), where increasing ν results in larger kA12 k
and kA21 k as well as (usually) kTxm um k and ρ(A22 ). This,
in turn, may result in numerical difficulties when one
constructs the AR model from recorded input-output data
(which is the subject of the next section).

Next, we show that there exists an AR model (3) whose
transfer function coincides with the manifest transfer function
if the latent subnetwork is acyclic.
Corollary 3.5: (Exact manifest transfer function
identification for acyclic latent subnetworks). Under the
assumptions of Theorem 3.2, further assume that the latent
subnetwork is acyclic, i.e., there exists τ 22 ∈ Z≥1 such that
Aτ2222 = 0nl ×nl . Then, the matrix sequence Ã∗0 , · · · , Ã∗τ 22
in (4) ensures Tx̃m um = Txm um .
The proof of the result follows by comparing (6) and (7),
and using the assumption that the latent subnetwork is
acyclic. Theorem 3.2 and Corollary 3.5 show that it is
possible to identify the transfer function of the manifest
subnetwork without any knowledge of the passive latent
nodes. However, (4) cannot be directly applied to determine
the auto-regressive model because its evaluation requires
knowledge of the adjacency matrix A of the whole network,
which is unknown. This problem can be circumvented by
nm
employing the measured data sequence {y}N
, as
1 ⊂ R
explained in the next section.
IV. I DENTIFICATION VIA LEAST- SQUARES ESTIMATION
In this section we employ least-squares estimation to
compute from data the sequence of matrices defining the
auto-regressive model. We show that the estimates resulting
from this method asymptotically converge in probability, as
the data length N and model order τ increase, to the optimal
matrix sequence identified in Theorem 3.2. Finally, we
particularize our discussion to the case of acyclic latent
subnetworks.
A. Least-squares auto-regressive estimation
nm
Given a vector sequence {y}N
, the problem of least1 ⊂R
squares auto-regressive (LSAR) model estimation with order
τ ∈ Z≥1 is to find a matrix sequence {Â}τ0 −1 ⊂ Rnm ×nm
−1
that minimizes the 2-norm of the residual sequence {e}N
⊂
τ
nm
R
defined by
e(k) = y(k + 1) −

τ
−1
X

Âi y(k − i),

(14)

i=0

for k ∈ {τ , . . . , N − 1}. Equation (14) can be written in
compact vector form as
~yN = Âτ ΦN + ~eN ,

(15)

6

Proof. For any quasi-stationary signal4 {s}, let

where

~yN = y(τ + 1)

~eN = e(τ ) e(τ

Âτ = Â0 Â1

y(τ )
y(τ − 1)

ΦN = 
..

.
y(1)

y(τ + 2)
+ 1)

y(N ) ∈ Rnm ×(N −τ ) ,

e(N − 1) ∈ Rnm ×(N −τ ) ,

···

···

· · · Âτ −1
y(τ + 1)
y(τ )
..
.
y(2)





∈ Rnm ×nm τ ,

···
···
..
.
···

N
1 X
E[s(i)s(i − j)T ].
N →∞ N
i=1

Rs (j) , lim

Using the Birkhoff’s Ergodic Theorem [28, Thm 7.2.1] (see
also [28, Thm 7.1.3]) and the fact that {y} is the output of a
stable system (and thus the effects of initial conditions
asymptotically vanish), we can show that


y(N − 1)
y(N − 2)

.
..

.
y(N − τ )

N
1 X
y(i)y(i − j)T = Ry (j).
N →∞ N i=1

plim

Using the square of the prediction error [26]
tr(~eN ~eTN ) = tr (~yN − Âτ ΦN )(~yN − Âτ ΦN )T


As a result,
bility and

as the cost function, we compute its gradient
∂ tr(~eN ~eTN )

= (~yN − Âτ ΦN )(−ΦTN ) = Âτ ΦN ΦTN − ~yN ΦTN .

∂ Âτ
Setting this to zero, we get a system of linear equations for
which a solution is guaranteed to exist (since the rows
of ~yN ΦTN belong to the row space of ΦN ΦTN , which is the
same as the row space of ΦTN ). By Assumption 2.3,
det(ΦN ΦTN ) 6= 0 and this solution is unique with probability
one.3 If det(ΦN ΦTN ) = 0, the minimum-norm solution can
be found as
Âτ =

~yN ΦTN (ΦN ΦTN )−1

=

~yN Φ+
N,

1
T
N ΦN ΦN

∈ Rnm τ ×nm τ also converges in proba-

1
RΦ , plim ΦN ΦTN
N →∞ N

Ry (0)
Ry (1)
···
 RyT (1)
R
(0)
···
y

=
..
..
..

.
.
.
RyT (τ − 1) RyT (τ − 2) · · ·
Define
ν(k) , y(k + 1) −

(16)

τ
−1
X


Ry (τ − 1)
Ry (τ − 2)

.
..

.
Ry (0)

Ã∗i y(k − i),

(18)

i=0

where (·)+ denotes the Moore-Penrose pseudo-inverse.
Since (16) is also valid for the nonsingular case, it is taken
as the solution to the LSAR estimation problem. In order to
indicate the dependency of the solution upon the measured
data sequence, we sometimes use the notation Âτ ({y}N
1 ).

and note that the transfer function from um to ν is
Tx̃−1
T
, where Txm um and Tx̃m um are given by (6)
m um xm um
and (7), respectively. Equation (18) can be written in
compact vector form as

B. Convergence in probability to manifest transfer function
Here we study the transfer function resulting from the
LSAR estimation method and characterize its convergence
properties, as the data length and the model order increase,
with respect to the transfer function of the manifest
subnetwork. Our first result establishes that the LSAR matrix
estimate (16) converges in probability to the optimal matrix
sequence identified in Theorem 3.2.
Proposition 4.1: (The LSAR estimate converges in
probability to optimal matrix sequence). Consider the LTI
network described by (2) where all latent nodes are passive.
Given the measured data sequence {y}N
1 generated from the
dynamics (2) stimulated by the white noise input {um }
according to Assumption 2.3 and any ρ̄ ∈ (ρ(A22 ), 1), there
exists β ∈ R>0 (depending only on the adjacency matrix A)
such that the LSAR estimate Âτ ({y}N
1 ) in (16) satisfies

with ~ν N , [ν(τ ) ν(τ + 1) · · · ν(N − 1)] ∈ Rnm ×(N −τ ) .
From (16) and (19), it follows that

∗
τ
k plim Âτ ({y}N
1 ) − Ãτ kmax ≤ βτ ρ̄ ,

(17)

N →∞



where Ã∗τ = Ã∗0 Ã∗1 · · · Ã∗τ −1 ∈ Rnm ×nm τ is the optimal matrix sequence given by (4).
−1
is because (each element of) {y}N
is an affine function of
1
T
and det(ΦN ΦN ) is a polynomial function of {y}1N −1 , so
N −2
det(ΦN ΦT
. Therefore, the level set
N ) is a polynomial function of {u}0
N −2
T
N = {{u}0
| det(ΦN ΦN ) = 0} has Lebesgue measure zero. Thus, by
Assumption 2.3, Pr(N ) = 0.
3 This

{u}0N −2 ,

~yN = Ã∗τ ΦN + ~ν N ,

(19)

1
1
~yN ΦTN ( ΦN ΦTN )−1
N
N
1
1
= Ã∗τ + plim ~ν N ΦTN ( ΦN ΦTN )−1 . (20)
N
N
N →∞

plim Âτ ({y}N
1 ) = plim

N →∞

N →∞

Moreover, Assumption 2.3 renders um (k) independent of
{y}k1 ,
which
further
implies
that
plimN →∞ N1 ~um,N ΦTN
=
0nm ×nm τ ,
where
~um,N , [um (τ ) um (τ + 1) · · · um (N − 1) ∈ Rnm ×(N −τ ) .
Therefore,
1
1
~ν N ΦTN = plim (~ν N − ~um,N )ΦTN = Ψ,
N
N →∞ N


where Ψ , Ψ1 Ψ2 · · · Ψτ ∈ Rnm ×nm τ , with
plim

(21)

N →∞

N −1
1 X
(ν(i) − um (i))y T (i − j + 1) ∈ Rnm ×nm .
N →∞ N i=τ

Ψj , plim

−1
Thus, using plimN →∞ ( N1 ΦN ΦTN )−1 = RΦ
, we have
−1
∗
plim Âτ ({y}N
1 ) − Ãτ = ΨRΦ .

N →∞

4 Basically, a signal is quasi-stationary if it has a well-defined covariance
function. See [26, Def 2.1] for a formal definition.

7

By the sub-additivity of the max norm, it holds for any j ∈
{1, . . . , τ } that
N −1
1 X
k(ν(i) − um (i))y T (i − j + 1)kmax
N →∞ N i=τ

kΨj kmax ≤ plim
(a)

N −1
ρ̄−τ X
(ν(i)−um (i))T (ν(i)−um (i))
N →∞ N i=τ

≤ plim

N −1
ρ̄τ X T
y (i − j + 1)y(i − j + 1)
N →∞ N i=τ

Given the result in Proposition 4.1, we next turn our
attention to the transfer function from e to y resulting from
the LSAR estimation (14), which we denote by
Tye ({y}N
1 , τ ). The next result shows that the H∞ -norm of
this transfer function is uniformly upper bounded with
respect to the model order τ .
Lemma 4.3: (H∞ -norm of Tye is uniformly upper
bounded). Under the assumptions of Proposition 4.1, there
exist positive scalars τ 0 and UT∞ye such that, for τ ≥ τ 0 ,

+ plim
= ρ̄

−τ

τ

tr(Rν−um (0)) + ρ̄ tr(Ry (0)),

∞
k plim Tye ({y}N
1 , τ )k∞ ≤ UTye .

(22)

where (a) follows from Lemma A.1 in the appendix with the
positive scalar M chosen as ρ̄τ . Using the fact that the transfer
function from um to ν −um is Tx̃−1
T
−Inm , we obtain
m um xm um

Proof. By definition of H∞ -norm, we have
k plim Tye ({y}N
1 , τ )k∞ =
N →∞

=

N
−1
X

1
E[(ν − um )(i)(ν − um )T (i)]
N →∞ N
Z π i=0
(a) 1
(T −1 Tx u (ω) − Inm )
=
2π −π x̃m um m m
· (Tx̃−1
T
(ω) − Inm )∗ dω
m u m xm u m

Rν−um (0) , lim



inf

−π≤ω≤π

sup

(23)

2
where γ̂ , γ̄ 2 1 + kA11 k + kA12 kkA21 kκ(1 − ρ̄)−1 is
constant, (a) follows from [29, eq. (9-193)], (b) follows by
the definition of H∞ -norm, (c) follows by the
sub-multiplicativity of induced norms, and (d) holds because
of Theorem 3.2 and the observation that, by Lemma 3.1,
kTx̃−1
k ≤ 1 + kA11 k + kA12 kkA21 kκ(1 − ρ̄)−1 .
m um ∞
We obtain from (22) and (23),
kΨj kmax ≤ ρ̄τ (γ̂nm + tr(Ry (0))),

N →∞

N →∞

τ
−1
X

z −i Âi

i=0

= Tx̃−1
(ω, τ ) −
m um

k plim Âτ ({y}N
1 )−
N →∞

−1
Ã∗τ kmax =kΨRΦ
kmax

−1
≤ nm τ kRΦ
kmax kΨkmax
−1
= nm τ kRΦ
kmax maxkΨj kmax ≤ βτ ρ̄τ ,
j

−1
where β = (γ̂n2m + tr(Ry (0))nm )kRΦ
kmax , as claimed.
When it is clear from context, we refer to
plimN →∞ Âi ({y}N
1 ) simply as Âi .
Remark 4.2: (Direct versus latent interactions – cont’d).
Proposition 4.1 shows that Âi converges in probability to Ã∗i
exponentially as the model order τ increases. Therefore,
within a margin of error that can be tuned as desired, we
deduce from the discussion in Remark 3.3 that the LSAR
estimate Â0 allows us to determine whether two manifest
nodes interact directly and the LSAR estimates {Âi }i≥1
allow us to determine whether two manifest nodes interact
indirectly through latent nodes with high probability as the
length of measurement data grows.


τ
−1
X

z −i (Âi − Ã∗i ), (26)

i=0

where z = ejω . However, for every ω ∈ [−π, π] and τ ∈ Z≥0 ,
τ
−1
X

z −i (Âi − Ã∗i ) ≤

τ
−1
X

(a)

kÂi − Ã∗i k ≤ nm

i=0

i=0

τ
−1
X

kÂi − Ã∗i kmax

i=0

(b)

≤ nm τ maxkÂi −
i

Ã∗i kmax

≤ nm βτ 2 ρ̄τ ,

where (a) follows from the fact that kAk ≤ nm kAkmax for
any matrix A ∈ Rnm ×nm and (b) follows from
Proposition 4.1. Therefore, using Weyl’s theorem for the
perturbation of singular values [30] in (26) and taking
inf −π≤ω≤π of both sides, we get

−1
inf σ min plim Tye
(ω, τ )
−π≤ω≤π

and from (20) and (21),

N →∞

−π≤ω≤π

(c)

≤ γ̂ρ̄2τ Inm ,



−1
−1
σ min plim Tye
(ω, τ )
. (25)

−1
plim Tye
(ω, τ ) = zInm −

T
− Inm k2∞ Inm
≤ kTx̃−1
m u m xm u m

(d)

σ max plim Tye (ω, τ )

Note that, for every ω ∈ [−π, π] and τ ∈ Z≥0 ,

(b)

≤ kTxm um − Tx̃m um k2∞ kTx̃−1
k2 I
m um ∞ nm

(24)

N →∞

≥
≥

N →∞

inf

τ
−1
X

σ min Tx̃−1
(ω,
τ
)
−
z −i (Âi − Ã∗i )
m um

inf

i=0

2 τ
σ min Tx̃−1
(ω,
τ
)
−
n
m βτ ρ̄ .
m um

−π≤ω≤π

−π≤ω≤π

In view of (10), let τ 0 be such that
nm βτ 2 ρ̄τ ≤

ϑ
,
2

∀τ ≥ τ 0 .

(27)

Then, the result follows from (25) with UT∞ye = ϑ2 .
We are finally ready to show that the transfer function Tye
obtained from the LSAR method converges in probability to
the transfer function Txm um of the manifest subnetwork.
Theorem 4.4: (The LSAR method consistently estimates
the manifest transfer function). Under the assumptions of
Proposition 4.1, for any ρ̄ ∈ (ρ(A22 ), 1), there exist positive
scalars β̄, γ̄ and τ 0 such that, for τ ≥ τ 0 ,
2
τ
k plim Tye ({y}N
1 , τ ) − Txm um k∞ ≤ (β̄τ + γ̄)ρ̄ .
N →∞

(28)

8

Consequently, plimN →∞,τ →∞ Tye ({y}N
1 , τ ) = Txm um .
Proof. We only need to prove (28) as it directly implies the
last equation in the statement. By the sub-additivity and submultiplicity of induced norms,
kTye (·, τ ) − Txm um k∞
≤ kTye (·, τ ) − Tx̃m um (·, τ )k∞ + kTx̃m um (·, τ ) − Txm um k∞
−1
≤ kTye (·, τ )k∞ kTx̃m um (·, τ )k∞ kTye
(·, τ ) − Tx̃−1
(·, τ )k∞
m um

+ kTx̃m um (·, τ ) − Txm um k∞ .

(29)

Next, by (9), Lemma 4.3, and Theorem 3.2, there exist positive
scalars τ 0 , UT∞ye and ϑ such that for τ ≥ τ 0 ,
k plim Tye (·, τ ) − Txm um k∞

(30)

N →∞

−1
≤ UT∞ye ϑ−1 k plim Tye
(·, τ ) − Tx̃−1
(·, τ )k∞ + γ̄ρ̄τ .
m um
N →∞

Finally, according to the definition of Tye (·, τ ) in (14) and
Tx̃m um (·, τ ) in (7), it follows that

C. Exact identification for acyclic latent subnetworks
Here we show that the transfer function of the manifest
subnetwork can be perfectly identified using the LSAR method
with a finite model order if the latent subnetwork is acyclic. We
start by refining the result in Proposition 4.1 and showing how,
in this case, the convergence of the LSAR matrix estimate (16)
to the optimal matrix sequence identified in Theorem 3.2 holds
in the mean-square sense.
Proposition 4.7: (The LSAR estimate converges in mean
square to optimal matrix sequence for acyclic latent
subnetworks). Consider the LTI network described by (2)
where all latent nodes are passive. Further assume that the
latent subnetwork is acyclic, i.e., there exists τ 22 ∈ Z≥1 such
that Aτ2222 = 0nl ×nl . Given the measured data sequence {y}N
1
generated from the dynamics (2) stimulated by the white
noise input {um } according to Assumption 2.3, the LSAR
estimate Âτ ({y}N
1 ) in (16) satisfies, for any τ ≥ τ 22 + 1,
∗ T
N
∗
lim E[(Âτ ({y}N
1 )− Ãτ ) (Âτ ({y}1 )− Ãτ )] = 0nm τ ×nm τ .

N →∞

τ
−1
X

Proof. If A22 is nilpotent, using Corollary 3.5, we deduce
z −i (plim Âi − Ã∗i )k∞ that the transfer function from um to ν defined in (18) is
N →∞
N →∞
i=0
Tx̃−1
T
= Inm . Consequently, the random vectors
m um xm um
−1
X
(b)
(a) τ
ν(k)’s
are
i.i.d.
with zero mean and finite second moment
k plim Âi − Ã∗i k ≤ nm βτ 2 ρ̄τ ,
(31)
≤
T
E[ν(k)ν
(k)]
=
I
nm . Define
N →∞

−1
(·, τ )k∞ = k
kplim Tye
(·, τ )−Tx̃−1
m um

i=0

where (a) holds by the sub-additivity and sub-multiplicity of
k·k and (b) follows by Proposition 4.1 and the fact that kAk ≤
nm kAkmax for any matrix A ∈ Rnm ×nm . Thus, we obtain (28)
for τ ≥ τ 0 , where β̄ , UT∞ye ϑ−1 nm β is a constant.
According to Theorem 4.4, when the length N of the measurement data is sufficiently large and the model order τ exceeds a certain threshold, the error kTye (τ ) − Txm um k∞ obtained by the LSAR method decreases exponentially with τ .
Remark 4.5: (Identification of manifest transfer function
requires higher-order models as stability margin of latent
subnetwork decreases). Even though an explicit expression of
the threshold τ 0 in Theorem 4.4 as a function of the network is
difficult to obtain, we can still make some useful observations.
From inequality (27) in the proof of Lemma 4.3, one can see
that τ 0 is an increasing function of ρ̄. Hence, as the latent
subnetwork becomes less stable (ρ(A22 ) gets closer to 1), the
corresponding τ 0 becomes larger, requiring the order of the
AR model to be higher to ensure exponential convergence. 
Remark 4.6: (Systems described by higher-order
difference equations – cont’d). As explained in Remark 3.4,
the AR representation of systems with order ν > 1 is
identical to the ν = 1 case, although they require larger AR
order τ . For large-scale systems (n  1), increasing τ
rapidly raises the number of parameters in (15), which leads
to over-parametrization of the LSAR identification. Our
simulations in Section V show how this can be overcome
both by increasing N (which is computationally costly) and
exponential regularization. Also, note that when ν > 1, the
(0)
(ν−1)
only member of the sequence of matrices A11 , . . . , A11
(denoting all current and past interactions among manifest
(0)
nodes) that is identifiable by the LSAR method is A11
(representing direct interactions among manifest states) while
the others are only identifiable in the aggregate form (5). 

1
(Âτ − Ã∗τ )ΦN ΦTN
N
(a) 1
(b) 1
=
(~ν N − ~eN )ΦTN = ~ν N ΦTN ,
N
N
where (a) follows from (15) and (19) and (b) follows from
the fact that the least-squares estimate Âτ in (16) renders
~eN ΦTN = 0nm ×nm τ . Combining the fact that the ν(k)’s are
i.i.d. and the fact that {y}k1 is a function of {ν}k−1
, we
1
deduce that ν(k) are independent of {y}k1 . This further
implies that E[ZN ] = 0nm ×nm τ . Furthermore,
ZN ,

1
E[ΦN ~ν TN ~ν N ΦTN ]
N2
1
= lim
RΦ = 0nm τ ×nm τ .
N →∞ N

T
lim E[ZN
ZN ] = lim

N →∞

N →∞

Therefore,
−1
limN →∞ E[Âτ − Ã∗τ ] = limN →∞ E[ZN ]RΦ
= 0nm ×nm τ
∗ T
and
limN →∞ E[(Âτ − Ãτ ) (Âτ − Ã∗τ )]
=
−1
−1
T
RΦ
limN →∞ E[ZN
ZN ]RΦ
= 0nm τ ×nm τ , as claimed.
We build on this result to show that the manifest transfer
function can be perfectly identified using the LSAR method
with a finite model order if the latent subnetwork is acyclic.
Theorem 4.8: (Exact manifest transfer function
identification for acyclic latent subnetworks). Under the
assumptions of Proposition 4.7, for any τ ≥ τ 22 + 1,
plim Tye ({y}N
1 , τ ) = Txm um .

N →∞

Proof.
We have plimN →∞ Âτ ({y}N
= Ã∗τ from
1 )
Proposition 4.7, which combined with (31) implies
−1
plim Tye
(τ ) = Tx̃−1
(τ ).
m um

N →∞

Moreover, from Corollary 3.5, we have Tx̃m um (τ ) = Txm um .
The statement then follows from (29) and Lemma 4.3.

9

In this section, we illustrate the performance of
least-squares auto-regressive estimation in identifying the
manifest transfer function in two examples, a deterministic
directed ring network and a group of Erdős–Rényi random
networks. We pay particular attention to the behavior
displayed as the length of measured data and the model
order change. In both examples, the input signal is a white
Gaussian process with unit variance.
Example 5.1: (Directed ring network). Consider a directed
ring network of 40 nodes with self-loops and all edge
weights equal to α = 0.25. The nodes with indices
{5, 23, 33, 34, 36} are manifest and the remaining 35 nodes
are passive latent. Fig. 2.(a) shows a 3D plot of the
identification error kTye − Txm um k∞ of the LSAR method,
with axes corresponding to length of measured data and
model order, respectively. We note that, when the measured
data length N is small, increasing the AR model order τ
does not provide better estimation of the manifest transfer
function. Similarly, when the model order τ is too low,
increasing the data length N is not helpful either. Instead,
when N and τ increase simultaneously, the LSAR method
provides good estimation of the manifest transfer function
without any knowledge of the latent nodes, as predicted by
Theorem 4.4. In Fig. 2.(b), we fix N = 106 and show that
the error of the model obtained by the LSAR method is
quite similar to the error kTx̃m um − Txm um k∞ of the ideal
AR model from Theorem 3.2. Note that the latter requires

in Â0 that are smaller than 0.1. We consider a fixed length
N = 106 of measured data and analyze the effect of varying
model order. Fig. 3 shows a 3D plot of the error in the
identification of the manifest transfer function by the LSAR
estimation, with axes corresponding to network index and
model order, respectively. One can see the improvement in

1.4
1.2
H∞ norm error

V. S IMULATIONS

1
0.8
0.6
0.4
0.2
2

2

4

6
Network index

4
6

8
10

Fig. 3. Illustration of the H∞ -norm error of the LSAR with respect to the
model order τ for the group of G(10, 0.35) Erdős–Rényi random networks
of Example 5.2. Performance improves as the model order τ increases for all
10 networks. The length of measured data is N = 106 .

performance as the model order increases for all 10
networks. Fig. 4 compares the identification error of the
LSAR method for the networks with indices 1, 6, 8, 10 in
Fig. 3 against the error of the optimal AR model from
Theorem 3.2. The latent subnetwork of network 6 is acyclic
Random network 1

Random network 6

H∞ norm error

2

H∞norm error

H∞ norm error

0.25
3
2

0.2

0.15

1

1
0.8
0.6
0.4
0.2

1.5
1
0.5

0.1

1

2

0.05

2
10

3

4

5

6

7

8

1

2

Random network 8

3

4

5

6

7

8

Random network 10

1.5
6 8

(a)

6

τ

4

2

1

2

3

4

5

6

7

8

9

10

τ

(b)

Fig. 2. H∞ -norm errors for the directed ring network of Example 5.1. (a)
The H∞ -norm error of the LSAR method as a function of data length N and
model order τ . Performance improves as N and τ increase. (b) Comparison of
the H∞ -norm errors of the LSAR method (blue dotted lines) and the optimal
AR model from Theorem 3.2 (red dashed lines) for N = 106 .

knowledge of the true adjacency matrix A, and we use it
here merely for comparison purposes.

Example 5.2: (Erdős–Rényi random network). Here we
consider a group of 10 Erdős–Rényi random networks [31].
Each network in the group is of type G(10, 0.35), with 5
manifest nodes chosen randomly and the remaining 5 nodes
are latent. Each pair of edges (i, j), (j, i), 1 ≤ i < j ≤ 10
has nonzero weights with probability 0.35 (we choose edges
in pairs so that, when plotting the graph, the edge direction
can be omitted). The weight of each edge has a uniform
distribution in {x ∈ R | 0.1 < x < 0.35} (note that (i, j) and
(j, i) can have different weights). Because of rounding errors
in the numerical computation, the estimated coefficient
matrices (16) of the AR model are usually full matrices. The
lower bound on the edge weights allows us to discard entries

H∞ norm error

4
log (N)

τ

8

1.5
1
1
0.5

0.5
1

2

3

4

5

τ

6

7

8

1

2

3

4

5

6

7

8

τ

Fig. 4. Comparison of the H∞ -norm errors of the LSAR method (red dashed
lines) and the optimal AR model from Theorem 3.2 (blue dotted lines) for
the Erdős–Rényi random networks with indices 1, 6, 8, 10 in Fig. 3. The estimation error for network 6 becomes 0 when the AR model has order higher
than 1 because the latent subnetwork is acyclic with τ 22 = 1. The length of
measured data is N = 106 .

(with A22 = 05×5 ), and the estimation error goes to 0 when
the AR model has order higher than τ 22 = 1, as predicted by
Theorem 4.8. To illustrate our observations in Remark 4.2
regarding the identification of manifest and latent
interactions, Fig. 5 shows on the left the networks with
indices 1, 6, 8, 10 of Fig. 3 and on the right the
corresponding reconstructions obtained with the LSAR
method. The indirect interactions represented by dashed
edges in these plots imply the presence of latent nodes. For
comparison, we have also used the brain connectivity
estimator technique called direct directed transfer function

10

(dDTF) measure [8], [25] from neuroscience to identify
direct connections between nodes. This technique is a
refinement of the directed transfer function (DTF) approach,
which instead cannot distinguish between direct and indirect
connections. We have employed the dynamical modeling
method within the Source Information Flow Toolbox
(SIFT) [32], [33] in EEGLAB [34], which is a widely used
open-source toolbox for EEG analysis. Fig. 6 shows the
interaction topology among the 5 manifest nodes in network
10 identified by SIFT using the dDTF measure. The dDTF
measure is in the frequency domain and can also be a
function of time (e.g., for time-varying networks). Since our
networks are time-invariant, the time axis can be ignored.
The plot shows that the dDTF identifies roughly equally
strong connections for (2, 4) (which is in reality mediated by
latent nodes) and (4, 5) (which is a true direct connection).
This is in contrast with the identification made with the
LSAR method presented in Fig. 5(d).


(a)

5
4
2

2

1

3

5
3

4
1

(b)
(a)

(b)
1

3

3

1

5

2

2
5

(c)

4

4

(d)

Fig. 5. Left: Erdős–Rényi random networks corresponding to the networks
with indices 1 (a), 6 (b), 8 (c), 10 (d) in Fig. 3, where red circles represent
manifest nodes and black circles represent latent nodes. Right: reconstructed
interaction graphs of the manifest subnetworks using the LSAR method. The
numbers next to these nodes indicate their indices. A blue solid edge represents direct interaction and a black dashed edge represents indirect interaction
through latent nodes. Note that the latent subnetwork of network 6 is acyclic.

Example 5.3: (Cortical brain network identification from
EEG data). In this example, we apply our method to a
multi-channel electroencephalogram (EEG) time-series
recorded from a human scalp during a selective visual
attention experiment in order to identify the manifest and
latent-mediated connections among the channels. The EEG
data is taken from the sample dataset available in the
EEGLAB MATLAB toolbox [34]. This dataset contains
recordings from 32 channels for more than 3 seconds with
Ts = 7.8 ms sampling time (128 Hz sampling frequency).
Channel locations are shown in Fig. 10(a) on a top (axial)
view of the skull. During the experiment, the subject is
asked to perform specific motor actions in response to
certain visual stimuli, requiring coordination among several
cortices. We take the first 13 EEG channels corresponding to
the fronto-temporal cortical areas (shown as blue squares in
Fig. 10(a)) as the manifest nodes and the remaining channels
as well as the truly hidden brain regions (the ones not
probed in the test) as the latent nodes. In the following, we
present the results of identifying the direct and indirect
connections among the manifest nodes using the LSAR

(c)
Fig. 6. (a) The interaction topology identified by the dDTF method for the
Erdős–Rényi network with index 10. (b and c) A zoom-in of the (indirect)
connection (2, 4) and the (direct) connection (4, 5), resp.

method as well as the dDTF algorithm [8], [25] and the S+L
algorithm of [22]. For each method, we only keep the edges
whose identified weights are above a certain threshold θ
(which we choose as a proportional constant α ∈ (0, 1) times
the largest edge weight in the network).
In neuroscience, the brain dynamics generating the EEG
data are usually approximated by a high-order AR model of
the form (12) (ν & 10). As mentioned in Remark 4.6, larger τ
and thus larger number of parameters are then required, which
may lead to over-parametrization. To prevent this, we use an
exponentially-regularized version of (16) by minimizing
tr(~eN ~eTN + γ Âτ P P T ÂTτ ),
−(τ −1)

(32)

where P = diag(1, ρ−1
)⊗Inm and, ideally, ρ0 =
0 , . . . , ρ0
ρ(A22 ) (in practice, it is found by trial and error). The role
of the exponential regularizer is to encourage the higher-order
AR terms to decay exponentially, as Ã∗i do. In the simulations
that follow, we have used γ = 10 and ρ0 = 0.9.
Fig. 7 shows the reconstructed manifest subnetwork with
direct and indirect connections using the LSAR method for
τ = 15 and different values of α. One can observe that the
sensitivity of the network structure to the threshold ratio α is

11

significant, showing that the majority of network links are
relatively weak with respect to the largest link (which is
usually a self-loop). This sensitivity, however, is smaller for
the indirect connections. Note that increasing α is a way of
enforcing sparsity among the manifest nodes similar (but not
equivalent) to [22]. Also, note that unlike [22], the manifest
subnetwork estimated by our method is directed (though
directions are not shown in Fig. 7 for simplicity).

FC5

F4

FC5

F4

Fz

Fz

FC1

FC1
F3

FC2

F3
FC2

FPz
FC6

FPz
FC6

T8

T8

T7

T7
Cz

FC5

F4

FC5

F4

Fz

Fz

FC1

C3

Direct
Indirect

F3

F3
FC5

FC2

Cz
C3

C4

(a) λ = 1, γλ = 0.1, l = 1

FC1

C4

(b) λ = 1, γλ = 0.5, l = 13

F4

FC2

FC5

F4

Fz

FPz

Fz

FC1

FPz

FC1
F3

FC6

Direct
Indirect

F3

FC6
FC2
T8

FC2

T8

T7

FPz

T7
Cz
C3

Cz
C3

C4

(a) α = 0.15

FC6

FPz
FC6

C4

T8

(b) α = 0.12

T7

T8
T7

Cz
FC5

F4

FC5

C3

F4

Fz

Fz

FC1

(c) λ = 1, γλ = 1, l = 18

FC1
F3

FC2

F3
FC2

FPz
FC6

C4

FPz

Cz
C3

C4

(d) λ = 0.1, γλ = 0.5, l = 50

Fig. 8. Reconstructed manifest subnetwork for Example 5.3 using the S+L
method in [22]. The direct (solid blue) and indirect (dashed black) connections are depicted for different values of weight parameters (λ, γ) and fixed
threshold ratio α = 0.01. l represents the estimated number of latent nodes.

FC6
T8
T7

T8
T7

Cz
C3

C4

(c) α = 0.09

Cz
C3

C4

(d) α = 0.06

Fig. 7. Reconstructed manifest subnetwork for the EEG data in Example 5.3
using our proposed method with the exponentially-regularized objective function (32) and γ = 10, ρ0 = 0.9, and τ = 15. The direct (solid blue) and
indirect (dashed black) connections are depicted for different values of threshold ratio α. For each value of α, the connections whose weights are smaller
than α times the largest network weight are removed.

For comparison, Fig. 8 shows the reconstructed manifest
subnetwork with direct and indirect connections using the
S+L method of [22] for n = 55 . Although the use of a
threshold value is not prescribed in [22], we have used a
fixed value of α = 0.01 for all values of (λ, γ), since the
absence of a threshold (α = 0) results in all nodes being
estimated to be (both directly and indirectly) connected. This
lack of sparsity occurs for all values of (λ, γ) (no matter
how large they are chosen), unless extremely large values are
employed, which results in a fully disconnected network.
From various plots, we see that even with the use of a
threshold value all the nodes are estimated to be indirectly
connected, with the sparsity of direct connections and the
estimated number of latent nodes being determined by
(λ, γ). This abundance of indirect connections and
parameter-based tuning of direct connectivity is similar to
our results in Fig. 7, even though the details of the
reconstructed networks do not exactly match.
5 n represents the model order in [22]. While the role of the model order
is not discussed in the reference, the use of higher-order models significantly
increases the computational cost of the algorithm. Also, note that there is no
one-to-one correspondence between the subfigures of Figs. 7-9.

Fig. 9 shows the result of applying both the Directed
Transfer Function (DTF) [35] and direct Directed Transfer
Function (dDTF) methods to the EEG channel data to
estimate the indirect and direct connections between the
manifest nodes, respectively, for different frequency bands.
Both methods are applied to the data using the EEGLAB
SIFT plugin for τ = 15 (selected based on SIFT Model
Order Selection). In all cases, a constant threshold ratio
α = 0.1 is used and the value of the threshold is computed
with respect to the largest off-diagonal link weight in the
same frequency. As can be seen, the connectivity pattern is
considerably different between lower and higher frequencies,
where several pairs are not even indirectly connected over
the δ-θ band. This is in contrast to the reconstructed
networks of Fig. 7 in which most pairs are at least indirectly
connected, even for threshold values as large as α = 0.15.
Nevertheless, a common feature of all the reconstructed
networks in Figs. 7-9 is that the density of direct connections
is higher in the fronto-central (FC) areas and lower in central
(C) areas and midline frontal pole (FPz). The independence
of this sparsity pattern from the employed reconstruction
method and parameter value suggests that it is a robust
feature of the actual brain connectivity among these areas.
Since the true network structure is unknown for this
example (and hence the methods are not directly
comparable), we validate our LSAR estimated connectivity
based on its ability to predict future (i.e., unseen) channel
activity. Thus, we used the first 80% of data for LSAR
estimation and the last 20% for evaluation, which is based on
PN 0
2
k=N +1 ke(k)k
2
,
(33)
R = 1 − PN
0
2
k=N +1 ky(k)k

12

FC5

F4

FC5

F4

Fz

Fz

FC1

Direct
Indirect

FC1
F3

FC2

F3
FC2

FPz
FC6

FPz
FC6

T8

T8

T7

T7
Cz
C3

Cz
C3

C4

(a) f = 1Hz (δ band)
FC5

C4

(b) f = 5Hz (θ band)

F4

FC5

F4

Fz

Fz

FC1

FC1
F3

FC2

F3
FC2

FPz
FC6

FPz
FC6

T8
T7

T8
T7

Cz
C3

C4

(c) f = 20Hz (β band)

Cz
C3

C4

(d) f = 50Hz (γ band)

Fig. 9. Reconstructed manifest subnetwork for Example 5.3 using the combination of DTF and dDTF estimation methods. The direct (solid blue) and
indirect (dashed black) connections are illustrated for different frequency values and fixed threshold ratio α = 0.1.

denoting the percentage of the future channel activity that is
correctly predicted by the model [26, §16.4], where
0
{y(k)}N
k=N +1 is the latter data sequence not used for
estimation. The blue curve in Fig. 10(c) shows the value of
R × 100% for the LSAR method as a function of model
order for the same selection of nodes as above (i.e.,
anterior)6 . This shows that the method is capable of
predicting more than 96.5% of unseen data with model
orders τ = 15 ∼ 20 (which is relatively low given the large
number of latent nodes and the high order of the underlying
brain dynamics). It should be noted that the R-value is not a
suitable measure for comparison among the networks
obtained by the LSAR, S+L, and dDTF methods. On the one
hand, the AR model underlying the dDTF method is almost
identical to the LSAR model used here, resulting in almost
identical R values, while the reconstructed networks are
considerably different (c.f. Figs. 7 and 9) due to different
interpretations of on the model implications for network
connectivity. On the other hand, the R value is not
well-defined for the S+L method since the right-hand side
of (33) is negative, i.e., the reconstructed AR model has
extremely poor prediction performance. This is not
surprising as the S+L method is aimed at maximizing the
entropy (and thus minimizing predictability).
Next, we analyzed the effect of the choice of manifest nodes
on the reconstructed network. In addition to selecting the 13
most anterior cortical nodes as above, we performed other
runs where we selected the 13 most posterior nodes and 13
random nodes to reconstruct the manifest network using the
LSAR method. We show in Fig. 10 these node choices (a),
6 Edge

values are not thresholded (α = 0) for computing R values.

the reconstructed network for the posterior (b) and random (d)
selections (α = 0.12), and (c) the R values for all three cases.
Interestingly, the density of direct connections is significantly
higher among the posterior nodes. Also, the LSAR prediction
performance is significantly lower in this case, suggesting less
conformity of the occipito-parietal cortex to the simplifying
assumptions of our AR model (linearity and passivity of latent
nodes). Consistently, the network density and R value of the
random case interpolates between the anterior and posterior
cases, as expected.
Finally, an interesting observation in Fig. 10(c) is that,
even an AR model with τ = 2 can predict about 95% of
unseen data in all cases. This, at first glance, questions the
need for any higher-order models as far as prediction is
concerned. Nevertheless, notice that even an AR model with
τ = 1, corresponding to an isolated manifest subnetwork,
can predict 90% of unseen data, while the visual
discrimination task performed by the subject heavily relies
on coordination between posterior (visual) and anterior
(motor planning and execution) areas. The reason why this
model can predict unseen data so well is in the strong
dominance of first-order local dynamics of every area (the
diagonal of Ã0 ) over the rest of network dynamics.7 Thus,
the prediction performance of a first-order model serves as a
baseline for higher orders, capturing the contribution of local
interactions to the overall brain dynamics. This enlightens
why the ∼ 1% improvement in prediction performance as we
go from τ = 2 to τ = 15 ∼ 20 is significant.
VI. C ONCLUSIONS
We have considered the problem of identifying the
interaction structure among a group of nodes, termed
manifest, that can be directly actuated and measured, and are
part of a larger linear-time invariant network containing an
unknown number of latent nodes. We have shown that, if
there are no inputs to the latent nodes, then the transfer
function of the manifest subnetwork can be approximated to
any degree of accuracy by means of an auto-regressive
model. We have proposed a least-squares estimation method
that uses measured data to generate estimates that converge
in probability to this AR model exponentially fast as the
length of data and the model order increase. The estimation
method does not require any knowledge of the number or the
states of the latent nodes. We have illustrated our results in a
directed ring network, a group of Erdős–Rényi random
graphs, and on a time-series of EEG data recorded from the
human brain. Future work will investigate the sensitivity of
the estimation’s performance to latent nodes, the
characterization of particular network structures which are
easier or more difficult to identify, the application of our
results to the analysis of brain data, and the extension of the
results to network models where, in addition to manifest and
latent, there are nodes that can be actuated but not measured,
and nodes that can be measured but not actuated.
7 This can be easily seen by inspecting the AR coefficients Ã estimated
i
from data, and is physiologically justified as each area is composed of millions
of neurons that are locally densely connected and serve specific purposes but
only (relatively) sparsely connected with remote areas.

13

EOG1

P8
Anterior
Posterior
Random

FP z

EOG2

P4
Pz

Direct
Indirect

PO7
P3

F3

F4

Fz

PO3
FC 5

FC 1

FC 6

FC 2

P7
T7

C3

C4

Cz

T8

POz
CP 1

CP 5

CP 2

CP 6

O2
Pz

P3

P4

P7

PO4

P8
PO 7

PO 3

PO z

PO 4

O1

Oz

O2

Oz

PO 8

PO8

(a)

O1

(b)
CP5

C4
T7

97

CP2

96

Fz

R × 100%

95

P3

94
93

F3

92

P4

91
Anterior
Posterior
Random

90
89
88
1

5

10

15

20

O2
P8
25

Oz
POz

τ

PO8

(c)
(d)
Fig. 10. Comparison between different selections of manifest nodes in Example 5.3: (a) Electrode locations. (b and d) The reconstructed network for
the 13 posterior nodes and 13 random nodes, resp. (α = 0.12). (c) Prediction
performance R for the three different choices of manifest nodes (reconstructed
network for anterior selection is given in Fig. 7(b)).

ACKNOWLEDGMENTS
The authors would like to thank Dr. John Iversen for
numerous discussions, helpful feedback on the contents of
the paper, and introducing us to the SIFT Matlab toolbox.
This work was supported by NSF Award CNS-1329619.
R EFERENCES
[1] Y. Zhao and J. Cortés, “Identification of linear networks with latent
nodes,” in American Cont. Conf., (Boston, MA), pp. 173–178, July 2016.
[2] Y. X. R. Wang and H. Huang, “Review on statistical methods for gene
network reconstruction using expression data,” Journal of Theoretical
Biology, vol. 362, pp. 53–61, 2014.
[3] A. Julius, M. Zavlanos, S. Boyd, and G. J. Pappas, “Genetic network
identification using convex programming,” IET Systems Biology, vol. 3,
no. 3, pp. 155–166, 2009.
[4] V. Sakkalis, “Review of advanced techniques for the estimation of
brain connectivity measured with eeg/meg,” Computers in Biology and
Medicine, vol. 41, no. 12, pp. 1110–1117, 2011.
[5] S. L. Bressler and A. K. Seth, “Wiener–Granger causality: a well established methodology,” Neuroimage, vol. 58, no. 2, pp. 323–329, 2011.
[6] J. R. Iversen, A. Ojeda, T. Mullen, M. Plank, J. Snider, G. Cauwenberghs, and H. Poizner, “Causal analysis of cortical networks involved
in reaching to spatial targets,” in Annual Int. Conf. of the IEEE Eng. in
Medicine and Biology Society, (Chicago, IL), pp. 4399–4402, 2014.
[7] A. Korzeniewska, C. M. Crainiceanu, R. Kuś, P. J. Franaszczuk, and
N. E. Crone, “Dynamics of event-related causality in brain electrical
activity,” Human Brain Mapping, vol. 29, no. 10, pp. 1170–1192, 2008.
[8] M. Kamiński, M. Ding, W. A. Truccolo, and S. L. Bressler, “Evaluating
causal relations in neural systems: Granger causality, directed transfer
function and statistical assessment of significance,” Biological Cybernetics, vol. 85, no. 2, pp. 145–157, 2001.
[9] C. D. Godsil and G. F. Royle, Algebraic Graph Theory, vol. 207 of
Graduate Texts in Mathematics. Springer, 2001.
[10] J. Sun, D. Taylor, and E. M. Bollt, “Causal network inference by optimal causation entropy,” SIAM Journal on Applied Dynamical Systems,
vol. 14, no. 1, pp. 73–106, 2015.
[11] M. Nabi-Abdolyousefi and M. Mesbahi, “Network identification via
node knockout,” IEEE Transactions on Automatic Control, vol. 57,
no. 12, pp. 3214–3219, 2012.

[12] S. Shahrampour and V. M. Preciado, “Topology identification of directed
dynamical networks via power spectral analysis,” IEEE Transactions on
Automatic Control, vol. 60, no. 8, pp. 2260–2265, 2015.
[13] M. Timme, “Revealing network connectivity from response dynamics,”
Physical Review Letters, vol. 98, no. 22, p. 224101, 2007.
[14] D. Materassi, G. Innocenti, L. Giarré, and M. V. Salapaka, “Model identification of a network as compressing sensing,” Systems & Control Letters, vol. 62, no. 8, pp. 664–672, 2013.
[15] D. Materassi and G. Innocenti, “Topological identification in networks of
dynamical systems,” IEEE Transactions on Automatic Control, vol. 55,
no. 8, pp. 1860–1871, 2010.
[16] M. J. Choi, V. Y. Tan, A. Anandkumar, and A. S. Willsky, “Learning
latent tree graphical models,” Journal of Machine Learning Research,
vol. 12, pp. 1771–1812, 2011.
[17] D. Materassi and M. V. Salapaka, “Network reconstruction of dynamical polytrees with unobserved nodes,” in IEEE Conf. on Decision and
Control, (Maui, Hawaii, USA), pp. 4629–4634, 2012.
[18] J. Gonçalves and S. Warnick, “Necessary and sufficient conditions for
dynamical structure reconstruction of LTI networks,” IEEE Transactions
on Automatic Control, vol. 53, no. 7, pp. 1670–1674, 2008.
[19] Y. Yuan, G. B. Stan, S. Warnick, and J. Goncalves, “Robust dynamical
network structure reconstruction,” Automatica, vol. 47, no. 6, pp. 1230–
1235, 2011. Special Issue on Systems Biology.
[20] Y. Yuan, K. Glover, and J. Goncalves, “On minimal realisations of dynamical structure functions,” Automatica, vol. 55, pp. 159–164, 2015.
[21] V. Chandrasekaran, P. A. Parrilo, and A. S. Willsky, “Latent variable
graphical model selection via convex optimization,” The Annals of Statistics, vol. 40, no. 4, pp. 2005–2013, 2012.
[22] M. Zorzi and R. Sepulchre, “AR identification of latent-variable graphical models,” IEEE Transactions on Automatic Control, vol. 61, no. 9,
pp. 2327–2340, 2016.
[23] C. I. Byrnes, S. V. Gusev, and A. Lindquist, “From finite covariance
windows to modeling filters: A convex optimization approach,” SIAM
Review, vol. 43, no. 4, pp. 645–675, 2001.
[24] E. Bullmore and O. Sporns, “Complex brain networks: graph theoretical
analysis of structural and functional systems,” Nature Reviews Neuroscience, vol. 10, no. 3, pp. 186–198, 2009.
[25] A. Korzeniewska, M. Mańczak, M. Kamiński, K. J. Blinowska, and
S. Kasicki, “Determination of information flow direction among brain
structures by a modified directed transfer function (dDTF) method,”
Journal of Neuroscience Methods, vol. 125, no. 1, pp. 195–207, 2003.
[26] L. Ljung, System Identification: Theory for the User. Prentice Hall
information and system sciences series, Prentice Hall, 1999.
[27] P. Henrici, Applied and Comp. Complex Analysis, Volume 1: Power Series Integration Conformal Mapping Location of Zero. Wiley, 1988.
[28] R. Durrett, Probability: Theory and Examples. Series in Statistical and
Probabilistic Mathematics, Cambridge University Press, 4th ed., 2010.
[29] A. Papoulis and S. U. Pillai, eds., Probability, Random Variables and
Stochastic Processes. McGraw-Hill, 2002.
[30] H. Weyl, “Das asymptotische verteilungsgesetz der eigenwerte linearer
partieller differentialgleichungen (mit einer anwendung auf die theorie
der hohlraumstrahlung),” Math. Annalen, vol. 71, pp. 441–479, 1912.
[31] B. Bollobás, Random Graphs. Cambridge University Press, 2 ed., 2001.
[32] T. Mullen, A. Delorme, C. Kothe, and S. Makeig, “An electrophysiological information flow toolbox for EEGLAB,” Biological Cybernetics,
vol. 83, pp. 35–45, 2010.
[33] A. Delorme, T. Mullen, C. Kothe, Z. A. Acar, N. Bigdely-Shamlo,
A. Vankov, and S. Makeig, “EEGLAB, SIFT, NFT, BCILAB, and ERICA: new tools for advanced EEG processing,” Computational Intelligence and Neuroscience, vol. 2011, p. 10, 2011.
[34] A. Delorme and S. Makeig, “EEGLAB: an open source toolbox for
analysis of single-trial EEG dynamics including independent component
analysis,” J. of Neuroscience Methods, vol. 134, no. 1, pp. 9–21, 2004.
[35] M. J. Kaminski and K. J. Blinowska, “A new method of the description
of the information flow in the brain structures,” Biological Cybernetics,
vol. 65, no. 3, pp. 203–210, 1991.

A PPENDIX
Lemma A.1: Given two vectors a, b ∈ Rn , it holds for any
M ∈ R>0 that kabT kmax ≤ M −1 aT a + M bT b.
Proof. By definition of the max norm,
n
X
2
2
kabT kmax = max |ai bj | ≤
(M −1 |ai | + M |bi | )
1≤i,j≤n

i=1

= M −1 aT a + M bT b.

