# Nonparametric Identification of Kronecker Networks

Mattia Zorzi<EMAIL_ADDRESS>Dipartimento di Ingegneria
dell’Informazione, Università degli studi di Padova, via Gradenigo 6/B, 35131
Padova, Italy

###### Abstract

We address the problem to estimate a dynamic network whose edges describe
Granger causality relations and whose topology has a Kronecker structure. Such
a structure arises in many real networks and allows to understand the
organization of complex networks. We proposed a kernel-based PEM method to
learn such networks. Numerical examples show the effectiveness of the proposed
method.

###### keywords:

Linear system identification; Sparsity inducing priors; Kernel-based methods;
Gaussian processes.

##  1 Introduction

In many applications we collect high dimensional data which describe systems
having a wide amount of variables. Thus, there is the need to understand how
those variables are “connected each other”, i.e. we also need to learn from
data the topology of the graphical model (Lauritzen, 1996) representing those
interactions. It is worth noting that a graphical model is meaningful only if
it is not complete, i.e. it has some degree of sparsity in terms of edges. It
is clear that the network identification paradigm depends on the definition of
“interactions”. The latter can describe conditional dependence relations
between the stochastic processes characterizing the variables of the system
and whose resulting network is undirected (Avventi et al., 2013; Zorzi &
Sepulchre, 2016; Ciccone et al., 2020; Alpago et al., 2018; Zorzi, 2019). In
the case of directed graphical models there is an edge from node $h$ to node
$j$ if we need the process at node $h$ in order to reconstruct the process at
node $j$. It is worth noting that such reconstruction can be both acausal
(Materassi & Innocenti, 2010; Materassi & Salapaka, 2020; Veedu et al., 2020;
Sepehr & Materassi, 2019; Doddi et al., 2020) and causal (Chiuso & Pillonetto,
2012; Zorzi & Chiuso, 2015). In the latter case, we say that edges describes
Granger causality relations (Granger, 1969). Finally, it is worth noting that
a closely related topic is the problem to learn modules of a network only
using data from neighbors nodes, see Ramaswamy et al. (2021); Ramaswamy &
Vandenhof (2020).

In this we paper we focus our attention on Kronecker networks (Leskovec &
Faloutsos, 2007) whose topology possesses some properties emerging in many
real networks, e.g. heavy-tailed degree distribution, small diameter, see
Leskovec et al. (2010) for more details. Such a structure is also useful to
understand the organization of complex networks (Leskovec, 2009).
Tsiligkaridis et al. (2013) consider a network identification problem whose
Kronecker structure is not only imposed on the topology of the network but
also on the parameters of the model. The latter constraint, however, could be
restrictive in some applications, e.g. for spatio-temporal MEG/EEG modeling
(Bijma et al., 2005). A possible way to overcome this restriction is to impose
the structure “sum of Kronecker products” on the parameters of the model
(Sinquin & Verhaegen, 2019; Tsiligkaridis & Hero, 2013) at the price that we
loose a meaningful definition of interaction, i.e. the resulting network is
fully connected. In Zorzi (2020a, 2018), instead, it has been considered the
problem to learn undirected dynamic networks having the Kronecker structure on
the topology, but not on the parameters of the model. Such a framework,
however, has not been extended to the case of directed dynamic networks yet.

The present paper considers the problem to learn directed dynamic networks
whose edges characterize Granger causality relations and whose topology has a
Kronecker structure.

A well established system identification paradigm is the so called prediction
error method (PEM), see Ljung (1999); Söderström & Stoica (1989). Within this
framework, candidate models are described in fixed parametric model
structures, e.g. ARMAX. Thus, the main difficulty is the correct choice of the
best model structure which is usually performed by AIC and BIC criteria
(Akaike, 1974). Kernel-based PEM methods aim to overcome this issue, see
Pillonetto & De Nicolao (2010); Chen et al. (2012); Chiuso (2016); Mu et al.
(2018); Ljung et al. (2020); Zorzi (2020b). More precisely, the candidate
model, described through the predictor impulse responses, is searched in an
infinite dimensional nonparametric model class; this is clearly an ill posed
problem because we have a finite set of measured data. However, it can be made
into a well posed one using a penalty term which favors models with specific
features. In the Bayesian viewpoint, this is equivalent to introduce an a
priori Gaussian probability (prior) on the unknown model. Hence, the prior
distribution is characterized by the covariance (i.e. kernel) function.

The contribution of this paper is the introduction of a kernel-based PEM
method to estimate Kronecker networks. The a priori information is that the
impulse responses must be Bounded Input Bounded Output (BIBO) stable and the
topology of the network respects the Kronecker structure. We derive the
corresponding kernel functions by using the maximum entropy principle.

The kernel functions are characterized by the decay rate of the predictor
impulse responses and by the number of edges in the network. These features
are not known and characterized by the so called hyperparameters. We estimate
them by minimizing the negative log-marginal likelihood of the measured data,
see Rasmussen & Williams (2006). Moreover, we show that the negative log-
marginal likelihood automatically sets some hyperparameters in such a way that
the a priori information is that the network has a sparse Kronecker structure.
The degree of sparsity depends on the variance of the noise process affecting
each node.

Finally, we provide some numerical experiments on synthetic data in order to
test the effectiveness of the proposed method. Moreover, we use the latter for
learning Granger causality relations in a bike sharing system.

The paper is organized as follows. In Section 2 we introduce the Kronecker
model corresponding to a network characterized in terms of Granger causality
relations. In Section 3 we introduce the Kronecker network identification
problem. In Section 4 we derive the maximum entropy kernels inducing BIBO
stability and a Kronecker structure for the topology of the network, while
Section 5 is devoted to the estimation of the hyperparameters. In Section 6 we
present some numerical examples. Finally, we draw the conclusions in Section
7.

### Notation

$\mathbb{N}$ is the set of natural numbers. Given a finite set $\mathcal{I}$,
$|\mathcal{I}|$ denotes its cardinality. $\mathbb{E}[\cdot]$ denotes the
expectation, while $\mathbb{E}[\cdot|\cdot]$ denotes the conditional mean.
Given three (possibly infinite dimensional) random vectors $\mathrm{a}$,
$\mathrm{b}$ and $\mathrm{c}$ we say that $\mathrm{a}$ is conditionally
independent of $\mathrm{b}$ given $\mathrm{c}$ if
$\mathbb{E}[\mathrm{a}|\mathrm{b},\mathrm{c}]=\mathbb{E}[\mathrm{a}|\mathrm{c}]$.
Given $G\in\mathbb{R}^{n\times p}$, $[G]_{i,j}$ denotes the entry of $G$ in
position $(i,j)$; $G>0$ ($G\geq 0$) means that $G$ is a positive
(semi-)definite matrix. $\mathcal{D}_{p}$ denotes the vector space of diagonal
matrices of dimension $p$. $\ell_{2}(\mathbb{N})$ denotes the space of
$\mathbb{R}$-valued infinite length sequences, which we think as infinite
dimensional column vectors $g:=[g_{1}\;g_{2}\;\ldots\;g_{j}\;\ldots]^{\top}$,
$g_{k}\in\mathbb{R}$, $k\in\mathbb{N}$, such that
$\|g\|_{2}:=\sqrt{\sum_{k=1}^{\infty}|g_{k}|^{2}}<\infty$. $\ell_{2}^{p\times
n}(\mathbb{N})$ is the space of matrices of sequences in
$\ell_{2}(\mathbb{N})$

$\displaystyle\Phi=\left[\begin{array}[]{ccc}(\phi^{[11]})^{\top}&\ldots&(\phi^{[1n]})^{\top}\\\
\vdots&\ddots&\vdots\\\
(\phi^{[p1]})^{\top}&\ldots&(\phi^{[pn]})^{\top}\end{array}\right]$ (4)

where $\phi^{[ij]}\in\ell_{2}(\mathbb{N})$, $i=1\ldots p$ and $j=1\ldots n$.
$\ell_{1}(\mathbb{N})$ denotes the space of $\mathbb{R}$-valued infinite
length sequences $g$ such that $\|g\|_{1}:=\sum_{k=1}^{\infty}|g_{k}|<\infty$.
$\ell_{1}^{p\times n}(\mathbb{N})$ is defined in similar way.
$\mathcal{S}_{2}(\mathbb{N})$ denotes the space of symmetric infinite
dimensional matrices $K$ such that
$\|K\|_{2}:=\sqrt{\sum_{i,j=1}^{\infty}|[K]_{i,j}|^{2}}<\infty$.
$\mathcal{S}_{2}^{p}(\mathbb{N})$ is the space of symmetric infinite
dimensional matrices

$\displaystyle
K=\left[\begin{array}[]{cccc}K^{[11]}&K^{[12]}&\ldots&K^{[1p]}\\\
K^{[12]}&K^{[22]}&&K^{[2p]}\\\ \vdots&&\ddots&\vdots\\\
K^{[1p]}&K^{[2p]}&\ldots&K^{[pp]}\end{array}\right]$ (9)

where $K^{[ij]}\in\mathcal{S}_{2}(\mathbb{N})$, $i,j=1\ldots p$. Given
$\Phi\in\ell_{2}^{p\times n}(\mathbb{N})$, $\Psi\in\ell_{2}^{m\times
n}(\mathbb{N})$ and $K\in\mathcal{S}_{2}^{n}(\mathbb{N})$, the products
$\Phi\Psi^{\top}$ and $\Phi K\Psi^{\top}$ are understood as $p\times m$
matrices whose entries are limits of infinite sequences (Jorgesen et al.,
2011). Given $g\in\ell_{2}^{n\times 1}(\mathbb{N})$ and
$K\in\mathcal{S}_{2}^{n}(\mathbb{N})$, $\|g\|^{2}_{K^{-1}}:=g^{\top}K^{-1}g$.
The definition is similar in the case that $g$ and $K$ have finite dimension.
With some abuse of notation the symbol $z$ will denote both the complex
variable as well as the shift operator $z^{-1}y(t):=y(t-1)$. Given a
stochastic process $y=\\{y(t)\\}_{t\in\mathbb{Z}}$, with some abuse of
notation, $y(t)$ will both denote a random vector and its sample value. From
now on the time $t$ will denote _present_ and we shall talk about _past_ and
_future_ with respect to time $t$. With this convention in mind,

$\displaystyle\mathrm{y}^{-}=\left[\begin{array}[]{ccc}y(t-1)^{\top}&y(t-2)^{\top}&\ldots\end{array}\right]^{\top}$
(11)

denotes the (infinite length) past data vector of $y$ at time $t$. We denote
as $\mathrm{supp}(G(z))$ the support of the transfer matrix $G(z)$, i.e. the
entries of $G(z)$ different from the null function correspond to entries equal
to one in $\mathrm{supp}(G(z))$, otherwise the latter are equal to zero.

##  2 Kronecker Models

We consider the following nonparametric model

$\displaystyle y(t)=$ $\displaystyle G(z)y(t)+F(z)u(t)+e(t)$ (12)

where $y(t)\in\mathbb{R}^{p_{1}p_{2}}$ is the output, $u(t)\in\mathbb{R}^{m}$
is the input, $e(t)$ is zero mean normalized white Gaussian noise (WGN) whose
components have variance $\sigma_{hk}^{2}$ with $h=1\ldots p_{1}$ and
$k=1\ldots p_{2}$; $G(z)=\sum_{t=1}^{\infty}G_{t}z^{-t}$ is a BIBO stable
transfer matrix of dimension $p_{1}p_{2}\times p_{1}p_{2}$, while
$F(z)=\sum_{t=1}^{\infty}F_{t}z^{-t}$ is a BIBO stable transfer matrix of
dimension $p_{1}p_{2}\times m$. The minimum variance one-step ahead predictor
of $y(t)$ based on the past data $\mathrm{y}^{-}$, denoted by
$\hat{y}(t|t-1)$, is given by

$\displaystyle\hat{y}(t|t-1)$ $\displaystyle=G(z)y(t)+F(z)u(t).$

Thus, $e(t)$ is the one-step-ahead prediction error

$e(t)=y(t)-\hat{y}(t|t-1).$

We assume that

$\displaystyle\begin{aligned} \mathrm{supp}(G(z))&=E_{1}\otimes E_{2}\\\
\mathrm{supp}(F(z))&=A_{1}\otimes A_{2}\end{aligned}$ (13)

where $E_{1}$, $E_{2}$, $A_{1}$ and $A_{2}$ denote the adjacency matrices of
dimension $p_{1}\times p_{1}$, $p_{2}\times p_{2}$, $p_{1}\times 1$ and
$p_{2}\times m$, respectively. Let $y_{hk}$ denote the scalar process of $y$
in position $(h-1)p_{2}+k$ with $h\in I_{1}:=\\{1\ldots p_{1}\\}$ and $k\in
I_{2}:=\\{1\ldots p_{2}\\}$. In a similar way, $\mathrm{y}^{-}_{hk}$ denotes
the past of $y_{hk}$ at time $t$. Finally, $u_{i}$ denotes the scalar
component of $u$ in position $i\in I_{u}:=\\{1\ldots m\\}$ and
$\mathrm{u}_{i}^{-}$ the $i$-th component of $\mathrm{u}^{-}$. Throughout the
paper we assume that $E_{1}$, $E_{2}$, $A_{1}$ and $A_{2}$ are different from
the null matrix.

It is possible to describe the structure of model (12)-(13) using a Bayesian
network (Lauritzen, 1996). The nodes correspond to the scalar processes
$y_{hk}$’s, with $h\in I_{1}$ and $k\in I_{2}$, while $u_{i}$’s, with $i\in
I_{u}$, represent exogenous variables affecting some nodes of the network.
Then, the connections among the nodes obey the following rules:

  * •

there is a directed link from node $y_{jl}$ to node $y_{hk}$ if

$\displaystyle\mathbb{E}[y_{hk}(t)|\mathrm{y}^{-},\mathrm{u}^{-}]\neq\mathbb{E}[y_{hk}(t)|\mathrm{y}_{\tilde{j}\tilde{l}}^{-},\tilde{j}\neq
j,\tilde{l}\neq l,\mathrm{u}^{-}],$

i.e. if $\mathrm{y}_{jl}^{-}$ is needed to predict $y_{hk}(t)$ given
$\mathrm{u}^{-}$ and $\mathrm{y}_{\tilde{j}\tilde{l}}^{-}$, for any
$\tilde{j}\in I_{1}\setminus\\{j\\}$ and $\tilde{l}\in I_{2}\setminus\\{l\\}$;
In this case, we shall say $y_{jl}$ conditionally Granger causes $y_{hk}$.

  * •

there is a directed link from $u_{i}$ to node $y_{hk}$ if

$\displaystyle\mathbb{E}[y_{hk}(t)|\mathrm{y}^{-},\mathrm{u}^{-}]\neq\mathbb{E}[y_{hk}(t)|\mathrm{y}^{-},\mathrm{u}^{-}_{\tilde{i}},\tilde{i}\neq
i]$

i.e. if $\mathrm{u}_{i}^{-}$ is needed to predict $y_{hk}(t)$ given
$\mathrm{y}^{-}$ and $\mathrm{u}_{\tilde{i}}^{-}$, for any $\tilde{i}\in
I_{u}\setminus\\{i\\}$; In this case, we shall say $u_{i}$ conditionally
Granger causes $y_{hk}$.

In such a network we can recognize $p_{1}$ modules containing $p_{2}$ nodes
and sharing the same topology described by $E_{2}$, while the interaction
among these $p_{1}$ modules is described by $E_{1}$. Finally, $A_{1}$
describes the topology among the input and the modules; $A_{2}$ describes the
topology among the input components and the nodes in a module which is the
same in each module. These facts are formalized by the next proposition.

######  Proposition 1.

Consider model (12)-(13). Let $y^{\star}_{h}$ be the process obtained by
stacking $y_{hk}$ with $k\in I_{2}$ and $y^{\dagger}_{k}$ be the process
obtained by stacking $h\in I_{1}$. Then:

  * •

$y_{j}^{\star}$ conditionally Granger causes $y^{\star}_{h}$ if and only if
$[E_{1}]_{h,j}=1$;

  * •

$y_{l}^{\dagger}$ conditionally Granger causes $y^{\dagger}_{k}$ if and only
if $[E_{2}]_{k,l}=1$;

  * •

$u$ conditionally Granger causes $y^{\star}_{h}$ if and only if
$[A_{1}]_{h}=1$;

  * •

$u_{i}$ conditionally Granger causes $y^{\dagger}_{k}$ if and only if
$[A_{2}]_{k,i}=1$.

######  Proof 2.1.

We only prove the first claim, the proof of the others is similar. Let
$g^{[hk,jl]}\in\ell_{1}(\mathbb{N})$ denote the impulse response in $G(z)$
having input $y_{jl}$ and output $y_{hk}$. By (13), we have

$\displaystyle g^{[hk,jl]}\neq 0,\,\hbox{ for some }k,l\in
I_{2}\iff[E_{1}]_{h,j}=1.$ (14)

From (12) we have that

$\displaystyle y_{hk}$ $\displaystyle(t)=\sum_{\begin{subarray}{c}j\in
I_{1}\\\ l\in I_{2}\end{subarray}}[G(z)]_{(h-1)p_{2}+k,(j-1)p_{2}+l}y_{jl}(t)$
$\displaystyle+\sum_{i\in I_{u}}[F(z)]_{(h-1)p_{2}+k,i}u_{i}(t)+e_{hk}(t).$

Thus, condition (14) is equivalent to

$\displaystyle\mathbb{E}[y_{hk}(t)|\mathrm{y}^{-},\mathrm{u}^{-}]\neq\mathbb{E}[y_{hk}(t)|\mathrm{y}_{\tilde{j}\tilde{l}}^{-},\tilde{j}\neq
j,\tilde{l}\neq l,\mathrm{u}^{-}],$ $\displaystyle\hbox{for some }k,l\in
I_{2}\iff[E_{1}]_{h,j}=1.$ (15)

By (2.1), we have

$\displaystyle\mathbb{E}[y_{hk}(t)|\mathrm{u}^{-},\mathrm{y}^{-}]\neq\mathbb{E}[y_{hk}(t)|\mathrm{y}^{\star-}_{\tilde{j}},\tilde{j}\neq
j,\mathrm{u}^{-}],$ $\displaystyle\hbox{for some }k\in
I_{2}\iff[E_{1}]_{h,j}=1.$ (16)

Staking condition (2.1) for any $k\in I_{2}$, we obtain

$\displaystyle\mathbb{E}[y_{h}^{\star}(t)|\mathrm{u}^{-},\mathrm{y}^{-}]\neq\mathbb{E}[y_{h}^{\star}(t)|\mathrm{u}^{-},\mathrm{y}^{\star-}_{\tilde{j}},\tilde{j}\neq
j]$ $\displaystyle\hskip 113.81102pt\iff[E_{1}]_{h,j}=1$

which proves the claim. ∎

Figure 1: A example of dynamic Kronecker network with $p_{1}=3$ modules (red,
blue and green) composed by $p_{2}=4$ nodes and $m=1$ exogenous component.

It is worth noting that condition (13) is weaker than $G(z)$ and $F(z)$ admit
a Kronecker decomposition, e.g. $G(z)=G_{1}(z)\otimes G_{2}(z)$. Therefore, we
do not constrain the dynamics in each module to be the same. An example of
Kronecker network is depicted in Figure 1. Here, we have

$y=[\,y_{11}\,y_{12}\,y_{13}\,y_{14}\,y_{21}\,y_{22}\,y_{23}\,y_{24}\,y_{31}\,y_{32}\,y_{33}\,y_{34}\,]^{\top}$

and the adjacency matrices are:

$\displaystyle E_{1}$ $\displaystyle=\left[\begin{array}[]{ccc}0&0&0\\\
1&0&0\\\ 0&1&0\end{array}\right],\,E_{2}=\left[\begin{array}[]{cccc}1&1&0&0\\\
0&0&0&1\\\ 0&1&0&0\\\
0&0&0&1\end{array}\right],\,A_{1}=\left[\begin{array}[]{c}1\\\ 1\\\
1\end{array}\right],\,A_{2}=\left[\begin{array}[]{c}0\\\ 1\\\ 0\\\
0\end{array}\right].$ (31)

###  2.1 Motivating examples

Dynamic spatio-temporal modeling. Consider a nonparametric time-varying model
of the form

$\displaystyle
x(s)=\tilde{G}_{s}(q)x(s)+\tilde{F}_{s}(q)w(s)+v(s),\;\;s\in\mathbb{Z}$ (32)

where $x(s)\in\mathbb{R}^{p2}$ is the output, $w(s)\in\mathbb{R}^{\tilde{m}}$
is the input, $v(t)$ is normalized WGN and
$\\{\tilde{G}_{s}(q),\,s\in\mathbb{Z}\\}$,
$\\{\tilde{F}_{s}(q),\,s\in\mathbb{Z}\\}$ are sequences of BIBO stable
(strictly causal) transfer matrices of dimension $p_{2}\times p_{2}$ and
$p_{2}\times\tilde{m}$, respectively, such that

$\displaystyle\tilde{G}_{s}(q)=\tilde{G}_{s+p_{1}}(q),\;\;\tilde{F}_{s}(q)=\tilde{F}_{s+p_{1}}(q)$
(33)

where $q^{-1}$ is the shift operator corresponding to time variable $s$, e.g.
$q^{-1}x(s)=x(s-1)$. Then, if we define

$\displaystyle y(t)$
$\displaystyle=[\,x((t-1)p_{1}+1)^{\top}\,\ldots\;x(tp_{1})^{\top}\,]^{\top}\in\mathbb{R}^{p_{1}p_{2}}$
$\displaystyle u(t)$
$\displaystyle=[\,w((t-1)p_{1}+1)^{\top}\,\ldots\;w(tp_{1})^{\top}\,]^{\top}\in\mathbb{R}^{p_{1}\tilde{m}}$
$\displaystyle e(t)$
$\displaystyle=[\,v((t-1)p_{1}+1)^{\top}\,\ldots\;v(tp_{1})^{\top}\,]^{\top}\in\mathbb{R}^{p_{1}p_{2}},$

we can write (32) as (12) where $G(z)$ and $F(z)$ are BIBO stable transfer
matrices of dimension $p_{1}p_{2}\times p_{1}p_{2}$ and $p_{1}p_{2}\times
p_{1}\tilde{m}$, respectively. Moreover,

$\displaystyle
G(z)=\left[\begin{array}[]{ccccc}G_{1,1}(z)&G_{1,2}(z)&\ldots&\ldots&G_{1,p_{1}}(z)\\\
zG_{2,1}(z)&G_{2,2}(z)&\ddots&&\vdots\\\ \vdots&\ddots&\ddots&\ddots&\vdots\\\
\vdots&&\ddots&\ddots&G_{p_{1}-1,p_{1}}(z)\\\
zG_{p_{1},1}(z)&\ldots&\ldots&zG_{p_{1},p_{1}-1}(z)&G_{p_{1},p_{1}}(z)\end{array}\right]$
(39)

where $G_{hj}(z)=\sum_{t=1}^{\infty}G_{hj,t}z^{-t}$ is a transfer matrix of
dimension $p_{2}\times p_{2}$ and $z^{-1}$ is the shift operator corresponding
to time variable $t$; $F(z)$ is defined in a similar way: it is composed by
the transfer matrices $F_{hj}(z)=\sum_{t=1}^{\infty}F_{hj,t}z^{-t}$,
$h,j=1\ldots p_{1}$, of dimension $p_{2}\times\tilde{m}$ and the blocks in the
strictly lower block triangular part are multiplied by $z$. In plain words, we
have written (32) as a time-invariant model through $y(t)$ and $u(t)$. It is
interesting to note that the time variables $s$ and $t$ play a different role;
indeed, there is a decimation relationship between them and the decimation
factor is $p_{1}$. If $G(z)$ and $F(z)$ satisfy condition (13), then by
Proposition 1 we have that $E_{2}$ is the matrix describing the conditional
Granger causality relations among the components of $x(s)$, while $A_{2}$ is
the matrix describing the conditional Granger causality relations among the
components of $w(s)$ and $x(s)$. Moreover, $E_{1}$ describes the conditional
Granger causality relations of $x(s)$ over the period $1\ldots p_{1}$, while
$A_{1}$ describes the conditional Granger causality relations among $x(s)$ and
$w(s)$ over the period $1\ldots p_{1}$. We conclude that (12) can be
understood as a dynamic spatio-temporal model. In particular, if we take

$\displaystyle G(z)=\bar{G},\;\;F(z)=\bar{F},$ (40)

i.e. constant matrices, then (12) describes a spatio-temporal model without
dynamic. The latter has been used to model magnetoencephalography (MEG)
measurements for mapping brain activity (Bijma et al., 2005):
$z(k,t)=[\,x_{k}((t-1)p_{1}+1)\,\ldots\,x_{k}(tp_{1})\,]^{\top}$ models the
measurements corresponding to the $k$-th brain region during the $t$-th trial
of length $p_{1}$, while
$s(k,t)=[\,w_{k}((t-1)p_{1}+1)\,\ldots\,w_{k}(tp_{1})\,]^{\top}$ is the
stimulus applied to the $k$-th region during the $t$-th trial. These trials
are assumed independent. Moreover, they are identically distributed because in
each trial the task performed by the patient is the same and the way that the
patient reacts to the stimulus is assumed to be the same. In our framework we
do not impose the restrictive condition (40) and thus we let dependence among
the trials. This means that the proposed model takes into account the fact
that the trials are taken in a sequential way, and thus the fact that the
behaviors of the patient in the next trial can be influenced by the previous
one. Finally, such a dynamic spatio-temporal model could be also used for
modeling the conditional Granger causality relations among the number rental
bikes in a bikesharing system and the corresponding weather information, see
Section 6.2.

Multi-task learning. Consider a network composed by $p_{2}\in\mathbb{N}$
agents (i.e. each node correspond to an agent). One agent is described by a
model having an input and output. Our aim is to model such a network under
$p_{1}\in\mathbb{N}$ heterogeneous conditions (i.e. tasks). The $k$-th agent
under the $h$-th task has input $u(t)$ and output $y_{hk(t)}$. Then, we could
model the network through the $p_{1}$ models:

$\displaystyle
y_{h}^{\star}(t)=G_{h}(z)y_{h}^{\star}+F_{h}(z)u(t)+e_{h}(t),\;\;h\in I_{1}$
(41)

where $y_{h}^{\star}(t)$ has been defined in Proposition 1. A more flexible
approach is to model all the models in (41) together in order to exploit
commonalities and differences across the tasks, see Yu et al. (2009). More
precisely, we consider the model

$\displaystyle y(t)=G(z)y(t)+F(z)u(t)$

where $y(t)$ is obtained by stacking $y_{h}^{\star}(t)$ with $h\in I_{1}$.
Under the assumption that (13) holds, and in view of Proposition 1, we have
that $E_{1}$ describes the conditional Granger causality relations among the
tasks; $E_{2}$ describes the conditional Granger causality relations among the
agents in each task; $A_{1}$ describes the conditional Granger causality
relations from the input to the tasks; $A_{2}$ describes the conditional
Granger causality relations from the input components to the agents in each
task.

##  3 Problem Formulation

Assume to measure the data $\\{y(t),u(t)\\}_{t=0,..,N}$ generated by (12). In
this section, we address the problem of estimating $G(z)$ and $F(z)$ which
admit the Kronecker product decomposition in (13).

The transfer matrices $G(z)$, $F(z)$ are parametrized in terms of their
impulse response coefficients $G_{t}$ and $F_{t}$. In particular, defining
$g^{[hk,jl]}\in\ell_{1}(\mathbb{N})$ to be the impulse response from input
$(j-1)p_{2}+l$ to output $(h-1)p_{2}+k$, we have:

$\displaystyle g^{[hk,jl]}:=$
$\displaystyle\left[\begin{array}[]{cc}[G_{1}]_{(h-1)p_{2}+k,(j+1)p_{2}+l}&\ldots\end{array}\right.$
(43) $\displaystyle\hskip
5.69046pt\left.\begin{array}[]{ccc}\ldots&[G_{s}]_{(h-1)p_{2}+k,(j+1)p_{2}+l}&\ldots\\\
\end{array}\right]^{\top};$ (45)

defining $f^{[hk,i]}\in\ell_{1}(\mathbb{N})$ to be the impulse response from
input $u_{i}$ to output $(h-1)p_{2}+k$, we have

$\displaystyle
f^{[hk,i]}:=\left[\begin{array}[]{cccc}[F_{1}]_{(h-1)p_{2}+k,i}&\ldots&[F_{s}]_{(h-1)p_{2}+k,i,i}&\ldots\\\
\end{array}\right]^{\top}.$ (47)

The coefficient vectors $\theta_{g}^{\top}\in\ell_{1}^{1\times
p_{1}^{2}p_{2}^{2}}(\mathbb{N})$ and $\theta_{f}^{\top}\in\ell_{1}^{1\times
p_{1}p_{2}m}(\mathbb{N})$ are defined as follows:

$\displaystyle\theta_{g}^{\top}=$
$\displaystyle\left[\begin{array}[]{ccc|c}(\theta_{g}^{[11]})^{\top}&\ldots&(\theta_{g}^{[1p_{2}]})^{\top}&\ldots\\\
\end{array}\right.$ (49) $\displaystyle\hskip
5.69046pt\left.\begin{array}[]{c|ccc}\ldots&(\theta_{g}^{[p_{1}1]})^{\top}&\ldots&(\theta_{g}^{[p_{1}p_{2}]})^{\top}\\\
\end{array}\right]$ (51) $\displaystyle\theta_{f}^{\top}=$
$\displaystyle\left[\begin{array}[]{ccc|c}(\theta_{f}^{[11]})^{\top}&\ldots&(\theta_{f}^{[1p_{2}]})^{\top}&\ldots\\\
\end{array}\right.$ (53) $\displaystyle\hskip
5.69046pt\left.\begin{array}[]{c|ccc}\ldots&(\theta_{f}^{[p_{1}1]})^{\top}&\ldots&(\theta_{f}^{[p_{1}p_{2}]})^{\top}\\\
\end{array}\right]$ (55)

where

$\displaystyle\theta_{g}^{[hk]^{\top}}=$
$\displaystyle\left[\begin{array}[]{ccc|c}(g^{[hk,11]})^{\top}&\ldots&(g^{[hk,1p_{2}]})^{\top}&\ldots\\\
\end{array}\right.$ (57) $\displaystyle\hskip
5.69046pt\left.\begin{array}[]{c|ccc}\ldots&(g^{[kh,p_{1}1]})^{\top}&\ldots&(g^{[hk,p_{1}p_{2}]})^{\top}\\\
\end{array}\right]$ (59) $\displaystyle\theta_{f}^{[hk]^{\top}}=$
$\displaystyle\left[\begin{array}[]{ccc}(f^{[hk,1]})^{\top}&\ldots&(f^{[hk,m]})^{\top}\\\
\end{array}\right].$ (61)

The measured data $y(1)\ldots y(N)$ are stacked in the vector $\mathrm{y}^{+}$
as follows

$\displaystyle\mathrm{y}^{+}=$
$\displaystyle\left[\begin{array}[]{ccc|c|ccc}\mathrm{y}_{11}^{+\top}&\ldots&\mathrm{y}_{1p_{2}}^{+\top}&\ldots&\mathrm{y}_{p_{1}1}^{+\top}&\ldots&\mathrm{y}_{p_{1}p_{2}}^{+\top}\\\
\end{array}\right]^{\top}$ (63)

where

$\displaystyle\mathrm{y}^{+}_{hk}=$
$\displaystyle\left[\begin{array}[]{ccc}y_{hk}(N)^{\top}&\ldots&y_{hk}(1)^{\top}\\\
\end{array}\right].$ (65)

The vector $\mathrm{e}^{+}$ is defined analogously. Let us also introduce the
Toeplitz matrices $\phi_{hk}\in\ell_{2}^{N\times 1}(\mathbb{N})$,
$\psi_{i}\in\ell_{2}^{N\times 1}(\mathbb{N})$:

$\displaystyle[\phi_{hk}]_{sn}:=y_{hk}(N-s-n+1)$
$\displaystyle[\psi_{i}]_{sn}:=u_{i}(N-s-n+1)$

with $s=1\ldots N$ and $n\in\mathbb{N}$. Then, we define the regression
matrices $\Phi\in\ell_{2}^{p_{1}p_{2}N\times p_{1}^{2}p_{2}^{2}}(\mathbb{N})$
and $\Psi\in\ell_{2}^{p_{1}p_{2}N\times p_{1}p_{2}m}(\mathbb{N})$ as:

$\displaystyle\Phi$
$\displaystyle=I_{p_{1}p_{2}}\otimes\left[\begin{array}[]{ccccccc}\phi_{11}&\ldots&\phi_{1p_{2}}&\ldots&\phi_{p_{1}1}&\ldots&\phi_{p_{1}p_{2}}\\\
\end{array}\right]$ (67) $\displaystyle\Psi$
$\displaystyle=I_{p_{1}p_{2}}\otimes\left[\begin{array}[]{ccccccc}\psi_{1}&\ldots&\psi_{m}\\\
\end{array}\right]$ (69)

so that, from (12) the vector $\mathrm{y}^{+}$ containing the measured output
data satisfies the linear regression model

$\displaystyle\mathrm{y}^{+}=\Phi\theta_{g}+\Psi\theta_{f}+\mathrm{e}^{+}$
(70)

and $\hat{\mathrm{y}}^{+}:=\Phi\theta_{g}+\Psi\theta_{f}$ is the one-step
ahead predictor of $\mathrm{y}^{+}$. It is worth noting that to construct
$\Phi$ and $\Psi$ we need the remote past of the output and the input which is
not available. Thus, model (70) has to be approximated truncating $\Phi$ and
$\Psi$ (and thus $\theta_{g}$ and $\theta_{f}$). This is equivalent to impose
zero initial conditions, which is a reasonable approximation given the decay,
as a function of $t$, of the coefficients $G_{t}$, $F_{t}$ by the BIBO
stability conditions.

######  Remark 2.

In the case that we consider the dynamic spatio-temporal model of Section 2.1,
the regression matrices in (70) are different because the strictly lower block
triangular parts of $G(z)$ and $F(z)$ do not contain a delay. More precisely,
$\Phi$ and $\Psi$ are $p_{1}\times p_{1}$ block diagonal matrices where the
block in position $(h,h)$, with $h\in I_{1}$, is, respectively,

$\displaystyle
I_{p_{2}}\otimes[\,\tilde{\phi}_{11}\,\ldots\,\tilde{\phi}_{h-1,p_{2}}\,\phi_{h,1}\,\ldots\phi_{p_{1},p_{2}}\,]$
$\displaystyle
I_{p_{2}}\otimes[\,\tilde{\psi}_{1}\,\ldots\,\tilde{\psi}_{(h-1)\tilde{m}}\,\psi_{(h-1)\tilde{m}+1}\,\ldots\psi_{p_{1}\tilde{m}}\,]$

and the Toeplitz matrices $\tilde{\phi}_{hk}\in\ell_{2}^{N\times
1}(\mathbb{N})$, $\tilde{\psi}_{i}\in\ell_{2}^{N\times 1}(\mathbb{N})$ are:

$\displaystyle[\tilde{\phi}_{hk}]_{sn}:=y_{hk}(N-s-n+2)$
$\displaystyle[\tilde{\psi}_{i}]_{sn}:=u_{i}(N-s-n+2).$

Therefore, our Kronecker identification problem can be formulated in terms of
PEM as follows.

######  Problem 3.

Assume to collect the measurements $\\{y(t),u(t)\\}_{t=0\ldots N}$ and that
the dimensions $p_{1}$ and $p_{2}$ are known. Find
$\theta_{g}^{\top}\in\ell_{1}^{1\times p_{1}^{2}p_{2}^{2}}(\mathbb{N})$ and
$\theta_{f}^{\top}\in\ell_{1}^{1\times p_{1}p_{2}m}(\mathbb{N})$ corresponding
to a Kronecker model minimizing the prediction error squared norm
$\|\mathrm{y}^{+}-\Phi\theta\|^{2}_{\Sigma^{-1}\otimes I_{N}}$ with
$\Sigma=\mathrm{diag}(\sigma_{11}^{2}\ldots\sigma_{1p_{2}}^{2}\ldots\sigma_{p_{1}1}^{2}\ldots\sigma_{p_{1}p_{2}})$.

Following the nonparametric Gaussian regression approach in Pillonetto et al.
(2011), we model $\theta_{g}$ and $\theta_{f}$ as zero-mean processes with
kernels $K_{g}\in\mathcal{S}_{2}^{p_{1}^{2}p_{2}^{2}}(\mathbb{N})$ and
$K_{f}\in\mathcal{S}_{2}^{p_{1}p_{2}m}(\mathbb{N})$, respectively. These
kernels may depend upon some tuning parameters, usually called hyperparameters
and denoted with $\zeta$ hereafter. As illustrated in Section 4, according to
the maximum entropy principle $\theta_{g}$ and $\theta_{f}$ will be modeled as
Gaussian and independent. In the following $\mathcal{H}_{K}$ denotes the
reproducing Hilbert space (Aronszajn, 1950) of deterministic functions on
$\mathbb{N}$, associated with $K$ and with norm denoted by
$\|\cdot\|_{K^{-1}}$. We assume that the past data $\mathrm{y}^{-}$ neither
affects the a priori probability on $\theta_{g}$ and $\theta_{f}$ nor carries
information on $\zeta$ and $\Sigma$ (Pillonetto & De Nicolao, 2010), that is

$\displaystyle\mathbf{p}$
$\displaystyle(\mathrm{y}^{+},\theta_{g},\theta_{f},\mathrm{y}^{-}|\zeta,\Sigma)$
$\displaystyle=\mathbf{p}(\mathrm{y}^{+}|\theta_{g},\theta_{f},\mathrm{y}^{-},\zeta,\Sigma)\mathbf{p}(\theta_{g},\theta_{f}|\mathrm{y}^{-},\zeta,\Sigma)\mathbf{p}(\mathrm{y}^{-}|\zeta,\Sigma)$
$\displaystyle\approx\mathbf{p}(\mathrm{y}^{+}|\theta_{g},\theta_{f},\mathrm{y}^{-},\zeta,\Sigma)\mathbf{p}(\theta_{g},\theta_{f}|\zeta,\Sigma)\mathbf{p}(\mathrm{y}^{-}).$
(71)

Let $\hat{\theta}_{g}=\mathbb{E}[\theta_{g}|\mathrm{y}^{+},\zeta,\Sigma]$ and
$\hat{\theta}_{f}=\mathbb{E}[\theta_{f}|\mathrm{y}^{+},\zeta,\Sigma]$ be the
minimum variance estimator, respectively, of $\theta_{g}$ and $\theta_{f}$
given $\mathrm{y}^{+}$, $\zeta$ and $\Sigma$. It is well known that
$\hat{\theta}_{g}$ and $\hat{\theta}_{f}$ are almost surely solution to the
following Tikhonov-type variational problem

$\displaystyle\underset{\begin{subarray}{c}\theta_{g}\in\mathcal{H}_{K_{g}}\\\
\theta_{f}\in\mathcal{H}_{K_{f}}\end{subarray}}{\arg\min}\|\mathrm{y}^{+}-\Phi\theta_{g}-\Psi\theta_{f}\|^{2}_{\Sigma^{-1}\otimes
I_{N}}+\|\theta_{g}\|^{2}_{K_{g}^{-1}}+\|\theta_{f}\|^{2}_{K_{f}^{-1}}.$

Moreover, almost surely:

$\displaystyle\begin{aligned} \hat{\theta}_{g}&=K_{g}\Phi^{\top}(\Phi
K_{g}\Phi^{\top}+\Psi K_{f}\Psi^{\top}+\Sigma\otimes I)^{-1}\mathrm{y}^{+}\\\
\hat{\theta}_{f}&=K_{f}\Psi^{\top}(\Phi K_{g}\Phi^{\top}+\Psi
K_{f}\Psi^{\top}+\Sigma\otimes I)^{-1}\mathrm{y}^{+}\end{aligned}.$

In what follows, $\hat{G}(z)$ and $\hat{F}(z)$ denote the transfer matrices
corresponding to $\hat{\theta}_{g}$ and $\hat{\theta}_{f}$, respectively. The
main task now is to design the kernels $K_{g}$ and $K_{f}$ in such a way that
$\hat{G}(z)$ and $\hat{F}(z)$ are almost surely BIBO stable while favoring a
sparse Kronecker decomposition of $\mathrm{sup}(\hat{G}(z))$ and
$\mathrm{sup}(\hat{F}(z))$.

##  4 Maximum entropy priors

The probability law for the joint process
$\theta=[\,\theta_{g}^{\top}\;\theta_{f}^{\top}]^{\top}$ under desired
constraints can be obtained by the maximum entropy principle. Indeed, maximum
entropy solutions rely on “information” arguments which essentially state that
the maximum entropy distribution is the one satisfying the given constraints
and containing the largest amount of “uncertainty” (Jaynes, 1982). In plain
words, this principle guarantees that this solution does not satisfy
additional constraints that are undesired and unnecessary. We shall make the
rather mild assumption that the process $\theta$ is zero-mean and absolutely
continuous with respect to the Lebesgue measure. We will see that the optimal
solution (i.e. maximizing the differential entropy) is a Gaussian process
where $\theta_{g}$ and $\theta_{f}$ are independent. Then, we will also
characterize the corresponding kernel functions.

We start with the constraints on $\theta_{g}$ inducing BIBO stability on
$\hat{G}(z)$. Let $P\in\mathcal{S}_{2}(\mathbb{N})$ be a strictly positive
definite kernel (in the sense of Moore) such that $[P]_{t,t}\leq\kappa
t^{\alpha}e^{-\beta t}$, $t\in\mathbb{N}$, with $\kappa,\beta>0$ and
$\alpha\in\mathbb{R}$. Let also $\vartheta$ be a zero-mean process which
satisfies the moment constraint

$\displaystyle\mathbb{E}[\|\vartheta\|_{P^{-1}}^{2}]\leq c$

where $c\geq 0$. Then, for any $\varepsilon>1$ there exists
$\bar{\kappa}_{\varepsilon}>0$ such that the covariance function (kernel) $K$
of $\vartheta$ satisfies (Zorzi & Chiuso, 2017, Proposition 4):

$\displaystyle[K]_{t,t}\leq\bar{\kappa}_{\varepsilon}t^{\alpha+\varepsilon}e^{-\beta
t},\;\;t\in\mathbb{N}.$

It is not difficult to see that the condition above on $P$ is satisfied by the
kernels usually employed in the identification of dynamical models (e.g.
stable spline, tuned/correlated and so on, see Pillonetto et al. (2014)).
Thus, we consider the constraints

$\displaystyle\mathbb{E}[\|g^{[hk,jl]}\|^{2}_{P^{-1}}]\leq\tilde{c}_{hk,jl},\;\forall\,h,j\in
I_{1},\;\forall\,k,l\in I_{2}$ (72)

where $c_{hk,jl}\geq 0$ and $P\in\mathcal{S}_{2}(\mathbb{N})$ as above. Then,
by (72), the covariance of the $k$-th element of $g^{[hk,jl]}$ decays
exponentially. Accordingly, its posterior mean is a BIBO stable transfer
function.

######  Remark 4.

The kernel $P$ induces a prior about the decay rate of the impulse responses
$g^{[hk,jl]}$. On the other hand, the hyperparameters of $P$ are estimated
from the measured data, see Section 5. Therefore, although $P$ could assign
high probability to a small subset of BIBO stable systems, the latter is a
“wise” set because it has been extrapolated from the data.

It remains to derive the constraints inducing a sparse Kronecker decomposition
as in (13), i.e.:

$\displaystyle[E_{1}]_{h,j}=0\iff g^{[hk,jl]}=0,\;\forall\,k,l\in I_{2}$ (73)
$\displaystyle[E_{2}]_{k,l}=0\iff g^{[hk,jl]}=0,\;\forall\,h,j\in I_{1}.$ (74)

We consider the constraints:

$\displaystyle\sum_{k,l\in I_{2}}\mathbb{E}[\|g^{[hk,jl]}\|^{2}_{P^{-1}}]\leq
c_{hj}^{\prime},\;\forall\,h,j\in I_{1};$ (75) $\displaystyle\sum_{h,j\in
I_{1}}\mathbb{E}[\|g^{[hk,jl]}\|^{2}_{P^{-1}}]\leq
c_{kl}^{\prime\prime},\;\forall\,k,l\in I_{2}$ (76)

where $c_{hj}^{\prime},c_{kl}^{\prime\prime}\geq 0$. First, notice that (75)
and (76) imply (72) provided that
$\max\\{c_{hj}^{\prime},c_{kl}^{\prime\prime}\\}\leq\tilde{c}_{hk,jl}$ for any
$h,k,j,l$. Accordingly, the posterior mean of $G(z)$, under the constraints
(75)-(76), is BIBO stable. Let $c_{hj}^{\prime}=0$, by (75) we have that
$g^{[hk,jl]}$ $\forall\,k,l\in I_{2}$ are the null sequence in mean square and
so are their posterior mean, i.e. the latter satisfy (73). Let
$c_{kl}^{\prime\prime}=0$, by (76) we have that $g^{[hk,jl]}$ $\forall\,h,j\in
I_{1}$ are the null sequence in mean square and so are their posterior mean,
i.e. the latter satisfy (74). In a similar way, condition (13) on $F(z)$ can
be written as:

$\displaystyle[A_{1}]_{h}=0\iff f^{[hk,i]}=0,\;\forall\,k\in
I_{2},\;\forall\,i\in I_{u}$ $\displaystyle[A_{2}]_{k,i}=0\iff
f^{[hk,i]}=0,\;\forall\,h\in I_{1}.$

Thus, the constraints inducing BIBO stability and a sparse Kronecker
decomposition are:

$\displaystyle\sum_{\begin{subarray}{c}k\in I_{2}\\\ i\in
I_{u}\end{subarray}}\mathbb{E}[\|f^{[hk,i]}\|^{2}_{R^{-1}}]\leq
d_{h}^{\prime},\;\forall\,h\in I_{1};$ (77) $\displaystyle\sum_{h\in
I_{1}}\mathbb{E}[\|f^{[hk,i]}\|^{2}_{R^{-1}}]\leq
d_{ki}^{\prime\prime},\;\forall\,k\in I_{2},\;\forall\,i\in I_{u}$ (78)

where $d_{h}^{\prime},d_{ki}^{\prime\prime}\geq 0$ and
$R\in\mathcal{S}_{2}(\mathbb{N})$ is a strictly positive kernel such that
$[R]_{t,t}\leq\kappa t^{\alpha}e^{-\beta t}$, $t\in\mathbb{N}$,
$\kappa,\beta>0$ and $\alpha\in\mathbb{R}$.

In order to build the desired prior distribution we make use of the Kolmogorov
extension Theorem, see Øksendal (1998), and work with finite vectors extracted
from process $\theta=[\,\theta_{g}^{\top}\;\theta_{f}^{\top}\,]^{\top}$. More
precisely, consider a finite index set
$\mathcal{I}=\mathcal{I}_{g}\times\mathcal{I}_{f}$ and
$\mathcal{I}_{g},\mathcal{I}_{f}\subset\mathbb{N}$. Let $\check{\theta}$ be
the random vector whose components are extracted, from processes $\theta_{g}$
and $\theta_{f}$ according to the index sets $\mathcal{I}_{g}$ and
$\mathcal{I}_{f}$, respectively. We denote by
$\mathbf{p}_{\mathcal{I}}(\check{\theta})$ the probability density of
$\check{\theta}$. In view of the Kolmogorov extension Theorem, process
$\theta$ can be characterized by specifying the joint probability density
$\mathbf{p}_{\mathcal{I}}$ for all finite sets
$\mathcal{I}\subset\mathbb{N}\times\mathbb{N}$. Thus, the maximum entropy
process $\theta$ can be constructed building all the marginals
$\mathbf{p}_{\mathcal{I}}$ using the maximum entropy principle, which can thus
be extended by the Kolmogorov extension theorem. Such principle states that
among all the probability densities $\mathbf{p}_{\mathcal{I}}$ satisfying the
desired constraints, the optimal one should maximize the differential entropy
(Cover & Thomas, 1991). Let $P_{\mathcal{I}_{g}}$ and $R_{\mathcal{I}_{f}}$ be
the kernel matrices whose entries are extracted from $P$ and $R$ according to
$\mathcal{I}_{g}$ and $\mathcal{I}_{f}$, respectively. Therefore, the
corresponding maximum entropy problem is

$\displaystyle\underset{\mathbf{p}_{\mathcal{I}}\in\mathcal{P}}{\max}$
$\displaystyle\mathbf{H}(\mathbf{p}_{\mathcal{I}})$ s.t. (75) (76) (77) (78)
hold for $\mathcal{I}$ (79)

where $\mathcal{P}$ is the class of probability densities in
$\mathbb{R}^{|\mathcal{I}|}$ which are bounded and taking positive values and
$\mathbf{H}(\mathbf{p}_{\mathcal{I}})$ denotes the differential entropy of
$\check{\theta}$. Moreover, the constraints in (4) must be understood as
follows: (75), (76), (77) and (78) hold with $\theta$, $P$ and $R$ replaced by
$\check{\theta}$, $P_{\mathcal{I}_{g}}$ and $R_{\mathcal{I}_{f}}$,
respectively.

######  Theorem 5.

Under the assumption that
$c_{hj}^{\prime},c_{kl}^{\prime\prime},d_{h}^{\prime},d_{ki}^{\prime\prime}>0$,
$h,j\in I_{1}$, $k,l\in I_{2}$ and $i\in I_{u}$, the unique optimal solution
to the maximum entropy problem (4) is such that
$\check{\theta}=[\,\check{\theta}_{g}^{\top}\;\check{\theta}_{f}^{\top}\,]^{\top}$
is Gaussian with zero mean. Moreover, $\check{\theta}_{g}$ and
$\check{\theta}_{f}$ are independent and with kernel matrix, respectively,

$\displaystyle\begin{aligned} \check{K}_{g}&=X_{g}\otimes
P_{\mathcal{I}_{g}}\\\ \check{K}_{f}&=X_{f}\otimes R_{\mathcal{I}_{f}}\\\
X_{g}&=(\Lambda\otimes\Gamma)(\Lambda\otimes
I_{p_{2}^{2}}+I_{p_{1}^{2}}\otimes\Gamma)^{-1}\\\
X_{f}&=(\Pi\otimes\Omega)(\Pi\otimes
I_{p_{2}m}+I_{p_{1}}\otimes\Omega)^{-1}\\\
\Lambda&=\mathrm{diag}(\lambda_{11}\ldots\lambda_{1p_{1}}\ldots\lambda_{p_{1}1}\ldots\lambda_{p_{1}p_{1}})\\\
\Gamma&=\mathrm{diag}(\gamma_{11}\ldots\gamma_{1p_{2}}\ldots\gamma_{p_{2}1}\ldots\gamma_{p_{2}p_{2}})\\\
\Pi&=\mathrm{diag}(\pi_{1}\ldots\pi_{p_{1}})\\\
\Omega&=\mathrm{diag}(\omega_{11}\ldots\omega_{1m}\ldots\omega_{p_{2}1}\ldots\omega_{p_{2}m})\end{aligned}$
(80)

where $\lambda_{hj},\gamma_{kl},\pi_{h},\omega_{ki}>0$.

######  Proof 4.1.

We prove the claim by using the duality theory. The Lagrangian is

$\displaystyle L(\mathbf{p}_{\mathcal{I}},$
$\displaystyle\tilde{\Lambda},\tilde{\Gamma},\tilde{\Pi},\tilde{\Omega})=\mathbf{H}(\mathbf{p}_{\mathcal{I}})$
$\displaystyle+\frac{1}{2}\sum_{h,j\in
I_{1}}\tilde{\lambda}_{hj}\left(c_{hj}^{\prime}-\sum_{k,l\in
I_{2}}\mathbb{E}[\|\check{g}^{[hk,jl]}\|^{2}_{P_{\mathcal{I}_{g}}^{-1}}]\right)$
$\displaystyle+\frac{1}{2}\sum_{k,l\in
I_{2}}\tilde{\gamma}_{kl}\left(c_{kl}^{\prime\prime}-\sum_{h,j\in
I_{1}}\mathbb{E}[\|\check{g}^{[hk,jl]}\|^{2}_{P_{\mathcal{I}_{g}}^{-1}}]\right)$
$\displaystyle+\frac{1}{2}\sum_{h\in
I_{1}}\tilde{\pi}_{h}\left(d_{h}^{\prime}-\sum_{\begin{subarray}{c}k\in
I_{2}\\\ i\in
I_{u}\end{subarray}}\mathbb{E}[\|\check{f}^{[hk,i]}\|^{2}_{R_{\mathcal{I}_{f}}^{-1}}]\right)$
$\displaystyle+\frac{1}{2}\sum_{\begin{subarray}{c}k\in I_{2}\\\ i\in
I_{u}\end{subarray}}\tilde{\omega}_{ki}\left(d_{ki}^{\prime\prime}-\sum_{h\in
I_{1}}\mathbb{E}[\|\check{f}^{[hk,i]}\|^{2}_{R_{\mathcal{I}_{f}}^{-1}}]\right)$
$\displaystyle=-\mathbb{E}[\log(\mathbf{p}_{\mathcal{I}})]+\frac{1}{2}\sum_{h,j\in
I_{1}}\tilde{\lambda}_{hj}c_{hj}^{\prime}+\frac{1}{2}\sum_{k,l\in
I_{2}}\tilde{\gamma}_{kl}c_{kl}^{\prime\prime}$
$\displaystyle-\frac{1}{2}\sum_{\begin{subarray}{c}h,j\in I_{1}\\\ k,l\in
I_{2}\end{subarray}}(\tilde{\lambda}_{hj}+\tilde{\gamma}_{kl})\mathbb{E}[\|\check{g}^{[hk,jl]}\|^{2}_{P_{\mathcal{I}_{g}}^{-1}}]$
$\displaystyle+\frac{1}{2}\sum_{h\in
I_{1}}\tilde{\pi}_{h}d_{h}^{\prime}+\frac{1}{2}\sum_{\begin{subarray}{c}k\in
I_{2}\\\ i\in I_{u}\end{subarray}}\tilde{\omega}_{ki}d_{ki}^{\prime\prime}$
$\displaystyle-\frac{1}{2}\sum_{\begin{subarray}{c}h\in I_{1},k\in I_{2}\\\
i\in
I_{u}\end{subarray}}(\tilde{\pi}_{h}+\tilde{\omega}_{ki})\mathbb{E}[\|\check{f}^{[hk,i]}\|^{2}_{R_{\mathcal{I}_{f}}^{-1}}]$
$\displaystyle=-\mathbb{E}[\log(\mathbf{p}_{\mathcal{I}})]+\frac{1}{2}\sum_{h,j\in
I_{1}}\tilde{\lambda}_{hj}c_{hj}^{\prime}+\frac{1}{2}\sum_{k,l\in
I_{2}}\tilde{\gamma}_{kl}c_{kl}^{\prime\prime}$
$\displaystyle+\frac{1}{2}\sum_{h\in
I_{1}}\tilde{\pi}_{h}d_{h}^{\prime}-\frac{1}{2}\mathbb{E}[\|\check{\theta}_{g}\|^{2}_{(\tilde{\Lambda}\otimes
I+I\otimes\tilde{\Gamma})\otimes P_{\mathcal{I}_{g}}^{-1}}]$
$\displaystyle+\frac{1}{2}\sum_{\begin{subarray}{c}k\in I_{2}\\\ i\in
I_{u}\end{subarray}}\tilde{\omega}_{ki}d_{ki}^{\prime\prime}-\frac{1}{2}\mathbb{E}[\|\check{\theta}_{f}\|^{2}_{(\tilde{\Pi}\otimes
I+I\otimes\tilde{\Omega})\otimes R_{\mathcal{I}_{f}}^{-1}}]$ (81)

where
$\tilde{\Lambda}=\mathrm{diag}(\tilde{\lambda}_{11}\ldots\tilde{\lambda}_{1p_{1}}\ldots\tilde{\lambda}_{p_{1}1}\ldots\tilde{\lambda}_{p_{1}p_{1}})$,
$\tilde{\Gamma}=\mathrm{diag}(\tilde{\gamma}_{11}\ldots\tilde{\gamma}_{1p_{2}}\ldots\tilde{\gamma}_{p_{2}1}\ldots\tilde{\gamma}_{p_{2}p_{2}})$,
$\tilde{\Pi}=\mathrm{diag}(\tilde{\pi}_{1}\ldots\ldots\tilde{\pi}_{p_{1}})$,
$\tilde{\Omega}=\mathrm{diag}(\tilde{\gamma}_{11}\ldots\tilde{\gamma}_{1m}\ldots\tilde{\gamma}_{p_{2}1}\ldots\tilde{\gamma}_{p_{2}m})$;
moreover,
$\tilde{\lambda}_{h,j},\tilde{\gamma}_{kl},\tilde{\pi}_{h},\tilde{\omega}_{hi}\geq
0$ are the Lagrange multipliers corresponding to (75), (76), (77) and (78),
respectively. It is not difficult to see that (4.1) is strictly concave in
$\mathcal{P}$. Moreover, its point of maximum exists under the assumptions
that $\tilde{\lambda}_{hj}+\tilde{\gamma}_{kl}>0$ and
$\tilde{\pi}_{h}+\tilde{\omega}_{ki}>0$, and it is Gaussian distributed:

$\displaystyle\mathbf{p}_{\mathcal{I}}(\check{\theta})=\mathbf{p}_{\mathcal{I}_{g}}(\check{\theta}_{g})\mathbf{p}_{\mathcal{I}_{f}}(\check{\theta}_{f})$
(82)

where

$\displaystyle\mathbf{p}_{\mathcal{I}_{f}}(\check{\theta}_{g})$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{p_{1}^{2}p_{2}^{2}|\mathcal{I}_{g}|}|\check{K}_{g}|}}\exp(-\frac{1}{2}\check{\theta}_{g}^{\top}\check{K}_{g}^{-1}\check{\theta}_{g})$
$\displaystyle\mathbf{p}_{\mathcal{I}_{f}}(\check{\theta}_{f})$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{p_{1}p_{2}m|\mathcal{I}_{f}|}|\check{K}_{f}|}}\exp(-\frac{1}{2}\check{\theta}_{f}^{\top}\check{K}_{f}^{-1}\check{\theta}_{f})$

and

$\displaystyle\check{K}_{g}$ $\displaystyle=(\tilde{\Lambda}\otimes
I+I\otimes\tilde{\Gamma})^{-1}\otimes P_{\mathcal{I}_{g}}$
$\displaystyle\check{K}_{f}$ $\displaystyle=(\tilde{\Pi}\otimes
I+I\otimes\tilde{\Omega})^{-1}\otimes R_{\mathcal{I}_{f}}.$

In the case that $\tilde{\lambda}_{hj}+\tilde{\gamma}_{kl}=0$ or
$\tilde{\pi}_{h}+\tilde{\omega}_{ki}=0$ for some $h,k,j,l,i$, then the maximum
cannot be attained in $\mathcal{P}$. Substituting the optimal form of
$\mathbf{p}_{\mathcal{I}}$ in the Lagrangian we obtain the dual function (up
to terms not depending on $\tilde{\Lambda}$, $\tilde{\Gamma}$, $\tilde{\Pi}$
and $\tilde{\Omega}$):

$\displaystyle
J(\tilde{\Lambda},\tilde{\Gamma},\tilde{\Pi},\tilde{\Omega})=J_{f}(\tilde{\Lambda},\tilde{\Gamma})+J_{g}(\tilde{\Pi},\tilde{\Omega})$

where

$\displaystyle J_{f}(\tilde{\Lambda},\tilde{\Gamma})=-$
$\displaystyle\frac{|\mathcal{I}_{g}|}{2}\sum_{\begin{subarray}{c}h,j\in
I_{1}\\\ k,l\in
I_{2}\end{subarray}}\log(\tilde{\lambda}_{hj}+\tilde{\gamma}_{kl})$
$\displaystyle\hskip 14.22636pt+\frac{1}{2}\sum_{h,j\in
I_{1}}\tilde{\lambda}_{hj}c_{hj}^{\prime}+\frac{1}{2}\sum_{k,l\in
I_{2}}\tilde{\gamma}_{kl}c_{kl}^{\prime\prime}.$ $\displaystyle
J_{g}(\tilde{\Pi},\tilde{\Omega})=-$
$\displaystyle\frac{|\mathcal{I}_{f}|}{2}\sum_{\begin{subarray}{c}h\in
I_{1},\,k\in I_{2}\\\ i\in
I_{u}\end{subarray}}\log(\tilde{\pi}_{h}+\tilde{\omega}_{ki})$
$\displaystyle\hskip 14.22636pt+\frac{1}{2}\sum_{h\in
I_{1}}\tilde{\pi}_{hj}d_{h}^{\prime}+\frac{1}{2}\sum_{\begin{subarray}{c}k\in
I_{1}\\\ i\in I_{u}\end{subarray}}\tilde{\omega}_{ki}d_{ki}^{\prime\prime}.$

Hence, we can minimize $(\tilde{\Lambda},\tilde{\Gamma})$ and
$(\tilde{\Pi},\tilde{\Omega})$ in an independent way. We start to minimize
$J_{f}(\tilde{\Lambda},\tilde{\Gamma})$ over the open and unbounded set

$\displaystyle\mathcal{C}=\\{\,(\tilde{\Lambda},\tilde{\Gamma})\in\mathcal{D}\hbox{
s.t. }\tilde{\Lambda}\geq 0,\,\tilde{\Gamma}\geq 0,\,\tilde{\Lambda}\otimes
I+I\otimes\tilde{\Gamma}>0\,\\}$

where $\mathcal{D}=\mathcal{D}_{p_{1}^{2}}\times\mathcal{D}_{p_{2}^{2}}$.

Next, we show that we can restrict the search of the minimum of $J_{f}$ over a
compact set and thus by the Weierstrass’ theorem, the dual problem admits
solution. Take a sequence
$(\tilde{\Lambda}^{(n)},\tilde{\Gamma}^{(n)})\in\mathcal{C}$,
$n\in\mathbb{N}$, such that $\|\tilde{\Lambda}^{(n)}\|\rightarrow\infty$
and/or $\|\tilde{\Gamma}^{(n)}\|\rightarrow\infty$ as $n\rightarrow\infty$.
Then, $J_{f}(\tilde{\Lambda}^{(n)},\tilde{\Gamma}^{(n)})\rightarrow\infty$
because the linear terms dominate the logarithmic ones. Therefore, such a
sequence cannot be an infimizing sequence. Accordingly, we can restrict the
set $\mathcal{C}$ as

$\displaystyle\mathcal{C}_{1}=\\{\,(\tilde{\Lambda}$
$\displaystyle,\tilde{\Gamma})\in\mathcal{D}\hbox{ s.t.
}0\leq\tilde{\Lambda}\leq\lambda_{MAX}I,$ $\displaystyle
0\leq\Gamma\leq\gamma_{MAX}I,\,\tilde{\Lambda}\otimes
I+I\otimes\tilde{\Gamma}>0\,\\}$

where $\lambda_{MAX},\gamma_{MAX}>0$ are constants taken sufficiently large.
Take a sequence
$(\tilde{\Lambda}^{(n)},\tilde{\Gamma}^{(n)})\in\mathcal{C}_{1}$,
$n\in\mathbb{N}$, such that
$\lambda^{(n)}_{hj}+\tilde{\gamma}_{kl}^{(n)}\rightarrow 0$ as
$n\rightarrow\infty$ for some $h,j\in I_{1}$ and $k,l\in I_{2}$. Then, we have

$\displaystyle
J_{f}(\tilde{\Lambda}^{(n)},\tilde{\Gamma}^{(n)})\geq-\frac{|\mathcal{I}_{g}|}{2}\sum_{\begin{subarray}{c}h,j\in
I_{1}\\\ k,l\in
I_{2}\end{subarray}}\log(\tilde{\lambda}_{hj}+\tilde{\gamma}_{kl})\rightarrow\infty.$

Accordingly, it cannot be an infimizing sequence, and thus we can restrict the
search to the set

$\displaystyle\mathcal{C}_{2}=\\{\,(\tilde{\Lambda}$
$\displaystyle,\tilde{\Gamma})\in\mathcal{D}\hbox{ s.t.
}0\leq\tilde{\Lambda}\leq\lambda_{MAX}I,$ $\displaystyle
0\leq\Gamma\leq\gamma_{MAX}I,\,\tilde{\Lambda}\otimes
I+I\otimes\tilde{\Gamma}\geq\varepsilon I\,\\}$

where $\epsilon$ is a constant taken sufficiently small. Thus,
$\mathcal{C}_{2}$ is closed and bounded and thus compact. We conclude that
$J_{f}$ admits minimum over $\mathcal{C}$. In a similar way, it is possible to
prove that $J_{g}$ admits minimum over the set

$\displaystyle\mathcal{C}=\\{\,(\tilde{\Pi},\tilde{\Omega})\in\tilde{\mathcal{D}}\hbox{
s.t. }\tilde{\Pi}\geq 0,\,\tilde{\Omega}\geq 0,\,\tilde{\Pi}\otimes
I+I\otimes\tilde{\Omega}>0\,\\}$

with $\tilde{\mathcal{D}}=\mathcal{D}_{p_{1}}\times\mathcal{D}_{p_{2}m}$. We
conclude that the dual problem admits solution. Accordingly, (82) is the
solution to the maximum entropy problem (4). Finally, by defining
$\Lambda=\tilde{\Lambda}^{-1}$, $\Gamma=\tilde{\Gamma}^{-1}$,
$\Pi=\tilde{\Pi}^{-1}$ and $\Omega=\tilde{\Omega}^{-1}$ we obtain (80). ∎

We assumed that constraints (75), (76), (77) and (78) are totally binding in
problem (4). However, we are interested in the limiting cases where
$c_{hj}^{\prime}$, $c_{kl}^{\prime\prime}$, $d_{h}^{\prime}$ and
$d_{ki}^{\prime\prime}$ might be equal to zero in order to obtain a posterior
corresponding to a Kronecker model. To include these scenarios, we consider
the limits as $c_{hj}^{\prime}\rightarrow 0$,
$c_{kl}^{\prime\prime}\rightarrow 0$, $d_{h}^{\prime}\rightarrow 0$,
$d_{ki}^{\prime\prime}\rightarrow 0$ and extend the maximum entropy solution
by continuity.

######  Proposition 6.

Let $\mathcal{E}_{1}:=\\{(h,j)\hbox{ s.t. }c_{hj}^{\prime}>0\\}$,
$\mathcal{E}_{2}:=\\{(k,l)\hbox{ s.t. }c_{kl}^{\prime\prime}>0\\}$. If the
sets $\mathcal{E}_{1},\mathcal{E}_{2}$ are nonempty, then the maximum entropy
solution extended by continuity is the probability density such that
$\check{\theta}_{g}$ is Gaussian, zero-mean and with kernel matrix as in (80)
where:

  * •

$(i,j)\notin\mathcal{E}_{1}$ $\implies$ $\lambda_{hj}=0$,

  * •

$(k,l)\notin\mathcal{E}_{2}$ $\implies$ $\gamma_{kl}=0$.

######  Proof 4.2.

First, recall that $\lambda_{hj}=\tilde{\lambda}_{hj}^{-1}$,
$\gamma_{kl}=\tilde{\gamma}_{kl}^{-1}$ and $\tilde{\lambda}_{hj}$,
$\tilde{\gamma}_{kl}$ are the original Lagrange multipliers. The maximum
entropy solution $\check{\theta}_{ME}$ is such that, by Theorem 5,
$\check{g}_{ME}^{[hk,jl]}$ is a zero-mean process with kernel matrix

$\displaystyle\frac{1}{\tilde{\lambda}_{hj}+\tilde{\gamma}_{kl}}P_{\mathcal{I}_{g}}.$

Moreover, if $\tilde{\gamma}_{kl}>0$ ($\lambda_{hj}>0$) then the maximum
entropy solution $\check{\theta}_{ME}$ satisfies the corresponding constraint
with equality. We consider two sequences
$\\{c_{hj}^{\prime(n)}\\}_{n\in\mathbb{N}}$, $c_{hj}^{\prime(n)}>0$, and
$\\{c_{kl}^{\prime\prime(n)}\\}_{n\in\mathbb{N}}$,
$c_{kl}^{\prime\prime(n)}>0$, such that $c_{hj}^{\prime(n)}\rightarrow
c_{hj}^{\prime}$ and $c_{kl}^{\prime\prime(n)}\rightarrow
c_{kl}^{\prime\prime}$ as $n\rightarrow\infty$. For any $c_{hj}^{\prime(n)}$,
we have

$\displaystyle c_{hj}^{\prime(n)}$ $\displaystyle\geq\sum_{k,l\in
I_{2}}\mathbb{E}[\|\check{g}_{ME}^{[hk,jl]}\|^{2}_{P_{\mathcal{I}_{g}}^{-1}}]$
$\displaystyle=\sum_{k,l\in I_{2}}\mathop{\rm
tr}\left(\mathbb{E}[\check{g}_{ME}^{[hk,jl]}\check{g}_{ME}^{[hk,jl]\top}]P^{-1}_{\mathcal{I}_{g}}\right)$
$\displaystyle=|\mathcal{I}_{g}|\sum_{k,l\in
I_{2}}\frac{1}{\tilde{\lambda}_{hj}^{(n)}+\tilde{\gamma}_{kl}^{(n)}}>0.$ (83)

In a similar way, we have

$\displaystyle c_{kl}^{\prime\prime(n)}$
$\displaystyle\geq|\mathcal{I}_{g}|\sum_{h,j\in
I_{1}}\frac{1}{\tilde{\lambda}_{hj}^{(n)}+\tilde{\gamma}_{kl}^{(n)}}>0.$ (84)

Let $(h,j)\notin\mathcal{E}_{1}$, then $c^{\prime(n)}_{hj}\rightarrow 0$.
Then, we consider the arbitrary subsequences $\tilde{\lambda}^{(n_{r})}_{hj}$
and $\tilde{\gamma}^{(n_{r})}_{kl}$, with $r\in\mathbb{N}$, which admit limit.
Taking into account (4.2), we may have two possible cases. First case.
$\tilde{\gamma}_{kl}^{(n_{r})}$ converges for some $k,l\in I_{2}$. By (4.2),
this implies that $\tilde{\lambda}_{hj}^{(n_{r})}\rightarrow\infty$ and thus
$\lambda_{hj}^{(n_{r})}\rightarrow 0$. Second case.
$\tilde{\gamma}_{kl}^{(n_{r})}\rightarrow\infty$ for any $k,l\in I_{2}$. This
means that there exists $\bar{r}\in\mathbb{N}$ sufficiently large for which
$\tilde{\gamma}_{kl}^{(n_{r})}>0$ for any $k,l\in I_{2}$ and $r\geq\bar{r}$.
Accordingly, the corresponding constraints in (84) are satisfied with equality
for $r\geq\bar{r}$ and $n=n_{r}$. Hence
$c_{kl}^{\prime\prime(n_{r})}\rightarrow 0$ for any $k,l\in I_{2}$ and thus
$\mathcal{E}_{2}=\emptyset$. This is not possible because
$\mathcal{E}_{2}\neq\emptyset$ by assumption. We conclude that
$c_{hj}^{\prime(n_{r})}\rightarrow 0$ implies
$\lambda_{hj}^{(n_{r})}\rightarrow 0$ for any arbitrary subsequence which
admits limit and thus $\lambda_{hj}^{(n)}\rightarrow 0$. Accordingly, in a
similar way, it is possible to prove that $c_{kl}^{\prime\prime(n)}\rightarrow
0$ implies $\gamma_{kl}^{(n)}\rightarrow 0$. ∎

######  Proposition 7.

Let $\mathcal{A}_{1}:=\\{(h)\hbox{ s.t. }d_{h}^{\prime}>0\\}$,
$\mathcal{A}_{2}:=\\{(k,i)\hbox{ s.t. }d_{ki}^{\prime\prime}>0\\}$. If the
sets $\mathcal{A}_{1},\mathcal{A}_{2}$ are nonempty, then the maximum entropy
solution extended by continuity is the probability density such that
$\check{\theta}_{f}$ is Gaussian, zero-mean and with kernel matrix as in (80)
where

  * •

$h\notin\mathcal{A}_{1}$ $\implies$ $\pi_{h}=0$,

  * •

$(k,l)\notin\mathcal{A}_{2}$ $\implies$ $\omega_{ki}=0$.

######  Proof 4.3.

The proof is similar to one of Proposition 6. ∎

Finally, in view of the Kolmogorov extension Theorem, from the probability
density of $\check{\theta}$ we can characterize the probability law of
$\theta$ maximizing the differential entropy.

######  Corollary 8.

Consider the zero-mean Gaussian process
$\theta=[\,\theta_{g}^{\top}\;\,\theta_{f}^{\top}]^{\top}$ with kernel,
respectively,

$\displaystyle\begin{aligned} K_{g}&=X_{g}\otimes P\\\ K_{f}&=X_{f}\otimes
R\\\ \end{aligned}$ (85)

where $X_{g}$ and $X_{f}$ are defined as in (80). For all finite sets
$\mathcal{I}\subset\mathbb{N}\times\mathbb{N}$, its joint probability density
is the extended solution to the maximum entropy problem (4).

It is worth noting that the maximum entropy solution is such that
$g^{[hk,jl]}$ and $f^{[hk,i]}$ are zero-mean Gaussian processes with kernel
function, respectively,

$\displaystyle
K_{g}^{[hk,jl]}:=\frac{\lambda_{hj}\gamma_{kl}}{\lambda_{hj}+\gamma_{kl}}P,\;\;K_{f}^{[hk,i]}:=\frac{\pi_{h}\omega_{ki}}{\pi_{h}+\omega_{ki}}R.$
(86)

In the case that both $\lambda_{hj}$ and $\gamma_{kl}$ are equal to zero, we
define by continuity $\lambda_{hj}\gamma_{kl}/(\lambda_{hj}+\gamma_{kl})=0$.
Accordingly, $K_{g}^{[hk,jl]}=0$ if $\lambda_{hj}=0$ or $\gamma_{kl}=0$.
Similarly, in the case that $\pi_{h}$ or $\omega_{ki}$ is equal to zero, we
define $\pi_{h}\omega_{ki}/(\pi_{h}+\omega_{ki})=0$.

In Chiuso & Pillonetto (2012) the transfer matrices in model (12) are modeled
as zero mean Gaussian processes with kernels

$\displaystyle
K_{g}^{[hk,jl]}:=\lambda_{hjkl}P,\;\;K_{f}^{[hk,i]}:=\pi_{hki}R$ (87)

where $\lambda_{hjkl},\pi_{hki}\geq 0$, with $h,j\in I_{1}$, $k,l\in I_{2}$
and $i\in I_{u}$, are the hyperparameters representing the so called scale
factors. The prior in (87) has been used to learn sparse networks. Finally, it
is worth noting that there are many more hyperparameters in (87) compared with
(86).

######  Remark 9.

It is worth noting that the constants $c_{hj}^{\prime}$’s,
$c_{kl}^{\prime\prime}$’s, $d_{h}^{\prime}$’s and $c_{ki}^{\prime\prime}$’s
were needed only to derive the kernel structure in (85). Indeed, we are not
interested in their particular values because, as we will show in the next
section, the hyperparameters of the kernel will be found through another
paradigm which exploits the measured data.

Hierarchical dynamic networks. An interesting scenario included in this
framework is the case of a system without input (i.e. $y(t)$ is a stationary
stochastic process which is characterized by the transfer matrix $G(z)$) with
$p_{1}=p_{2}$ and $E_{1}=E_{2}$. It is then natural to take

$\displaystyle K_{g}=(\Lambda\otimes\Lambda)(\Lambda\otimes
I_{p_{1}^{2}}+I_{p_{1}^{2}}\otimes\Lambda)^{-1}\otimes P$

and thus $g^{[hk,jl]}$ is modeled is a Gaussian process with kernel function

$\displaystyle
K_{g}^{[hk,jl]}:=\frac{\lambda_{hj}\lambda_{kl}}{\lambda_{hj}+\lambda_{kl}}P.$
(88)

This model corresponds to a network which is hierarchically organized into
communities (clusters), the communities then grow, creating miniature copies
of themselves. This idea can be applied recursively at the price that the
dimension of $G(z)$ increases exponentially. Thus, $\Lambda$ is the Kronecker
initiator for the prior on $G(z)$. These models have been proposed in Leskovec
et al. (2010): the main difference with respect to our approach is the fact
that we do not impose the Kronecker decomposition on $G(z)$, but on
$\mathrm{supp}(G(z))$; moreover $y(t)$ is not white noise (i.e. we consider a
process with dynamics).

##  5 Hyperparameters Estimation

In order to compute $\hat{\theta}_{g}$ and $\hat{\theta}_{f}$, we have to
estimate $\Sigma$ and $\zeta$. The latter contains the hyperparameters of
$P,R$ and $\xi:=\\{\Lambda,\Gamma,\Pi,\Omega\\}$. The noise covariance matrix
$\Sigma$ can be estimated using a low-bias ARX-model as suggested in (Goodwin
et al., 1992). To estimate the hyperparameters of $P$ and $R$ we consider a
nonparametric unstructured model, i.e. we consider model (12) but we do not
require that $G(z)$ and $F(z)$ satisfy condition (13). Then, the
hyperparameters of $P$ can be estimated by minimizing the negative log-
marginal likelihood of $\mathrm{y}^{+}$ computed using the unstructured model,
see Pillonetto et al. (2011). Throughout the paper, kernel $P$ is chosen as
Stable Spline (SS) with a modulation factor (Zorzi & Chiuso, 2018):

$\displaystyle[P]_{t,s}=\left[\frac{e^{-\beta(t+s)}e^{-\beta\max(t,s)}}{2}-\frac{e^{-3\beta\max(t,s)}}{6}\right]\cos(\omega_{0}(t-s))$

where $\beta\in(0,1)$ and $\omega_{0}\in(0,\pi)$ tune, respectively, the a
priori information about the decay rate and the frequency content of the
impulse responses. Kernel $R$ is chosen likewise. As explained in Pillonetto
et al. (2011), kernels of this type represent a good prior in the case that
the impulse responses of the predictor have “fast dynamics”.

Then, an estimate of $\xi$ is given by minimizing the negative log-marginal
likelihood $\ell$ of $\mathrm{y}^{+}$ under model (12), (85) with $\Sigma$,
$P$ and $R$ fixed as above. Under the assumptions in (3), we have (up to
constant terms)

$\displaystyle\ell(\mathrm{y}^{+},\xi)$ $\displaystyle=\log\det
V+\frac{1}{2}(\mathrm{y}^{+})^{\top}V^{-1}\mathrm{y}^{+}$ $\displaystyle V$
$\displaystyle=\Phi K_{g}\Phi^{\top}+\Psi K_{f}\Psi^{\top}+\Sigma\otimes I$
(89)

where $K_{g}$ and $K_{f}$ have been defined in (85) and $\xi$ belongs to the
compact set

$\displaystyle\Xi=\\{(\Lambda,\Gamma,\Pi,\Omega)\in\mathcal{D}\hbox{ s.t.
}0\leq\Lambda,\Gamma,\Pi,\Omega\leq\kappa I\\}$ (90)

where
$\mathcal{D}=\mathcal{D}_{p_{1}^{2}}\times\mathcal{D}_{p_{2}^{2}}\times\mathcal{D}_{p_{1}}\times\mathcal{D}_{p_{2}m}$
and $\kappa>0$ is a constant taken sufficiently large. Notice that the choice
of $\Xi$ is not restrictive. Indeed, if we take: $\Gamma\geq 0$ such that
$\gamma_{kl}>0$ for some $k,l\in I_{2}$; $\Lambda\geq 0$ such that
$\lambda_{hj}\rightarrow\infty$ for some $h,j\in I_{1}$. Then,
$\lambda_{hj}\gamma_{kl}/(\lambda_{hj}+\gamma_{kl})\rightarrow\infty$ and thus
$\ell(\mathrm{y}^{+},\xi)\rightarrow\infty$. This means that the minimum
cannot be attained for $\lambda_{hj}>\kappa$. In the case that $\Gamma=0$,
then we have
$\ell(\mathrm{y}^{+},\\{\Lambda,\Gamma,\Pi,\Omega\\})=\ell(\mathrm{y}^{+},\\{\Lambda+\tilde{\Lambda},\Gamma,\Pi,\Omega\\})$
for any $\tilde{\Lambda}\in\mathcal{D}_{p_{1}^{2}}$ such that
$\Lambda+\tilde{\Lambda}\geq 0$; accordingly, we can restrict the values of
$\Lambda$ as in (90). A similar reasoning holds for $\Gamma$, $\Pi$ and
$\Omega$. Since $\ell$ is a continuous function over the compact set $\Xi$, it
follows that $\ell$ admits a point of minimum over $\Xi$.

The minimization of (5) with respect to $\xi$ in the numerical examples
presented in the next section is performed through the nonlinear optimization
Matlab solver fmincon.m which uses a gradient based method and the gradient,
i.e. $\partial\;\ell(\mathrm{y}^{+},\xi)/\partial\xi$, is supplied. The
computation of the gradient can be done in a similar way to that explained in
(Chen & Ljung, 2013, Section 5). In the case one is interested to make more
efficient the algorithm, it is possible to perform the minimization by means
of the scaled gradient projection method proposed in Bonettini et al. (2015).
Finally, it is worth noting that only a local minimum can be found because
$\ell(\mathrm{y}^{+},\xi)$ is not convex with respect to $\xi$.

Next, we prove that the posterior mean of $\theta_{g}$ and $\theta_{f}$ under
the prior defined through kernels in (85) leads to a sparse Kronecker network.
The analysis is performed by drawing inspiration from Aravkin et al. (2014)
which deals with the case of a sparse Bayesian network (Chiuso & Pillonetto,
2012), i.e. the prior is defined through kernels in (87), and indeed our prior
is closely related to Auto Relevance Determination (ARD) method (MacKay,
1994). First, notice that the estimates $\hat{\theta}_{g}$ and
$\hat{\theta}_{f}$ describe a sparse Kronecker network only if
$\hat{\xi}=\\{\hat{\Lambda},\hat{\Gamma},\hat{\Pi},\hat{\Omega}\\}$ is such
that $\hat{\Lambda},\Gamma,\hat{\Pi},\hat{\Omega}$ are spare matrices. We
assume that the noise variance at each node is the same, i.e.
$\sigma_{hk}^{2}=\sigma^{2}$ with $h\in I_{1}$ and $k\in I_{2}$, and that

$\displaystyle(I\otimes P^{1/2})\Phi^{\top}\Phi(I\otimes P^{1/2})=NI,$
$\displaystyle(I\otimes P^{1/2})\Phi^{\top}\Psi(I\otimes R^{1/2})=0,$
$\displaystyle(I\otimes R^{1/2})\Psi^{\top}\Psi(I\otimes R^{1/2})=NI.$ (91)

Moreover $\sigma^{2}$, $P$ are $R$ are fixed. First, we consider the case that
$g^{[hk,jl]}$’s and $f^{[hk,i]}$’s are finite dimensional vectors, say $T$
their dimension. Accordingly, the corresponding transfer matrices $G(z)$ and
$F(z)$ are finite impulse responses (FIR) and $T$ represents their practical
length. Moreover, we assume that $N\geq T$. We define

$\displaystyle\hat{\theta}_{g,LS}$
$\displaystyle=(\Phi^{\top}\Phi)^{-1}\Phi^{\top}\mathrm{y}^{+}=\frac{1}{N}(I\otimes
P)\Phi^{\top}\mathrm{y}^{+},$ $\displaystyle\;\hat{\theta}_{f,LS}$
$\displaystyle=(\Psi^{\top}\Psi)^{-1}\Psi^{\top}\mathrm{y}^{+}=\frac{1}{N}(I\otimes
R)\Psi^{\top}\mathrm{y}^{+}$

which are the least squares estimators of $\theta_{g}$ and $\theta_{f}$. Thus,
we have

$\displaystyle\hat{g}^{[hk,jl]}_{LS}=\frac{1}{N}P\phi^{\top}_{jl}\mathrm{y}^{+}_{hk},\;\;\hat{f}^{[hk,i]}_{LS}=\frac{1}{N}R\psi^{\top}_{i}\mathrm{y}^{+}_{hk}.$
(92)

######  Proposition 10.

Let $\lambda_{hj}$, $\gamma_{k,l}$, $\pi_{h}$ and $\omega_{ki}$ be the entries
in the main diagonal of $\Lambda$, $\Gamma$, $\Pi$, $\Omega$, respectively.
The following facts hold:

  * •

if

$\displaystyle\frac{T}{N}\sigma^{2}-\|\hat{g}_{LS}^{[hk,jl]}\|^{2}_{P^{-1}}\geq
0,\;\;\forall\,k,l\in I_{2},$ (93)

then the subset

$\Xi_{1}=\\{(\Lambda,\Gamma,\Pi,\Omega)\hbox{ s.t.
}\Lambda,\Gamma,\Pi,\Omega\geq 0,\;\lambda_{hj}=0\\}$

contains a local minimum for $\ell$;

  * •

if

$\displaystyle\frac{T}{N}\sigma^{2}-\|\hat{g}_{LS}^{[hk,jl]}\|^{2}_{P^{-1}}\geq
0,\;\;\forall\,h,j\in I_{1},$

then the subset

$\Xi_{2}=\\{(\Lambda,\Gamma,\Pi,\Omega)\hbox{ s.t.
}\Lambda,\Gamma,\Pi,\Omega\geq 0,\;\gamma_{kl}=0\\}$

contains a local minimum for $\ell$;

  * •

if

$\displaystyle\frac{T}{N}\sigma^{2}-\|\hat{f}_{LS}^{[hk,i]}\|^{2}_{R^{-1}}\geq
0,\;\;\forall\,k\in I_{2},\;i\in I_{u},$

then the subset

$\Xi_{3}=\\{(\Lambda,\Gamma,\Pi,\Omega)\hbox{ s.t.
}\Lambda,\Gamma,\Pi,\Omega\geq 0,\;\pi_{h}=0\\}$

contains a local minimum for $\ell$;

  * •

if

$\displaystyle\frac{T}{N}\sigma^{2}-\|\hat{f}_{LS}^{[hk,i]}\|^{2}_{R^{-1}}\geq
0,\;\;\forall\,h\in I_{1},$ (94)

then the subset

$\Xi_{4}=\\{(\Lambda,\Gamma,\Pi,\Omega)\hbox{ s.t.
}\Lambda,\Gamma,\Pi,\Omega\geq 0,\;\omega_{ki}=0\\}$

contains a local minimum for $\ell$.

######  Proof 5.1.

Let

$\displaystyle X:=\left[\begin{array}[]{cc}X_{g}\otimes I&0\\\ 0&X_{f}\otimes
I\end{array}\right],\;\;\tilde{K}:=\left[\begin{array}[]{cc}I\otimes P&0\\\
0&I\otimes R\end{array}\right],$ (99)

thus

$\displaystyle K:=\left[\begin{array}[]{cc}K_{g}&0\\\
0&K_{f}\end{array}\right]=X\tilde{K}=\tilde{K}X=\tilde{K}^{1/2}X\tilde{K}^{1/2}$
(102)

and condition (5) can be written as

$\displaystyle\tilde{K}^{1/2}\left[\begin{array}[]{c}\Phi^{\top}\\\
\Psi^{\top}\end{array}\right]\left[\begin{array}[]{cc}\Phi&\Psi\end{array}\right]\tilde{K}^{1/2}=NI.$
(106)

In what follows, the symbol $\equiv$ means “equal to up to terms not depending
on $\xi$”. Next, we rewrite $\ell(\mathrm{y}^{+},\xi)$ using (92) and (106).
First,

$\displaystyle\log\det(\Phi K_{g}\Phi^{\top}+\Psi
K_{f}\Psi^{\top}+\sigma^{2}I)$
$\displaystyle\equiv\log\det\left(\frac{1}{\sigma^{2}}\left[\begin{array}[]{cc}\Phi&\Psi\end{array}\right]K\left[\begin{array}[]{c}\Phi^{\top}\\\
\Psi^{\top}\end{array}\right]+I\right)$ (110)
$\displaystyle\equiv\log\det\left(\frac{1}{\sigma^{2}}\left[\begin{array}[]{cc}\Phi&\Psi\end{array}\right]\tilde{K}^{1/2}X^{1/2}X^{1/2}\tilde{K}^{1/2}\left[\begin{array}[]{c}\Phi^{\top}\\\
\Psi^{\top}\end{array}\right]+I\right)$ (114)
$\displaystyle\equiv\log\det\left(\frac{1}{\sigma^{2}}X^{1/2}\tilde{K}^{1/2}\left[\begin{array}[]{c}\Phi^{\top}\\\
\Psi^{\top}\end{array}\right]\left[\begin{array}[]{cc}\Phi&\Psi\end{array}\right]\tilde{K}^{1/2}X^{1/2}+I\right)$
(118)
$\displaystyle=\log\det\left(\frac{N}{\sigma^{2}}\left[\begin{array}[]{cc}X_{g}\otimes
I_{T}&0\\\ 0&X_{f}\otimes I_{T}\end{array}\right]+I\right)$ (121)
$\displaystyle=\log\det\left(\frac{N}{\sigma^{2}}X_{g}\otimes
I_{T}+I\right)+\log\det\left(\frac{N}{\sigma^{2}}X_{f}\otimes I_{T}+I\right)$
$\displaystyle=T\sum_{\begin{subarray}{c}h\in I_{1}\\\ k\in
I_{2}\end{subarray}}\left[\sum_{\begin{subarray}{c}j\in I_{1}\\\ l\in
I_{2}\end{subarray}}\log\det\left(\frac{N}{\sigma^{2}}\frac{\lambda_{hj}\gamma_{kl}}{\lambda_{hj}+\gamma_{kl}}+1\right)\right.$
$\displaystyle\left.\hskip 11.38092pt+\sum_{\begin{subarray}{c}i\in
I_{u}\end{subarray}}\log\det\left(\frac{N}{\sigma^{2}}\frac{\pi_{h}\omega_{ki}}{\pi_{h}+\omega_{ki}}+1\right)\right]$
$\displaystyle(\mathrm{y}^{+})^{\top}(\Phi K_{g}\Phi^{\top}+\Psi
K_{f}\Psi^{\top}+\sigma^{2}I)^{-1}\mathrm{y}^{+}$
$\displaystyle=(\mathrm{y}^{+})^{\top}\left(\left[\begin{array}[]{cc}\Phi&\Psi\end{array}\right]K\left[\begin{array}[]{c}\Phi^{\top}\\\
\Psi^{\top}\end{array}\right]+\sigma^{2}I\right)^{-1}\mathrm{y}^{+}$ (125)
$\displaystyle\equiv-\frac{1}{\sigma^{2}}(\mathrm{y}^{+})^{\top}\left[\begin{array}[]{cc}\Phi&\Psi\end{array}\right]\left(\left[\begin{array}[]{c}\Phi^{\top}\\\
\Psi^{\top}\end{array}\right]\left[\begin{array}[]{cc}\Phi&\Psi\end{array}\right]+\sigma^{2}K^{-1}\right)^{-1}\left[\begin{array}[]{c}\Phi^{\top}\\\
\Psi^{\top}\end{array}\right]\mathrm{y}^{+}$ (132)
$\displaystyle=-\frac{N^{2}}{\sigma^{2}}\left[\begin{array}[]{cc}\hat{\theta}_{g,LS}^{\top}&\hat{\theta}_{f,LS}^{\top}\end{array}\right]\tilde{K}^{-1/2}\left(\tilde{K}^{1/2}\left[\begin{array}[]{c}\Phi^{\top}\\\
\Psi^{\top}\end{array}\right]\left[\begin{array}[]{cc}\Phi&\Psi\end{array}\right]\tilde{K}^{1/2}\right.$
(137) $\displaystyle\left.\hskip
11.38092pt+\sigma^{2}X^{-1}\right)^{-1}\tilde{K}^{-1/2}\left[\begin{array}[]{c}\hat{\theta}_{g,LS}\\\
\hat{\theta}_{f,LS}\end{array}\right]$ (140)
$\displaystyle=-\frac{N^{2}}{\sigma^{2}}\left[\begin{array}[]{cc}\hat{\theta}_{g,LS}^{\top}&\hat{\theta}_{f,LS}^{\top}\end{array}\right]\tilde{K}^{-1/2}\left(NI+\sigma^{2}X^{-1}\right)^{-1}\tilde{K}^{-1/2}\left[\begin{array}[]{c}\hat{\theta}_{g,LS}\\\
\hat{\theta}_{f,LS}\end{array}\right]$ (144)
$\displaystyle=-\frac{N^{2}}{\sigma^{2}}\sum_{\begin{subarray}{c}h\in I_{1}\\\
k\in I_{2}\end{subarray}}\left[\sum_{\begin{subarray}{c}j\in I_{1}\\\ l\in
I_{2}\end{subarray}}\frac{\|\hat{g}^{[hk,jl]}_{LS}\|_{P^{-1}}^{2}}{N+\sigma^{2}\frac{\lambda_{hj}+\gamma_{kl}}{\lambda_{hj}\gamma_{kl}}}+\sum_{i\in
I_{u}}\frac{\|\hat{f}^{[hk,i]}_{LS}\|_{R^{-1}}^{2}}{N+\sigma^{2}\frac{\pi_{h}+\omega_{ki}}{\pi_{h}\omega_{ki}}}\right].$

Thus,

$\displaystyle\ell(\mathrm{y}^{+},\xi)\equiv\sum_{\begin{subarray}{c}h\in
I_{1}\\\ k\in I_{2}\end{subarray}}\left[T\sum_{\begin{subarray}{c}j\in
I_{1}\\\ l\in
I_{2}\end{subarray}}\log\det\left(\frac{N}{\sigma^{2}}\frac{\lambda_{hj}\gamma_{kl}}{\lambda_{hj}+\gamma_{kl}}+1\right)\right.$
$\displaystyle\hskip 11.38092pt+T\sum_{\begin{subarray}{c}i\in
I_{u}\end{subarray}}\log\det\left(\frac{N}{\sigma^{2}}\frac{\pi_{h}\omega_{ki}}{\pi_{h}+\omega_{ki}}+1\right)$
$\displaystyle\left.-\frac{N^{2}}{\sigma^{2}}\sum_{\begin{subarray}{c}j\in
I_{1}\\\ l\in
I_{2}\end{subarray}}\frac{\|\hat{g}^{[hk,jl]}_{LS}\|_{P^{-1}}^{2}}{N+\sigma^{2}\frac{\lambda_{hj}+\gamma_{kl}}{\lambda_{hj}\gamma_{kl}}}-\frac{N^{2}}{\sigma^{2}}\sum_{i\in
I_{u}}\frac{\|\hat{f}^{[hk,i]}_{LS}\|_{R^{-1}}^{2}}{N+\sigma^{2}\frac{\pi_{h}+\omega_{ki}}{\pi_{h}\omega_{ki}}}\right].$

Next, we prove the first statement regarding $\lambda_{hj}$. The expression of
the partial derivative of $\ell$ with respect to $\lambda_{hj}$ depends on the
value of $\gamma_{kl}$’s. Given $\gamma_{kl}\geq 0$ for $k,l\in I_{2}$, let
$\tilde{J}$ denotes the subset of $I_{2}\times I_{2}$ containing the indexes
$(k,l)$ such that $\gamma_{kl}>0$. Then, we have:

$\displaystyle\frac{\partial\ell(\mathrm{y}^{+},\xi)}{\partial\lambda_{hj}}$
$\displaystyle=\sum_{\begin{subarray}{c}(k,l)\in\tilde{J}\end{subarray}}\left(\frac{\gamma_{kl}}{\lambda_{hj}+\gamma_{kl}}\right)^{2}\frac{T\left[\frac{\lambda_{hj}\gamma_{kl}}{\lambda_{hj}+\gamma_{kl}}+\frac{\sigma^{2}}{N}\right]-\|\hat{g}_{LS}^{[hk,jl]}\|^{2}_{P^{-1}}}{\left[\frac{\lambda_{hj}\gamma_{kl}}{\lambda_{hj}+\gamma_{kl}}+\frac{\sigma^{2}}{N}\right]^{2}},$

thus

$\displaystyle\left.\frac{\partial\ell(\mathrm{y}^{+},\xi)}{\partial\lambda_{hj}}\right|_{\begin{subarray}{c}\lambda_{hj}=0\end{subarray}}=\sum_{\begin{subarray}{c}(k,l)\in\tilde{J}\end{subarray}}\frac{\frac{T}{N}\sigma^{2}-\|\hat{g}_{LS}^{[hk,jl]}\|^{2}_{P^{-1}}}{\sigma^{4}/N^{2}}.$
(145)

If condition (93) holds, and in view of (145), then

$\displaystyle\left.\frac{\partial\ell(\mathrm{y}^{+},\xi)}{\partial\lambda_{hj}}\right|_{\lambda_{hj}=0}\geq
0.$ (146)

Notice that $\ell$ is continuous in the compact set $\Xi_{1}$ and thus $\ell$
admits a point of minimum in $\Xi_{1}$. Therefore, if (146) holds, then
$\Xi_{1}$ contains a local minimum for $\ell$. The remaining statements can be
proved in a similar way. ∎

######  Remark 11.

In the nonparametric case, i.e. the limit case $T\rightarrow\infty$, the above
result holds provided that the length of the data diverges in the same way.
More precisely, assume that $T\rightarrow\infty$, $N\rightarrow\infty$ and
$T/N\rightarrow 1$, then Proposition 10 holds where conditions (93)-(94) are
replaced by

$\displaystyle\sigma^{2}-\|\hat{g}_{LS}^{[hk,jl]}\|^{2}_{P^{-1}}\geq
0,\;\;\forall\,k,l\in I_{2}$
$\displaystyle\sigma^{2}-\|\hat{g}_{LS}^{[hk,jl]}\|^{2}_{P^{-1}}\geq
0,\;\;\forall\,h,j\in I_{1}$
$\displaystyle\sigma^{2}-\|\hat{f}_{LS}^{[hk,i]}\|^{2}_{R^{-1}}\geq
0,\;\;\forall\,k\in I_{2},\;i\in I_{u}$
$\displaystyle\sigma^{2}-\|\hat{f}_{LS}^{[hk,i]}\|^{2}_{R^{-1}}\geq
0,\;\;\forall\,h\in I_{1}.$

Let $\hat{G}_{LS}(z)$ and $\hat{F}_{LS}(z)$ be the transfer matrices
corresponding to the least squares estimates in (92). Then, the corresponding
model is

$\displaystyle y(t)=\hat{G}_{LS}(z)y(t)+\hat{F}_{LS}(z)+e(t)$

The first statement of Remark 11 means that if the $\ell_{2}$ norms of
$\hat{g}_{LS}^{[hk,jl]}$’s describing the conditional Granger causality
relations from $y^{\star}_{j}$ to $y^{\star}_{h}$ do not dominate the variance
of the noise process, then $\hat{\lambda}_{hj}=0$111Here, we implicitly assume
that the local minimum in $\Xi_{1}$ is chosen as optimal solution. and thus
the posterior mean of $g^{[hk,jl]}$’s belonging to the block $(h,j)$ are null
functions. Similar reasonings hold for the other statements.

##  6 Numerical experiments

In this section we test the performance of the proposed estimator.

###  6.1 Synthetic examples

We consider five Monte Carlo studies of 100 runs. For each run we generate an
ARMAX model of the form

$\displaystyle y(t)=G_{\mathrm{TR}}(z)y(t)+F_{\mathrm{TR}}(z)u(t)+e(t)$ (147)

as follows:

  * •

the rational transfer matrices $G_{\mathrm{TR}}(z)$ and $F_{\mathrm{TR}}(z)$
are strictly casual and generated through the MATLAB function drmodel.m where
the order of the system is equal to $20$; moreover, their poles are restricted
to have absolute value less than 0.95.

  * •

some entries of $G_{\mathrm{TR}}(z)$ and $F_{\mathrm{TR}}(z)$ are set equal to
the null function in such a way that they satisfy condition (13); matrices
$E_{1}$, $E_{2}$, $A_{1}$ and $A_{2}$ are randomly generated in such a way
such that the fraction of nonnull entries is equal to 0.6.

For each run we generate a dataset of length $N=500$ if not otherwise
specified; the input $u$ has uncorrelated components which are generated
through the MATLAB function idinput.m as a realization from a random Gaussian
noise with band [0, 0.4].

We compare the following estimators:

  * •

K: this is the estimator based on model (12) where the transfer matrices are
modeled as zero mean Gaussian processes with kernels (85);

  * •

S: this is the estimator based on model (12) where the transfer matrices are
modeled as zero mean Gaussian processes with kernels (87);

  * •

SS: this is the estimator based on model (12) where the transfer matrices are
modeled as zero mean Gaussian processes with kernels $\lambda
I_{p_{1}^{2}p_{2}^{2}}\otimes P$ and $\pi I_{p_{1}p_{2}m}\otimes R$ where the
hyperparameters $\lambda,\pi\geq 0$ represent the scale factors.

The impulse responses are truncated with length $T=50$, i.e. we make the
approximation
$G(z)=\sum_{t=1}^{\infty}G_{t}z^{-t}\approx\sum_{t=1}^{T}G_{t}z^{-t}$. Such a
truncation does not affect the estimates because both $P$ and $R$ force the
estimated impulse responses to decay exponentially, i.e. it does not involve
any kind of trade-off between bias and variance (Pillonetto & De Nicolao,
2010). In practice, it is only important to take the practical length very
large.

The performance of these estimators is assessed through the following indexes:

  * •

Average impulse response fit

$\displaystyle\mathrm{AIRF}=\frac{fit(G_{\mathrm{TR}},\hat{G})}{2}+\frac{fit(F_{\mathrm{TR}},\hat{F})}{2}$

where

$\displaystyle
fit(G_{\mathrm{TR}},\hat{G})=100\left(1-\sqrt{\frac{\frac{1}{T}\sum_{t=1}^{T}\|G_{\mathrm{TR},t}-\hat{G}_{t}\|^{2}}{\frac{1}{T}\sum_{t=1}^{T}\|G_{\mathrm{TR},t}-\bar{G}_{\mathrm{TR}}\|^{2}}}\right)$
$\displaystyle\bar{G}_{\mathrm{TR}}=\frac{1}{T}\sum_{t=1}^{T}G_{\mathrm{TR},t}$

and $\hat{G}(z)$, $\hat{F}(z)$ are the estimated transfer matrices;

  * •

Fraction of misspecified edges

$\displaystyle\mathrm{ERR}=\frac{\|E_{1}\otimes
E_{2}-\hat{E}_{1}\otimes\hat{E}_{2}\|_{0}}{2p_{1}^{2}p_{2}^{2}}+\frac{\|A_{1}\otimes
A_{2}-\hat{A}_{1}\otimes\hat{A}_{2}\|_{0}}{2p_{1}p_{2}m}$

where $\mathrm{supp}(\hat{G}(z))=\hat{E}_{1}\otimes\hat{E}_{2}$,
$\mathrm{supp}(\hat{F}(z))=\hat{A}_{1}\otimes\hat{A}_{2}$, while $E_{1}$,
$E_{2}$, $A_{1}$ and $A_{2}$ correspond to the support of $G_{\mathrm{TR}}(z)$
and $F_{\mathrm{TR}}(z)$, respectively.

First Monte Carlo study. We consider the case in which model (147) is without
input $u$, i.e. we have a dynamic network corresponding to an ARMA stochastic
process. In this case we only have to estimate $G(z)$, accordingly
$\mathrm{AIRF}=fit(G_{\mathrm{TR}},\hat{G})$ and $\mathrm{ERR}=\|E_{1}\otimes
E_{2}-\hat{E}_{1}\otimes\hat{E}_{2}\|_{0}/(p_{1}^{2}p_{2}^{2})$. Besides the
previous estimators, we also consider:

  * •

PEM: this is the classical PEM approach that uses BIC for model order
selection. The orders of the polynomials in the ARMA model ranges from 1 to 15
and are not allowed to be different from each other since this would lead to a
combinatorial explosion of the number of competitive models.

  * •

PEM+OR: this is the same as PEM with an additional oracle knowing which
impulse responses are null functions.

The length of the dataset is $N=3900$. It is worth noting the latter has been
chosen so large in such a way that PEM with polynomial orders equal to 15 has
enough data to estimate all the model parameters whose amount is 3848.

Figure 2: First Monte Carlo study “ARMA case” with $p_{1}=4$, $p_{2}=4$. Left
panel. Average impulse response fit. Right panel. Fraction of misspecified
edges.

Figure 2 shows the performance of the estimators. K outperforms all the other
estimators in terms of $\mathrm{AIRF}$ while both PEM+OR, even though it owns
the correct network topology, and PEM are the worst. Finally, K is better than
S in terms of $\mathrm{ERR}$. Notice that we did not plot $\mathrm{ERR}$ for
PEM and SS because they always provide a full network. In particular, in
regard to SS, the entries of $G(z)$ share the same prior (and thus the same
scale factor $\lambda$ whose estimated value is always strictly positive). We
also considered Monte Carlo studies where $p_{1}$ and $p_{2}$ are different:
the results, in terms of performance, do not change.

Second Monte Carlo study. The output dimensions are $p_{1}=p_{2}=3$, while the
input dimension is $m=2$.

Figure 3: Second Monte Carlo study with $p_{1}=3$, $p_{2}=3$ and $m=2$. Left
panel. Average impulse response fit. Right panel. Fraction of misspecified
edges.

Figure 3 shows the performance of the three estimators: K outperforms all the
other estimators in terms of $\mathrm{AIRF}$, while SS is the worst. Finally,
K is better than S in terms of $\mathrm{ERR}$.

Third Monte Carlo study. We consider the case in which $p_{1}$ and $p_{2}$ are
different; more precisely, we have $p_{1}=3$, $p_{2}=4$ and $m=2$.

Figure 4: Third Monte Carlo study with $p_{1}=3$, $p_{2}=4$ and $m=2$. Left
panel. Average impulse response fit. Right panel. Fraction of misspecified
edges.

Figure 4 shows the obtained results. Also in this case K is the best estimator
both in terms of $\mathrm{AIRF}$ and $\mathrm{ERR}$, while SS is the worst.

Fourth Monte Carlo study. We consider $p_{1}=4$, $p_{2}=3$ and $m=2$; this
Monte Carlo study differs from the last one because the output dimensions are
swapped.

Figure 5: Fourth Monte Carlo study with $p_{1}=4$, $p_{2}=3$ and $m=2$. Left
panel. Average impulse response fit. Right panel. Fraction of misspecified
edges.

The situation substantially does not change, see Figure 5.

Fifth Monte Carlo study. We consider the case in which model (147) is without
input $u$ and corresponds to a hierarchical dynamic network (i.e.
$E_{1}=E_{2}$). We set $p_{1}=p_{2}=4$. Besides the previous estimators, we
also consider:

  * •

H: this is the estimator based on model (12), without input, where the
transfer matrix in $G(z)$ is modeled as zero mean Gaussian process with kernel
(88).

$\mathrm{AIRF}$ and $\mathrm{ERR}$ are defined as in the first Monte Carlo
study.

Figure 6: Fifth Monte Carlo study “ARMA case - hierarchical network” with
$p_{1}=4$, $p_{2}=4$. Left panel. Average impulse response fit. Right panel.
Fraction of misspecified edges.

In Figure 6 the performance of the estimators is depicted. S and SS are worse
than H and K. The latter perform in a similar way in terms of $\mathrm{AIRF}$.
H, however, outperforms K in terms of $\mathrm{ERR}$.

######  Remark 12.

The negative log-marginal likelihood function is not convex with respect to
the hyperparameters even in the simplest case where the hyperparameters are
only the scale factors, see Chen et al. (2014). Accordingly, the kernel-based
PEM methods could suffer about local minima in the estimation of the
hyperparameters. In particular, both S and K-H have the same issue. The main
difference is the number of hyperparameters to estimate: K-H have to estimate
very less hyperparameters than S. More precisely, in K we have to estimate
$p_{1}^{2}+p_{2}^{2}+p_{1}+p_{2}m$ ($p_{1}^{2}+p_{2}^{2}$ in the case of
systems without input) hyperparameters, in H we have to estimate $p_{1}^{2}$
hyperparameters (only applicable in the case of systems without input and
$p_{1}=p_{2}$), while in S we have to estimate
$p_{1}^{2}p_{2}^{2}+p_{1}p_{2}m$ ($p_{1}^{2}p_{2}^{2}$ in the case of systems
without input) hyperparameters. Table 1 summarizes the number of
hyperparameters needed in the previous Monte Carlo studies.

Table 1: Number of hyperparameters to estimate. M.C. study | S | K | H
---|---|---|---
First | 99 | 27 | -
Second | 168 | 36 | -
Third | 168 | 35 | -
Fourth | 256 | 32 | -
Fifth | 256 | 32 | 16

Therefore, K and H alleviate the aforementioned issue since the search space
for optimizing the marginal likelihood is small in respect to the one of S.

###  6.2 Learning a bike sharing system

Bike sharing systems aim to make automatic bike rentals, in particular the
rental and the return back. We consider the number of users and some variables
describing the weather, that is wind speed, humidity and apparent temperature.
These quantities have been collected in the period 1 January 2011 - 14 May
2012 (500 days in total) from the Capital Bike Sharing (CBS) system at
Washington, D.C., USA. The sampling time is equal to 1 hour. For more details
see Fanaee-T & Gama (2014). These data describe a four dimensional stochastic
process $x(s)=[\,x_{1}(s)\;x_{2}(s)\;x_{3}(s)\;x_{4}(s)\,]^{\top}$, where
$x_{1}(s)$ is the number of users, $x_{2}(s)$ is the wind speed, $x_{3}(s)$ is
the humidity and $x_{4}(s)$ is the apparent temperature. Process $x(s)$ is
non-stationary during a day: for instance, we expect that the number of users
during the night is smaller than the ones during the day. Accordingly, we
consider the process

$\displaystyle y(t)=[\,x(24(t-1)+1)^{\top}\;\ldots\;\,x(24t)^{\top}\,]^{\top}$

taking values in $\mathbb{R}^{96}$ and the corresponding sampling time is 1
day. Thus, the corresponding dataset has length $N=500$. After detrending it,
we model $y$ as the dynamic spatio-temporal model of Section 2.1 where
$p_{1}=24$, $p_{2}=4$. Then, we estimate $G(z)$ by K and the regression
matrices are the ones in Remark 2. The estimated network describing the
conditional Granger causality relations among the number of users, wind speed,
humidity and apparent temperature is depicted in Figure 7.

Figure 7: Estimated network describing conditional Granger causality relations
among: the number of users (U), wind speed (W), humidity (H) and apparent
temperature (T).

As expected, the number of users does not conditionally Granger causes the
weather variables. Moreover, the number of users is conditionally Granger
caused by the apparent temperature.

Figure 8: Weighted matrix describing the estimated Granger causality relations
among hours.

Figure 8 shows the weighted matrix describing the conditional Granger
causality relations among hours. More precisely, the entry in position $(h,j)$
is the $\ell_{1}$ norm of the $4\times 4$ impulse response, formed by
$\hat{g}^{[hk,jl]}$ with $k,l\in\\{1,2,3,4\\}$, which describes the strength
of the conditional Granger causality relation from hour $j$ to hour $h$. We
can notice that all the hours are mainly conditionally Granger caused by the
previous hour.

##  7 Conclusions

We have considered the problem to estimate a dynamic network describing
conditional Granger causality relations and having a Kronecker structure. We
have proposed a kernel-based PEM method for estimating such a model from data.
The kernel functions have been derived from the maximum entropy principle.
Numerical simulation showed that the proposed estimator outperforms the one
estimating dynamic networks without structure. Finally, we have used the
proposed method to learn the conditional Granger causality relations in a bike
sharing system.

## References

  * Øksendal (1998) Øksendal, B. (1998).  Stochastic differential equations: an introduction with applications.  Universitext.  Berlin, Heidelberg, New York: Springer.
  * Akaike (1974) Akaike, H. (1974).  A new look at the statistical model identification.  IEEE Transactions on Automatic Control, 19, 716–723.
  * Alpago et al. (2018) Alpago, D., Zorzi, M., & Ferrante, A. (2018).  Identification of sparse reciprocal graphical models.  IEEE Control Systems Letters, 2, 659–664.
  * Aravkin et al. (2014) Aravkin, A., Burke, J. V., Chiuso, A., & Pillonetto, G. (2014).  Convex vs non-convex estimators for regression and sparse estimation: the mean squared error properties of ARD and GLasso.  The Journal of Machine Learning Research, 15, 217–252.
  * Aronszajn (1950) Aronszajn, N. (1950).  Theory of reproducing kernels.  Trans. Amer. Math. Soc., 68, 337–404.
  * Avventi et al. (2013) Avventi, E., Lindquist, A. G., & Wahlberg, B. (2013).  ARMA identification of graphical models.  IEEE Transactions on Automatic Control, 58, 1167–1178.
  * Bijma et al. (2005) Bijma, F., De Munck, J. C., & Heethaar, R. M. (2005).  The spatiotemporal MEG covariance matrix modeled as a sum of kronecker products.  NeuroImage, 27, 402–415.
  * Bonettini et al. (2015) Bonettini, S., Chiuso, A., & Prato, M. (2015).  A scaled gradient projection methods for Bayesian learning in dynamical systems.  SIAM Journal on Scientific Computing, 37, A1297–A1318.
  * Chen et al. (2014) Chen, T., Andersen, M., Ljung, L., Chiuso, A., & Pillonetto, G. (2014).  System identification via sparse multiple kernel-based regularization using sequential convex optimization techniques.  IEEE Transactions on Automatic Control, 59, 2933–2945.
  * Chen & Ljung (2013) Chen, T., & Ljung, L. (2013).  Implementation of algorithms for tuning parameters in regularized least squares problems in system identification.  Automatica, 49, 2213–2220.
  * Chen et al. (2012) Chen, T., Ohlsson, H., & Ljung, L. (2012).  On the estimation of transfer functions, regularizations and gaussian processes-revisited.  Automatica, 48, 1525–1535.
  * Chiuso (2016) Chiuso, A. (2016).  Regularization and Bayesian learning in dynamical systems: Past, present and future.  Annual Reviews in Control, 41, 24–38.
  * Chiuso & Pillonetto (2012) Chiuso, A., & Pillonetto, G. (2012).  A Bayesian approach to sparse dynamic network identification.  Automatica, 48, 1553–1565.
  * Ciccone et al. (2020) Ciccone, V., Ferrante, A., & Zorzi, M. (2020).  Learning latent variable dynamic graphical models by confidence sets selection.  IEEE Transactions on Automatic Control, 65, 5130–5143.
  * Cover & Thomas (1991) Cover, T., & Thomas, J. (1991).  Information Theory.  New York: Wiley.
  * Doddi et al. (2020) Doddi, H., Deka, D., Talukdar, S., & Salapaka, M. (2020).  Estimating linear dynamical networks of cyclostationary processes.  arXiv preprint arXiv:2009.12667, .
  * Fanaee-T & Gama (2014) Fanaee-T, H., & Gama, J. (2014).  Event labeling combining ensemble detectors and background knowledge.  Progress in Artificial Intelligence, 2, 113–127.
  * Goodwin et al. (1992) Goodwin, G., Gevers, M., & Ninness, B. (1992).  Quantifying the error in estimated transfer functions with application to model order selection.  IEEE Transactions on Automatic Control, 37, 913–928.
  * Granger (1969) Granger, C. (1969).  Investigating causal relations by econometric models and cross-spectral methods.  Econometrica, 37, 424–438.
  * Jaynes (1982) Jaynes, E. (1982).  On the rationale of maximum-entropy methods.  Proceedings of the IEEE, 70, 939–952.
  * Jorgesen et al. (2011) Jorgesen, P., Kornelson, K., & Shuman, K. (2011).  Iterated Function Systems, Moments, and Transformations of Infinite Matrices.  American Mathematical Society.
  * Lauritzen (1996) Lauritzen, S. (1996).  Graphical Models.  Oxford: Oxford University Press.
  * Leskovec (2009) Leskovec, J. (2009).  Networks, communities and Kronecker products.  In Proceedings of the 1st ACM international workshop on Complex networks meet information & knowledge management (pp. 1–2).
  * Leskovec et al. (2010) Leskovec, J., Chakrabarti, D., Kleinberg, J., Faloutsos, C., & Ghahramani, Z. (2010).  Kronecker graphs: An approach to modeling networks.  Journal of Machine Learning Research, 11, 985–1042.
  * Leskovec & Faloutsos (2007) Leskovec, J., & Faloutsos, C. (2007).  Scalable modeling of real graphs using kronecker multiplication.  In Proceedings of the 24th international conference on Machine learning (pp. 497–504).
  * Ljung (1999) Ljung, L. (1999).  System Identification: Theory for the User.  New Jersey: Prentice Hall.
  * Ljung et al. (2020) Ljung, L., Chen, T., & Mu, B. (2020).  A shift in paradigm for system identification.  International Journal of Control, 93, 173–180.
  * MacKay (1994) MacKay, D. J. (1994).  Bayesian nonlinear modeling for the prediction competition.  ASHRAE transactions, 100, 1053–1062.
  * Materassi & Innocenti (2010) Materassi, D., & Innocenti, G. (2010).  Topological identification in networks of dynamical systems.  IEEE Transactions on Automatic Control, 55, 1860–1871.
  * Materassi & Salapaka (2020) Materassi, D., & Salapaka, M. V. (2020).  Signal selection for estimation and identification in networks of dynamic systems: A graphical model approach.  IEEE Transactions on Automatic Control, 65, 4138–4153.
  * Mu et al. (2018) Mu, B., Chen, T., & Ljung, L. (2018).  On asymptotic properties of hyperparameter estimators for kernel-based regularization methods.  Automatica, 94, 381–395.
  * Pillonetto et al. (2011) Pillonetto, G., Chiuso, A., & De Nicolao, G. (2011).  Prediction error identification of linear systems: A nonparametric gaussian regression approach.  Automatica, 47, 291–305.
  * Pillonetto & De Nicolao (2010) Pillonetto, G., & De Nicolao, G. (2010).  A new kernel-based approach for linear system identification.  Automatica, 46, 81–93.
  * Pillonetto et al. (2014) Pillonetto, G., Dinuzzo, F., Chen, T., De Nicolao, G., & Ljung, L. (2014).  Kernel methods in system identification, machine learning and function estimation: A survey.  Automatica, 50, 657–682.
  * Ramaswamy et al. (2021) Ramaswamy, K. R., Bottegal, G., & Van den Hof, P. M. (2021).  Learning linear modules in a dynamic network using regularized kernel-based methods.  Automatica, 129, 109591.
  * Ramaswamy & Vandenhof (2020) Ramaswamy, K. R., & Vandenhof, P. M. J. (2020).  A local direct method for module identification in dynamic networks with correlated noise.  IEEE Transactions on Automatic Control, (pp. 1–1).
  * Rasmussen & Williams (2006) Rasmussen, C., & Williams, C. (2006).  Gaussian Processes for Machine Learning.  The MIT Press.
  * Sepehr & Materassi (2019) Sepehr, F., & Materassi, D. (2019).  Blind learning of tree network topologies in the presence of hidden nodes.  IEEE Transactions on Automatic Control, 65, 1014–1028.
  * Sinquin & Verhaegen (2019) Sinquin, B., & Verhaegen, M. (2019).  Quarks: Identification of large-scale kronecker vector-autoregressive models.  IEEE Transactions on Automatic Control, 64, 448–463.
  * Söderström & Stoica (1989) Söderström, T., & Stoica, P. (1989).  System Identification.  Hemel Hempstead, UK: Prentice-Hall International.
  * Tsiligkaridis & Hero (2013) Tsiligkaridis, T., & Hero, A. O. (2013).  Covariance estimation in high dimensions via Kronecker product expansions.  IEEE Transactions on Signal Processing, 61, 5347–5360.
  * Tsiligkaridis et al. (2013) Tsiligkaridis, T., Hero III, A. O., & Zhou, S. (2013).  On convergence of Kronecker graphical lasso algorithms.  IEEE transactions on signal processing, 61, 1743–1755.
  * Veedu et al. (2020) Veedu, M. S., Doddi, H., & Salapaka, M. V. (2020).  Topology learning of linear dynamical systems with latent nodes using matrix decomposition.
  * Yu et al. (2009) Yu, K., Lafferty, J., Zhu, S., & Gong, Y. (2009).  Large-scale collaborative prediction using a nonparametric random effects model.  In Proceedings of the 26th Annual International Conference on Machine Learning (pp. 1185–1192).
  * Zorzi (2018) Zorzi, M. (2018).  Learning quasi-Kronecker product graphical models.  In IEEE Conference on Decision and Control (pp. 1427–1432).
  * Zorzi (2019) Zorzi, M. (2019).  Graphical model selection for a particular class of continuous-time processes.  Kybernetika, 55, 782–801.
  * Zorzi (2020a) Zorzi, M. (2020a).  Autoregressive identification of kronecker graphical models.  Automatica, 119, 109053.
  * Zorzi (2020b) Zorzi, M. (2020b).  A new kernel-based approach for spectral estimation.  In European Control Conference (ECC) (pp. 534–539).
  * Zorzi & Chiuso (2015) Zorzi, M., & Chiuso, A. (2015).  A Bayesian approach to sparse plus low rank network identification.  In Proceedings of the IEEE Conference on Decision and Control (pp. 7386–7391).  Osaka.
  * Zorzi & Chiuso (2017) Zorzi, M., & Chiuso, A. (2017).  Sparse plus low rank network identification: A nonparametric approach.  Automatica, 76, 355–366.
  * Zorzi & Chiuso (2018) Zorzi, M., & Chiuso, A. (2018).  The harmonic analysis of kernel functions.  Automatica, 94, 125–137.
  * Zorzi & Sepulchre (2016) Zorzi, M., & Sepulchre, R. (2016).  AR identification of latent-variable graphical models.  IEEE Transactions on Automatic Control, 61, 2327–2340.