# On the Spatial Pattern of Input-Output Metrics for a Network Synchronization
Process

Subir Sarker and Sandip Roy The authors acknowledge support from ARPA-E,
Pacific Northwest National Laboratory, and National Science Foundation.The
authors are with the School of Electrical Engineering and Computer Science,
Washington State University, Pullman, WA 99163, USA.
<EMAIL_ADDRESS>
###### Abstract

A graph-theoretic analysis is undertaken for a compendium of input-output
(transfer) metrics of a standard discrete-time linear synchronization model,
including $l_{p}$ gains, frequency responses, frequency-band energy, and
Markov parameters. We show that these transfer metrics exhibit a spatial
degradation, such that they are monotonically nonincreasing along vertex
cutsets away from an exogenous input. We use this spatial analysis to
characterize signal-to-noise ratios (SNRs) in diffusive networks driven by
process noise, and to develop a notion of propagation stability for dynamical
networks. Finally, the formal results are illustrated through an example.

##  I Introduction

The literature on dynamical networks has predominantly focused on emergence
and control of global internal properties such as synchronization or consensus
[1, 2]. Recently, several studies considered input-to-output or transfer
characteristics of dynamical networks. In particular, input-output analyses
have been motivated by questions regarding the disturbance responses of
networks as they are scaled [3, 4, 5]; these efforts generalize notions of
string stability [6] toward general network structures. In a separate track,
input-to-output analyses have also been motivated by controller design needs
in infrastructure networks with sparse actuation and measurement capabilities
[7, 8, 9, 10]. Additionally, questions related to security and estimability of
network processes have motivated input-output analyses [11, 12].

Prior input-to-output analyses of dynamical networks have been concerned with
particular metrics arising from application contexts, such as $l_{\infty}$
gains or transfer function zeros, and have primarily been focused on tying the
input-output metrics to global graph properties (e.g. [10]). Relative to these
efforts, our main contribution in this short article is to demonstrate that a
compendium of input-output metrics exhibit a special spatial pattern, for a
canonical network synchronization/consensus model. In particular, we study a
standard discrete-time model for consensus/synchronization among nodes with
scalar states [13]. For this model, we consider several input-to-output
metrics between node pairs, including $l_{p}$ gains, frequency responses,
frequency-band energy, and Markov parameters. The main outcome of our analyses
is to show that these metrics satisfy a decrescence property, wherein they are
nonincreasing in value across cutsets away from an input node. Additionally,
two implications of the spatial pattern analysis are developed in brief
vignettes. First, the analysis is used to give some insights into signal-to-
noise ratios in network measurements, which bear on estimation/detection and
sensor placement. Second, the analysis is used to define a notion of input-
output or propagation stability for networks, and to verify propagation
stability for the canonical model regardless of model scale (size).

##  II Model and Goals

A dynamical network with $n$ agents or nodes is considered. Each node
$i=1,\ldots,n$ has a scalar state $x_{i}(k)$ which evolves in discrete time
($k=0,1,2,\ldots$), via linear diffusive interactions with other nodes. Here,
our primary interest is in characterizing transfer characteristics in the
network, hence we also model an exogenous input $u(k)$ applied at one source
node $s$ and an output $y(k)$ taken at a target node $i$. In particular, the
following single-input single-output dynamics is considered:

$\displaystyle X(k+1)=AX(k)+e_{s}u(k)$ (1) $\displaystyle y(k)=e_{i}^{T}X(k)$

where $X(k)=\begin{bmatrix}x_{1}[k]\\\ \vdots\\\ x_{n}[k]\end{bmatrix}$,
$A=[a_{ij}]$ is assumed to be row stochastic or substochastic (i.e.
$a_{ij}\geq 0$ and $\sum_{j}a_{ij}\leq 1$), and the notation $e_{b}$ is used
for a $0--1$ indicator vector whose $b$th entry is equal to $1$. The zero-
state response of the system (1) is of interest in this study. The model (1)
is a standard model for network synchronization or consensus (for stochastic
$A$), and also is representative of other diffusive processes in networks.

A digraph $\Gamma$ is used to represent pairwise interactions among the nodes
in the network. Specifically, the graph $\Gamma$ is defined to have $n$
vertices labeled $1,\ldots,n$, which correspond with the $n$ network nodes. A
directed edge is drawn from vertex $j$ to vertex $l$ if $a_{lj}>0$, which
indicates that the state of vertex $j$ directly influences that of vertex $l$.
The vertices $s$ and $i$ are referred to as the source and target vertices,
respectively. To simplify the development, we assume throughout that $\Gamma$
is strongly connected, i.e., there is a directed path from each vertex to
every other vertex. The strongly-connected case is of primary interest when
considering transfer characteristics; although the results can be simply
generalized for non-strongly connected cases, these details are omitted.

Our primary aim in this study is to compare input-output metrics for the
network model (1) for different target vertex locations $i$ in $\Gamma$, so as
to characterize the spatial pattern of the input-output response. The specific
metrics that we consider for the system (1) for a particular target location
$i$ are:

  * •

The $l_{p}$ gain $G_{p}(i,k_{f})$ over a time horizon $k_{f}$, defined as the
maximum of the quantity
$\left[\sum_{k=0}^{k_{f}}|y(k)|^{p}\right]^{\frac{1}{p}}$ over all inputs
$u(0),\ldots,u(k_{f}-1)$ such that
$\left[\sum_{k=0}^{k_{f}}|u(k)|^{p}\right]^{\frac{1}{p}}=1$.

  * •

The frequency response $H_{i}(e^{j\Omega})$, where $\Omega$ is the frequency
of the (discrete time) input and $H_{i}()$ is the system’s transfer function.
Additionally, the responses to arbitrary periodic inputs, and the signal
content in a frequency band, are also characterized.

  * •

The Markov parameters $M_{i}(k)={\bf e}_{i}^{T}A^{k}{\bf e}_{s}$ for
$k=0,1,2,\ldots$.

These metrics together form the basis for the external stability analysis of
linear systems, and modulate estimation and feedback controller design.

Conceptually, the diffusive structure of the dynamics (1) suggests that inputs
should have a localized sphere of influence, and hence the input-output
metrics should exhibit a spatial degradation with distance from the source.
Our goal is to provide a formal characterization of this spatial falloff, and
to develop implications of this spatial analysis.

##  III Main Results: Spatial Pattern Analysis

Graph theoretic analyses are developed for the input-output metrics defined
for the system (1). The main results show that the metrics ($l_{p}$ gain,
frequency response, Markov parameters) for different target locations falloff
monotonically along graph cutsets away from the source location. To formalize
these notions, it is convenient to define the notion of a separating cutset
for a graph. To do so, let us consider a network graph $\Gamma$ with a source
vertex $s$ and a particular target vertex $i=q^{*}$. A set of vertices
$C=\\{c(1),c(2),...,c(m)\\}$ (which does not contain $s$ and $q^{*}$) is said
to be a separating cutset, if any directed path from $s$ to $q^{*}$ passes
through at least one vertex in $C$. The concept is illustrated in Figure 1. We
also find it convenient to use the notation $Z$ for the partition of $\Gamma$
containing the target vertex, upon removal of the separating cutset. We refer
to $Z$ as the target partition.

Figure 1: Illustration of a vertex cutset.

The spatial pattern analyses of the input-output metrics depend on a
reformulation of system (1), wherein the target node’s state dynamics are
expressed in terms of the states of nodes corresponding to a separating
cutset, rather than directly in terms of the input. To develop this
reformulation, let us define the cutset state vector $X_{c}(k)$ as containing
the time-$k$ states of the network nodes corresponding to the vertices in $C$.
Likewise, we define the target partition state vector $X_{z}(k)$ as containing
the time-$k$ states corresponding to the vertices in $Z$ (including $q^{*}$).
The dynamical update of the target partition state vector can be expressed in
terms of only the current cutset state and target partition state:

$X_{z}(k+1)=A_{z}X_{z}(k)+BX_{c}(k)$ (2)

where $A_{z}$ is the principal submatrix of $A$ formed by selecting the rows
and columns of $A$ indicated in $Z$; and $B$ is a submatrix of $A$ formed by
selecting the rows indicated in $Z$ and the columns indicated in $C$.

The main graph-theoretic results build on the following characterization of
the reformulated state dynamics (2):

######  Lemma 1.

Assume that the network model (1) is initially relaxed and an exogenous input
$u(k)$ is applied at the source node $s$. Consider any separating cutset $C$
and corresponding target partition $Z$. The target partition state vector over
the interval $k=0,1,\ldots,k_{f}$ can be expressed in terms of the cutset
state vector as follows:

${\small\begin{bmatrix}X_{z}(0)\\\ \vdots\\\
X_{z}(k_{f})\end{bmatrix}}=Q_{k_{f}}{\small\begin{bmatrix}X_{c}(k_{f})\\\
\vdots\\\ X_{z}(0)\end{bmatrix}},$ (3)

where

$Q_{k}={\small\begin{bmatrix}0&\dots&0&0\\\ 0&\dots&0&B\\\ 0&\dots&B&A_{z}B\\\
\vdots&\ddots&\vdots&\vdots\\\
0&\dots&A_{z}^{k-2}B&A_{z}^{k-1}B\end{bmatrix}}.$ (4)

Further, the matrix $Q_{k_{f}}$ has nonnegative entries, and the sum of the
entries in each row is at most $1$.

###### Proof.

The target partition state vector can be computed in terms of the cutset state
vector by solving (2). This gives:

$X_{z}(k)=\sum_{j=0}^{k-1}A_{z}^{j}BX_{c}(k-j-1)$ (5)

Assembling the responses for $k=0,\ldots,k_{f}$ then immediately yields
Equations (3) and (4) in the lemma statement.

The entries in $Q_{k_{f}}$ are nonnegative, as $A_{z}$ and $B$ are
nonnegative. To characterize the row sums for $Q_{k_{f}}$, let us first
consider the last block of rows $\hat{Q}_{k_{f}}=[0\quad B\ldots
A_{z}^{k_{f}-2}B\quad A_{z}^{k_{f}-1}B]$. We characterize the row sums of
$\hat{Q}_{k_{f}}$ by considering powers of the following matrix $F$:

$F={\small\begin{bmatrix}A_{z}&B\\\ \mathbf{0}&I_{m}\end{bmatrix}}$

The $k_{f}$th power of matrix $F$ is given by:

$F^{k_{f}}={\small\begin{bmatrix}A^{k_{f}}_{z}&\sum_{k=0}^{k_{f}-1}A^{k}_{z}B\\\
\mathbf{0}&I_{m}\end{bmatrix}}$

Since $[A_{z}\quad B]$ is a submatrix of $A$, it is immediate that $F$ is a
stochastic or substochastic matrix. Consequently, the row sums of $F^{k_{f}}$
are also at most $1$. From the expression for $F^{k_{f}}$, it follows that the
matrix $\sum_{k=0}^{k=k_{f}-1}A_{z}^{k}B$ has row sums of at most $1$.
However, the row sums of $\sum_{k=0}^{k=k_{f}-1}A_{z}^{k}B$ and the row sums
of the matrix $\hat{Q}_{k_{f}}$ are identical. Hence, the row sums of
$\hat{Q}_{k_{f}}$ are upper bounded by $1$. Further, the sum of each row of
the matrix $Q_{k_{f}}$ is upper bounded by one of the row sums of
$\hat{Q}_{k_{f}}$. The result thus follows.

∎

Our first main result is a graph-theoretic characterization of the $l_{p}$
gains of the input-output model (1):

######  Theorem 1.

Consider the $l_{p}$ gain of the network model (1) for a particular target
node $i=q^{*}$, i.e. $G_{p}(q^{*},k_{f})$. Also consider the $l_{p}$ gains
when the target node is alternately on a separating cutset
$C=\\{c(1),c(2),...,c(m)\\}$, i.e. $G_{p}(c(i),k_{f})$. For any time horizon
$k_{f}$, it holds that $G_{p}(q^{*},k_{f})\leq G_{p}(c(i),k_{f})$ for some
$i=1,2,\ldots,m$.

###### Proof.

Theorem 1 is verified for finite $p$ (i.e., $1\leq p<\infty$), and then for
$p\to\infty$.

######  Case 1.

$1\leq p<\infty$

The state of any node $q$ other than the source $s$ evolves as:

$x_{q}(k+1)=a_{qq}x_{q}(k)+\sum_{j\in N(q)}a_{qj}x_{j}(k)$ (6)

where $N(q)$ contains the upstream neighbors of vertex $q$ in $\Gamma$ (the
vertices with directed edges to $q$). Notice that $a_{qj}>0$ for $j\in{\cal
N}(q)$ and $\sum_{j=1}^{n}a_{qj}\leq 1$.

For any $1\leq p<\infty$, the function $f(z)=|z|^{p}$ is convex over all $z$.
Considering $z=x_{q}(k+1)$, it follows immediately that

$\begin{split}|x_{q}(k+1)|^{p}\leq a_{qq}|x_{q}(k)|^{p}+\sum_{j\in
N(q)}a_{qj}|x_{j}(k)|^{p}\end{split}$ (7)

Summing the inequality for $k=0,\ldots,k_{f}$ yields:

$\sum_{k=0}^{k_{f}}|x_{q}(k+1)|^{p}\leq
a_{qq}\sum_{k=0}^{k_{f}}|x_{q}(k)|^{p}+\sum_{j\in
N(q)}a_{qj}\sum_{k=0}^{k_{f}}|x_{j}(k)|^{p}$ (8)

Rewriting the leftmost sum and using $x_{q}(0)=0$, we get:

$|x_{q}(k_{f}+1)|^{p}+\sum_{k=0}^{k_{f}}|x_{q}(k)|^{p}\leq
a_{qq}\sum_{k=0}^{k_{f}}|x_{q}(k)|^{p}+\\\ \sum_{j\in
N(q)}a_{qj}\sum_{k=0}^{k_{f}}|x_{j}(k)|^{p}$ (9)

As $|x_{q}(k_{f}+1)|^{p}$ is nonnegative, we then obtain:

$\sum_{k=0}^{k_{f}}|x_{q}(k)|^{p}\leq\sum_{j\in
N(q)}\frac{a_{qj}}{1-a_{qq}}\sum_{k=0}^{k_{f}}|x_{j}(k)|^{p}$ (10)

Let us define $W_{j}=\frac{a_{qj}}{1-a_{qq}}$. Notice that $W_{j}>0$ and
$\sum_{j=1}^{n}W_{j}\leq 1$. In this notation, the equation (10) can be
written as:

$\sum_{k=0}^{k_{f}}|x_{q}(k)|^{p}\leq\sum_{j\in
N(q)}W_{j}\sum_{k=0}^{k_{f}}|x_{j}(k)|^{p}$ (11)

The term $\sum_{k=0}^{k_{f}}|x_{j}(k)|^{p}$ is the $p$th power of the $p$-norm
of the signal $x_{j}(k)$ over the interval $[0,k_{f}]$. Since equation (11)
holds for any input, it follows that the $l_{p}$ gains satisfy the following
for any node $q$ other than the source $s$:

$(G_{p}(q,k_{f}))^{p}\leq\sum_{j\in N(q)}W_{j}(G_{p}(j,k_{f}))^{p}$ (12)

Now consider the target vertex $q^{*}$. From equation (12),
$G_{p}(q^{*},k_{f})\leq G_{p}(j,k_{f})$ for some $j\in N(q^{*})$, and further
the inequality is strict unless $G_{p}(j,k_{f})=G_{p}(q^{*},k_{f})$ for all
$j\in N(q^{*})$. Then choose a neighbor $\hat{j}$ for which $G_{p}(j,k_{f})$
is maximum. Repeating this argument for $\hat{j}$, we see that there is a
vertex other than $q^{*}$ and $\hat{j}$, say $\overline{j}$, such that
$G_{p}(\overline{j},k_{f})\leq G_{p}(\hat{j},k_{f})$. Iterating, we finally
get that $G_{p}(q^{*},k_{f})\leq G_{p}(c(i),k_{f})$ for some $i=1,2,\ldots,m$.

######  Case 2.

$p\rightarrow\infty$

The proof follows readily from Lemma 1. Specifically, consider the separating
cutset $C$ and the resulting target partition $Z$. Applying Equation (3) in
Lemma 1 and then taking the infinity norm of both sides of Equation (3), we
obtain:

$\|{\small\begin{bmatrix}X_{z}(0)\\\ \vdots\\\
X_{z}(k_{f})\end{bmatrix}}\|_{\infty}\leq\|Q_{k_{f}}\|_{\infty}\|{\small\begin{bmatrix}X_{c}(k_{f})\\\
\vdots\\\ X_{c}(0)\end{bmatrix}}\|_{\infty}$ (13)

The maximum row sum of $Q_{k_{f}}$ is $1$ and hence
$\|Q_{k_{f}}\|_{\infty}\leq 1$. As a result, we get the following inequality:

$\|{\small\begin{bmatrix}X_{z}(0)\\\ \vdots\\\
X_{z}(k_{f})\end{bmatrix}}\|_{\infty}\leq\|{\small\begin{bmatrix}X_{c}(k_{f})\\\
\vdots\\\ X_{c}(0)\end{bmatrix}}\|_{\infty}$ (14)

Equation (14) holds for all inputs $u(k)$, including the input that maximizes
the $\infty$-norm when the target node is $q^{*}$. Since $q^{*}\in Z$, it
follows immediately from (14) that $G_{p}(q^{*},k_{f})$ is upper bounded by
the maximum magnitude entry in $\begin{bmatrix}X_{c}(k_{f})\\\ \vdots\\\
X_{c}(0)\end{bmatrix}$ for this input. Therefore, $G_{p}(q^{*},k_{f})\leq
G_{p}(c(i),k_{f})$ for some $i=1,2,\ldots,m$.

∎

######  Remark 1.

The spatial pattern on the $l_{p}$ gains also holds in the infinite-horizon
limit ($G_{p}(q^{*},k_{f})$ for $k_{f}\rightarrow\infty$), provided that the
gain is finite. A finite gain is achieved if $A$ is strictly substochastic
(i.e., at least one row sum is strictly less than $1$), since the matrix $A$
is irreducible (the graph $\Gamma$ is strongly connected) by assumption.

######  Remark 2.

Per Theorem 1, the spatial degradation of the $l_{p}$ norms holds for model
(1) for any input signal. It follows that the spatial result of the $l_{p}$
gains also holds for the closed-loop system, when a feedback is applied from
the target to the source node.

######  Remark 3.

The proof also trivially extends to mixed-norm gains, since the spatial
pattern holds for any input.

Next, the response of the input-output model (Equation 1) to a periodic input,
and hence the frequency response, is also verified to exhibit a spatial
degradation:

######  Theorem 2.

Consider the zero-state (i.e. $X(0)=0$) response of the network input-output
model (1), under the assumption that the matrix $A$ is strictly substochastic
(i.e. at least one row sum is strictly less than $1$). Consider a particular
target vertex $i=q^{*}$. Also, consider any separating cutset
$C=\\{c(1),c(2),...,c(m)\\}$. Then the following statements are true:

  1. (i)

For a periodic input $u(k)$, the response at the target vertex $q^{*}$,
$x_{q^{*}}(k)\leq\max\\{0,\max_{\widehat{k}=0,1,\ldots,k}x_{c(i)}(\widehat{k})\\}$
for some $i=1,\ldots,m$. Also,
$x_{q^{*}}(k)\geq\min\\{0,\min_{\widehat{k}=0,1,\ldots,k}x_{c(j)}(\widehat{k})\\}$
for some $j=1,\ldots,m$.

  2. (ii)

At each frequency $\Omega$,
$|H_{q^{*}}(e^{j\Omega})|\leq|H_{c(i)}(e^{j\Omega})|$ for some $i=1,\ldots,m$,
where $H_{b}(j\omega)$ is the frequency response when the output is taken at
vertex $b$.

###### Proof.

Since $A$ is substochastic and also irreducible ($\Gamma$ is strongly
connected) by assumption, the system is asymptotically stable in the sense of
Lyapunov. Hence, the response to a periodic input is periodic in the
asymptote, and bounded for all time.

Let us consider the solution of equation (2) at time $k$:

$X_{z}(k)={\small\begin{bmatrix}0&B&A_{z}B&\ldots
A_{z}^{k-1}B\end{bmatrix}}{\small\begin{bmatrix}X_{c}(k)\\\ \vdots\\\
X_{c}(0)\end{bmatrix}}$ (15)

Since the target vertex $q^{*}$ is in the target vertex set $X_{z}(k)$, we can
write the corresponding response as follows:

$x_{q^{*}}(k)={\small\begin{bmatrix}0&B&A_{z}B&\ldots
A_{z}^{k-1}B\end{bmatrix}}_{q^{*}}{\small\begin{bmatrix}X_{c}(k)\\\ \vdots\\\
X_{c}(0)\end{bmatrix}}$ (16)

where ${\small\begin{bmatrix}0&B&A_{z}B&\ldots
A_{z}^{k-1}B\end{bmatrix}}_{q^{*}}$ is the row of
${\small\begin{bmatrix}0&B&A_{z}B&\ldots A_{z}^{k-1}B\end{bmatrix}}$
corresponding to the target vertex.

Notice that the row sum of ${\small\begin{bmatrix}0&B&A_{z}B&\ldots
A_{z}^{k-1}B\end{bmatrix}}$ is less than or equal to unity. Therefore, it is
immediate that for a periodic input,
$x_{q^{*}}(k)\leq\max\\{0,\max_{\widehat{k}=0,1,\ldots,k}x_{c(i)}(\widehat{k})\\}$
for all $k$ and some $i=1,\ldots,m$. Similarly, the lower bound of the target
state can be characterized in terms of cutset state as
$x_{q^{*}}(k)\geq\min\\{0,\min_{\widehat{k}=0,1,\ldots,k}x_{c(i)}(\widehat{k})\\}$,
for some $i=1,\ldots,m$.

To prove part (ii), consider the response of system (1) for a sinusoidal input
$u(k)=\cos{\Omega k}$. The response at the each node has a sinusoidal steady-
state component at the driving frequency $\Omega$, plus a transient component.
The sinusoidal component of the response at the target node is given by

$x_{q^{*},SS}(k)=|H_{q^{*}}(e^{j\Omega})|\cos({\Omega
k+\angle{H_{q^{*}}(e^{j\Omega})}}).$ (17)

Similarly, the sinusoidal component of the response at each separating cutset
node is given by:

$x_{c(i),SS}(k)=|H_{c(i)}(e^{j\Omega})|\cos({\Omega
k+\angle{H_{c(i)}(e^{j\Omega})}}).$ (18)

Now consider finding the response at the target node using Equation (2), when
the driving signal $X_{c}(k)$ is set to contain the sinusoidal components of
the cutset nodes’ state $x_{c(i),SS}(k)$ for the sinusoidal input. The
response at target node computed in this way has a transient and a sinusoidal
steady-state component. Importantly, the sinusoidal component is identical to
the sinusoidal response of the original network model (1) for the input
$u(k)=\cos{\Omega k}$, i.e. the response is
$x_{q^{*},SS}(k)=|H_{q^{*}}(e^{j\Omega})|\cos({\Omega
k+\angle{H_{q^{*}}(e^{j\Omega})}})$.

Finally, using the same argument as for other periodic inputs, one finds that
$x_{q^{*},SS}(k)\leq\max_{\widehat{k}=0,1,\ldots,k}x_{c(i),SS}(\widehat{k})$
for some $i=1,\ldots,m$, for all sufficiently large $k$. However, this is only
possible if $|H_{q^{*}}(e^{j\Omega})|\leq|H_{c(i)}(e^{j\Omega})|$ for some
$i=1,\ldots,m$.

∎

######  Remark 4.

The result in Theorem 2(i) is presented in a generalized form, to account for
any periodic input. If the maximum of the responses over the cutset vertex set
is positive, then the comparison with $0$ can be excluded in the statement for
the maximum response. Similarly, if the minimum of the responses over the
cutset vertices is negative, the comparison with 0 can be excluded in the
statement for the minimum.

######  Remark 5.

The zero ($0$) terms in the comparisons in Theorem 2(i) can be excluded, if
the state matrix $A$ is stochastic and further the asymptotic response is
considered. That is, the upper and lower bound of the target is dependent only
on the response of the cutset vertices in this case; this can be verified by
using the fact that the row sums of ${\small\begin{bmatrix}0&B&A_{z}B&\ldots
A_{z}^{k-1}B\end{bmatrix}}$ approach unity asymptotically.

The frequency-response analysis carries through to the case that $A$ is
stochastic, i.e. its row sums are unity, with only a couple of slight
limitations. This formalized in the following theorem:

######  Theorem 3.

Consider the network model (1) in the case that $A$ is an ergodic stochastic
matrix. Also, consider a target node $q^{*}$, and any separating cutset
$C=\\{c(1),c(2),...,c(m)\\}$. Then, for each frequency $\Omega>0$,
$|H_{q^{*}}(e^{j\Omega})|\leq|H_{c(i)}(e^{j\Omega})|$ for some $\Omega>0$. .

###### Proof.

Since $A$ is assumed ergodic and stochastic, it has a single eigenvalue at
unity with remaining eigenvalues strictly inside the unit circle. It follows
that the response of the system (1) to a sinusoidal input at a non-zero
frequency $\Omega$ produces a bounded sinusoidal response at the same
frequency, i.e. the frequency response is finite and well-defined at $\Omega$.
With this observation, the remainder of the proof is identical to that of
Theorem 2.

∎

Theorems 2 and 3 have characterized the spatial response properties of the
system (1) at a particular frequency. The energy in the response over a
frequency band can also be shown to exhibit a spatial falloff, by relying on
Parseval’s theorem and the spatial analysis of the two-norm response. This
result is formalized in the following theorem:

######  Theorem 4.

Consider the system (1), under the assumption that $A$ is substochastic. Also,
consider a separating cutset $C=c(1),c(2),...,c(m)$ (where the target node is
$q^{*}$). Then,
$\int_{\Omega_{1}}^{\Omega_{2}}|H_{q^{*}}(e^{j\Omega})|^{2}d\Omega\leq\int_{\Omega_{1}}^{\Omega_{2}}|H_{c(i)}(e^{j\Omega})|^{2}d\Omega$
for some $i=1,2,\ldots,m$, and any frequency band
$[\Omega_{1}\quad\Omega_{2}]$.

###### Proof.

The result is proved based on the norm inequality (11), for the special case
where $p=2$. Equation (11) in this case is give by:

$\sum_{k=0}^{k_{f}}|x_{q}(k)|^{2}\leq\sum_{j\in
N(q)}W_{j}\sum_{k=0}^{k_{f}}|x_{j}(k)|^{2}$ (19)

Now, consider the target vertex $q^{*}$ and the vertex cutset
$C=c(1),c(2),\ldots,c(m)$. Using similar arguments to those used in the proof
of theorem 1, we can obtain the following inequality:

$\sum_{k=0}^{k_{f}}|x_{q^{*}}(k)|^{2}\leq\sum_{k=0}^{k_{f}}|x_{c(i)}(k)|^{2},$
(20)

for some $i=1,\ldots,m$. The inequality holds for any input $u(k)$, and any
$k_{f}$. Considering $k_{f}\to\infty$ and then applying Parseval’s theorem
(which holds since the system is asymptotically stable), Equation (20) can be
written in the frequency domain as follows:

$\frac{1}{2\pi}\int_{-\pi}^{\pi}|X_{q}(e^{j\Omega})|^{2}d\Omega\leq\frac{1}{2\pi}\int_{-\pi}^{\pi}|X_{c(i)}(e^{j\Omega})|^{2}d\Omega$
(21)

where we use the notation $X_{i}(e^{j\Omega})$ as the Fourier transform of the
signal $x_{i}(k)$. Using the relationship
$|X_{i}(e^{j\omega})|=|H_{i}(e^{j\Omega})||U(e^{j\Omega})|$. Equation (21) can
be rewritten as:

$\int_{-\pi}^{\pi}|H_{q}(e^{j\Omega})|^{2}|U(e^{j\Omega})|^{2}d\Omega\leq\int_{-\pi}^{\pi}|H_{c(i)}(e^{j\Omega})|^{2}|U(e^{j\Omega})|^{2}d\Omega$
(22)

As Equation (22) holds for any input, let us choose the input to be
$u(k)=\frac{1}{\pi k}\sin(Wk)e^{jW_{0}k}$. The Fourier transform of this
input, $U_{1}(e^{j\Omega})$, is unity for $|\Omega-W_{0}|<W$ and zero
otherwise. Choosing $W_{0}-W=\Omega_{1}$ and $W_{0}+W=\Omega_{2}$ and
substituting into (22), we find:

$\int_{\Omega_{1}}^{\Omega_{2}}|H_{q^{*}}(e^{j\Omega})|^{2}d\Omega\leq\int_{\Omega_{1}}^{\Omega_{2}}|H_{c(i)}(e^{j\Omega})|^{2}d\Omega$
(23)

where $c(i)$ is a vertex in vertex cutset $C$.

∎

######  Remark 6.

Since Equation (22) is valid for any input, the spatial degradation pattern
also holds for any frequency-domain signal which is filtered by the network
dynamics. Specifically,
$\int_{\Omega_{1}}^{\Omega_{2}}|H_{q}(e^{j\Omega})|^{2}|F(e^{j\Omega})|^{2}d\Omega\leq\int_{\Omega_{1}}^{\Omega_{2}}|H_{c(i)}(e^{j\Omega})|^{2}|F(e^{j\Omega})|^{2}d\Omega$
for the target node $q^{*}$ and some $c(i)$ in vertex cutset $C$, where
$F(e^{j\Omega})$ is an arbitrary function which is filtered by the network.

Now, we show that the Markov parameters also exhibit a spatial degradation:

######  Theorem 5.

Consider the Markov parameters for the system (1), in the cases that the
output is taken at the target vertex $i=q^{*}$ and at vertices on a separating
cutset $C=\\{c(1),c(2),...,c(m)\\}$. The Markov parameters satisfy the
following inequality: for all $k$,
$M_{q^{*}}(k)\leq\max_{\widehat{k}=0,1,\ldots,k}M_{c(i)}(\widehat{k})$ for
some $i=1,2,\ldots,m$.

###### Proof.

The proof of theorem 3 can be derived from the proof of theorem 2. Notice that
equation (16) holds for any input. So, for any input it follows that
$x_{q^{*}}(k)\leq\max_{k}x_{c(i)}(k)$ for all $k$ and some $i=1,2,\ldots,m$ in
vertex cutset $C$. If the system is initially relaxed and driven by an impulse
input, the response of target node $q^{*}$ can be written as
$x_{q^{*}}(k)=e_{q^{*}}A^{k-1}e_{s}=M_{q^{*}}(k-1)$. Similarly, the response
in a vertex cutset node $c(i)$ can be expressed as
$x_{q^{*}}(k)=e_{c(i)}A^{k-1}e_{s}=M_{c(i)}(k-1)$. Form this, it is immediate
that $M_{q^{*}}(k)\leq\max_{k}M_{c(i)}(k)$ for some $i=1,2,\ldots,m$.

∎

We have focused in this short article on the single-input single-output system
(1), with the goal of understanding pairwise transfer characteristics in
network synchronization or diffusion processes. However, it turns out that the
spatial pattern of responses in the synchronization model holds even when
inputs are applied at multiple network nodes. To formalize this notion, we
consider a modified model where inputs are applied at multiple nodes:

$\displaystyle X(k+1)=AX(k)+BU(k)$ (24) $\displaystyle y(k)=e_{i}^{T}X(k)$

where $X(k)$, $A$, and $e_{i}$ are defined as before in Equation (1),
$U(k)=\begin{bmatrix}u_{1}(k)\\\ \vdots\\\ u_{\hat{m}}(k)\end{bmatrix}$ is a
vector of inputs, and $B$ is a $n\times\hat{m}$ matrix whose columns are $0--
1$ indicator vectors for a set of input nodes in the network (input vertices
in the network’s graph).

As an illustration, we characterize the $l_{p}$ gain for the multi-input case.
This requires a redefinition of the gain concept:

######  Definition 1.

The $l_{p}$ gain $G_{p}(i,k_{f})$ for the system (24) over the time horizon
$k_{f}$, with output taken at node $i$, is defined as the maximum of the
quantity $\left[\sum_{k=0}^{k_{f}}|y(k)|^{p}\right]^{\frac{1}{p}}$ over all
inputs $U(0),\ldots,U(k_{f}-1)$ such that
$\left[\sum_{k=0}^{k_{f}}\sum_{j=0}^{\hat{m}}|u_{j}(k)|^{p}\right]^{\frac{1}{p}}=1$.

The $l_{p}$ gains of the system in equation (24) also show a spatial falloff
defined by cutsets in the network graph, as formalized in the following
theorem:

######  Theorem 6.

Consider the $l_{p}$ gain $G_{p}(q^{*},k_{f})$ of the system (24) for a target
node $q^{*}$ and a time horizon $k_{f}$. Also, assume that
$C=\\{c(1),c(2),...,c(m)\\}$ is a separating cutset of the network, in the
sense that all directed paths in $\Gamma$ from each input vertex to the target
vertex pass through at least one of these vertices. Then, the $l_{p}$ gain of
a cutset node majorizes the $l_{p}$ gain of the target, i.e.
$G_{p}(q^{*},k_{f})\leq G_{p}(c(i),k_{f})$ for some $i=1,2,\ldots,m$.

###### Proof.

The proof of Theorem 4 is nearly identical to that of Theorem 1, hence we omit
the details.

∎

##  IV Simulations

The spatial fall-off in the $l_{p}$ gains (Theorem 1) is illustrated in an
example. A network with $9$ nodes with the following substochastic state
matrix is considered:

${\small A=\begin{bmatrix}0&0&0&0&0&0&0.2&0.3&0.25\\\
0&0&0&0&0&0&0.35&0.25&0.2\\\ 0&0&0&0&0&0&0.15&0.25&0.45\\\
0.2&0.4&0.35&0&0&0&0&0&0\\\ 0.2&0.15&0.25&0&0&0&0&0&0\\\
0.2&0.45&0.15&0&0&0&0&0&0\\\ 0&0&0&0.25&0.25&0.15&0&0&0\\\
0&0&0&0.3&0.35&0.25&0&0&0\\\ 0&0&0&0.25&0.35&0.1&0&0&0\\\ \end{bmatrix}}$

The network’s digraph, which has 27 directed edges, is shown in Figure 2. The
source vertex (Vertex 1) is also indicated.

Figure 2: The digraph of the example network.

To illustrate Theorem 1, the $l_{1}$ gain is computed when each node in the
network is selected as the target. These gains are plotted on top of the
network’s graph in Figure 3. The plot shows the decrescence in the gains along
cutsets in the network graph. As one example, the $l_{1}$ gain when vertex $9$
is the target is upper-bounded by the maximum among the $l_{1}$ gains when
vertices $4$, $5$, and $6$ are chosen as targets. Since these vertices form a
separating cutset, the example matches with the formal result.

Figure 3: The $l_{1}$ gains for different target nodes.

The decrescence of $l_{p}$ gains along cutsets also immediately implies that
the maximum gain among target vertices at a certain distance from the source
is a non-increasing function of the distance. This is true because the set of
vertices at each distance forms a separating cutset between the source vertex
and more distant vertices. Thus, it is insightful to calculate maximum $l_{p}$
gains for target vertices at different distances. We have done so for the $1$,
$2$, and $\infty$ gains in Table I. Each $l_{p}$ gain is highest at the source
(distance $0$), and decreases with distance from source.

Distance | Max $l_{1}$ Gain | Max $l_{2}$ Gain | Max $l_{\infty}$ Gain
---|---|---|---
0 | 1.2122 | 1.1676 | 1.2122
1 | 0.4098 | 0.3654 | 0.4098
2 | 0.3328 | 0.2994 | 0.3328
3 | 0.2348 | 0.2113 | 0.2348
TABLE I: Decrescence of $l_{p}$ gains with distance from the source.

##  V Applications of the Spatial Pattern Analysis

Two brief vignettes are presented, to illustrate potential applications of the
spatial analyses developed in the short paper. The first is concerned with
signal recovery from remote measurements in a network, while the second is
concerned with definitions of (external) propagation stability in networks.

###  V-A Signal Recovery and Signal-to-Noise Ratios in Diffusive Networks

The problem of estimating or detecting an external input to a network
synchronization or diffusion process from noisy remote measurements arises in
several settings[14]. For instance, it may be necessary to locate and
characterize a pollution source in a complex environment, based on localized
measurements of the diffusing pollutant. Similarly, it may be of interest to
identify a stubborn or manipulative agent in a network opinion dynamics, based
on imperfect measurements of certain agents’ opinions. Standard machinery for
detection/estimation (e.g., hypothesis testing, Wiener filtering) can be
brought to bear for these problems. In turn, the statistical performance of
the recovery technique can be characterized, and used to develop algorithms
for sensor placement. However, these analyses – particularly ones for
performance analysis and sensor placement – can become challenging to develop
and/or computationally infeasible, particularly when the network is high
dimensional or incompletely known. For these reasons, graph-theoretic insights
into network signal estimation/detection performance are valuable.

Broadly, the performance of estimators and detectors is often closely tied
with the signal-to-noise ratio (SNR) in the measured signal, i.e. the ratio
between the energy/power of the signal of interest and that of the noise. The
characterizations of input-output responses developed in this article imply
that signal-to-noise ratios also degrade spatially in a network away from a
signal source, and hence suggests that detection/estimation performance
degrades monotonically away from the source.

To present the concept more explicitly, let us consider the network dynamics
(1) when the input $u(k)$ at the source node is an unknown signal, and
measurements at a remote target node are subject to an additive zero-mean
white noise with variance $\sigma^{2}$. That is, we assume that measurements
taken at a node $i$ are of the form $\overline{y}(k)=x_{i}(k)+v(k)$, where
$v(k)$ is a zero-mean white noise signal.

In the case where the the input $u(k)$ is a finite-energy or transient signal,
it is natural to define the SNR at the target node as the ratio between the
energy (or other norm) of the response signal $x_{i}(k)$ and the noise
intensity, e.g. $SNR=\frac{||x_{i}(k)||_{p}}{\sigma}$. From Theorem 1, it is
immediate that the SNR for any input signal (and the best-case SNR among all
inputs with a certain energy/norm) exhibits the spatial degradation pattern
defined by graph cutsets. That is, as compared to the SNR at a target node,
the SNR will be higher for at least one node on a cutset separating the source
and target. Thus, estimation or detection performance should be improved on
the separating cutset, i.e. at nodes closer to the source.

In the case where the input $u(k)$ is a persistent signal, the SNR is
typically defined in terms of the power of the response signal $x_{i}(k)$ and
the noise. A typical case may be that the input and hence the the response
signal $x_{i}(k)$ is a stochastic band-limited signal. In this case, the total
power in $x_{i}(k)$ can be determined as the integral of the power spectrum
over the band. Thus, the SNR can be found as:
$SNR=\frac{\int_{\Omega}|H_{i}(e^{j\Omega})|^{2}S_{UU}(\Omega)|d\Omega}{\sigma^{2}}$,
where $S_{UU}(\Omega)$ is the power spectrum of the input signal $u(k)$. From
Theorem 3 and the following remark (Remark 2), we immediately see that the SNR
in this case also exhibits a spatial degradation along cutsets away from the
source. Thus, estimation or detection is again seen to be improved at cutsets
near the source node.

We stress that the spatial degradation of signal strengths and SNRs holds not
only for the $2$-norm but for all other $p$-norms, which are relevant for a
number of signal reconstruction and detection problems [15].

###  V-B Network Propagation Stability Analysis

The stability analysis of network synchronization models has been primarily
focused on internal (Lyapunov) stability notions, which capture emergence of a
synchronous state [2, 1, 16]. However, in many application areas, the
responses of network processes to exogenous disturbances is of substantial
interest [17, 18, 19]. With this motivation, a number of studies have
considered external stability of network synchronization models [20, 21].
These studies generally define external stability as a bounded-input bounded-
output (BIBO) stability notion, with the assumption that exogenous
(disturbance) inputs are present at any node or a set of leader nodes, and the
output is the full network’s state or projections thereof.

For many network processes, a primary concern is whether a disturbance at one
node can amplify and propagate across the network, or whether alternately it
remains localized. For instance, power-system operators recognize that
oscillatory disruptions originating from generator controls can couple with
the system’s natural modes to cause system-wide oscillations. The propagation
or amplification of a disturbance across a network is distinct from both the
internal stability and external stability concepts. This distinction motivates
alternative definitions of stability for networks, concerned with disturbance
propagation.

In fact, propagation of disturbances in cascaded systems and bidirectional
line networks has been extensively studied, under the heading of string
stability [22, 23, 6]. Several subtly different definitions of string
stability have been proposed. Broadly, the definitions require that a
disturbance originating at one location in the string can only be amplified by
a specified finite gain for any output location, regardless of the length of
the string. A stronger notion, referred to as strict string stability,
requires that the disturbance response is diminished (scaled down) at each
subsequent node in the string [24, 25]. The string stability notion has been
extended for other network topologies in a couple of ways. First, there is a
considerable literature on mesh stability, which generalizes the string-
stability concept to general directional networks [26, 27]. Also, an important
recent study by Studli and co-workers [3] has considered generalization of
string stability concepts to more general networks. This study provides a
parallel definition to string stability for general networks, as well as a
parallel to strict string stability for tree-like networks, although without
characterizations for particular network models. Other input-to-state
stability concepts with a similar flavor, which generalize mesh stability,
have also recently been developed for networks with general graph topologies
and nonlinear nodal dynamics [28].

We argue that our spatial input-output analysis immediately suggests a
definition for strict stability with respect to disturbance propagation, for
general network processes. Specifically, a network dynamics can be viewed as
strictly propagation stable, if responses degrade along cutsets away from the
source. This notion is formalized in the following definition:

######  Definition 2.

Consider a dynamical network process where an input ${\bf u}(k)$ is applied at
a source node $s$, and the state responses ${\bf y}_{q^{*}}(k)$ are considered
at target nodes $q^{*}$ in the network. The network is ${\cal L}_{p,t}$
strictly propagation stable if, for every source node $s$, target node
$q^{*}$, separating cutset $C=\\{c(1),c(2),...,c(m)\\}$, and finite-$p$-norm
input ($||u(k)||_{p}<\infty$), it holds that $||{\bf
y}_{q^{*}}(k)||_{t}\leq||{\bf y}_{c(i)}(k)||_{t}$ for some $i=1,\ldots,m$.

The definition for propagation stability can easily be rephrased in terms of
paths in the network’s graph. In particular, a network is ${\cal L}_{p.t}$
strictly propagation stable if and only if, for every source vertex $s$,
target vertex $q^{*}$, and finite-$p$-norm input, there exists a path from $s$
to $q^{*}$ in the network graph such that the $t$-norm of the response at the
corresponding network nodes is decreasing along the path. From this phrasing,
it is evident that the concept reduces to the standard notions for strict
stability for strings and tree graphs. The statement also clarifies that the
definition is flexible to the complex propagations that may occur in networks,
in that decrescence is only required along one path between a source and
target node (not all paths).

From the proof of Theorem 1, the scalar network synchronization model (1) is
immediately seen to be ${\cal L}_{p,t}$ strictly propagation stable. Thus, we
see that disturbances in this model are localized/diminishing in the sense of
the strict propagation stability definition, regardless of the specific
network topology. For more complex synchronization models involving
multivariate nodal dynamics, one anticipates that the node or subsystem model
and the network graph will together determine whether or not strict
propagation stability holds.

##  VI Conclusions

This short note has focused on the input-to-output analysis of a canonical
discrete time network synchronization model. Our key finding is the a family
of input-to-output metrics (e.g. $l_{p}$ gains, frequency responses,
frequency-band energy, Markov parameters) show a spatial degradation across
cutsets away from the input location. Two applications of the analysis have
been briefly developed: 1) characterization of signal-to-noise ratios in
network processes, and 2) definition of propagation stability notions for
dynamical networks. Of interest, the spatial analyses carry through to a
number of time-varying and nonlinear processes. We expect to address these
cases in future work.

## References

  * [1] L. M. Pecora and T. L. Carroll, “Master stability functions for synchronized coupled systems,” _Physical review letters_ , vol. 80, no. 10, p. 2109, 1998\.
  * [2] W. Ren, R. W. Beard, and E. M. Atkins, “Information consensus in multivehicle cooperative control,” _IEEE Control systems magazine_ , vol. 27, no. 2, pp. 71–82, 2007.
  * [3] S. Stüdli, M. M. Seron, and R. H. Middleton, “From vehicular platoons to general networked systems: String stability and related concepts,” _Annual Reviews in Control_ , vol. 44, pp. 157–172, 2017.
  * [4] B. Besselink and S. Knorn, “Scalable input-to-state stability for performance analysis of large-scale networks,” _IEEE Control Systems Letters_ , vol. 2, no. 3, pp. 507–512, 2018.
  * [5] M. Mirabilio, A. Iovine, E. De Santis, M. D. Di Benedetto, and G. Pola, “Scalable mesh stability of nonlinear interconnected systems,” _IEEE Control Systems Letters_ , vol. 6, pp. 968–973, 2021.
  * [6] D. Swaroop and J. K. Hedrick, “String stability of interconnected systems,” _IEEE transactions on automatic control_ , vol. 41, no. 3, pp. 349–357, 1996\.
  * [7] M. Pirani, J. W. Simpson-Porco, and B. Fidan, “System-theoretic performance metrics for low-inertia stability of power networks,” in _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_. IEEE, 2017, pp. 5106–5111.
  * [8] K. Koorehdavoudi, M. Hatami, S. Roy, V. Venkatasubramanian, P. Panciatici, F. Xavier, and J. A. Torres, “Input-output characteristics of the power transmission network’s swing dynamics,” in _2016 IEEE 55th Conference on Decision and Control (CDC)_. IEEE, 2016, pp. 1846–1852.
  * [9] R. N. Mahia and D. M. Fulwani, “On some input–output dynamic properties of complex networks,” _IEEE Transactions on Circuits and Systems II: Express Briefs_ , vol. 65, no. 2, pp. 216–220, 2017.
  * [10] J. Abad Torres and S. Roy, “Graph-theoretic characterisations of zeros for the input–output dynamics of complex network processes,” _International Journal of Control_ , vol. 87, no. 5, pp. 940–950, 2014.
  * [11] A. Vosughi, C. Johnson, M. Xue, S. Roy, and S. Warnick, “Target control and source estimation metrics for dynamical networks,” _Automatica_ , vol. 100, pp. 412–416, 2019.
  * [12] F. Pasqualetti, F. Dörfler, and F. Bullo, “Cyber-physical attacks in power networks: Models, fundamental limitations and monitor design,” in _2011 50th IEEE Conference on Decision and Control and European Control Conference_. IEEE, 2011, pp. 2195–2201.
  * [13] L. Xiao and S. Boyd, “Fast linear iterations for distributed averaging,” _Systems & Control Letters_, vol. 53, no. 1, pp. 65–78, 2004.
  * [14] L. Ye, S. Roy, and S. Sundaram, “Optimal sensor placement for kalman filtering in stochastically forced consensus networks,” in _2018 IEEE Conference on Decision and Control (CDC)_. IEEE, 2018, pp. 6686–6691.
  * [15] V. R. S. Banjade, C. Tellambura, and H. Jiang, “Performance of $p$-norm detector in awgn, fading, and diversity reception,” _IEEE Transactions on Vehicular Technology_ , vol. 63, no. 7, pp. 3209–3222, 2014.
  * [16] Z. Li and G. Chen, “Global synchronization and asymptotic stability of complex dynamical networks,” _IEEE Transactions on Circuits and Systems II: Express Briefs_ , vol. 53, no. 1, pp. 28–33, 2006.
  * [17] J. Lu and D. W. Ho, “Stabilization of complex dynamical networks with noise disturbance under performance constraint,” _Nonlinear Analysis: Real World Applications_ , vol. 12, no. 4, pp. 1974–1984, 2011.
  * [18] B. Liu, C. Xu, and D. Liu, “Input-to-state-stability-type comparison principles and input-to-state stability for discrete-time dynamical networks with time delays,” _International Journal of Robust and Nonlinear Control_ , vol. 23, no. 4, pp. 450–472, 2013.
  * [19] M. Siami and N. Motee, “Fundamental limits and tradeoffs on disturbance propagation in linear dynamical networks,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 12, pp. 4055–4062, 2016.
  * [20] Z. Li, Z. Duan, and G. Chen, “On ${H}_{\infty}$ and ${H}_{2}$ performance regions of multi-agent systems,” _Automatica_ , vol. 47, no. 4, pp. 797–803, 2011.
  * [21] J. Wang, Y. Tan, and I. Mareels, “Robustness analysis of leader-follower consensus,” _Journal of systems science and complexity_ , vol. 22, no. 2, pp. 186–206, 2009.
  * [22] L. Peppard, “String stability of relative-motion pid vehicle control systems,” _IEEE Transactions on Automatic Control_ , vol. 19, no. 5, pp. 579–581, 1974.
  * [23] B. Besselink and K. H. Johansson, “String stability and a delay-based spacing policy for vehicle platoons subject to disturbances,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 9, pp. 4376–4391, 2017.
  * [24] J. Ploeg, N. Van De Wouw, and H. Nijmeijer, “Lp string stability of cascaded systems: Application to vehicle platooning,” _IEEE Transactions on Control Systems Technology_ , vol. 22, no. 2, pp. 786–793, 2013.
  * [25] J. A. Rogge and D. Aeyels, “Vehicle platoons through ring coupling,” _IEEE Transactions on Automatic Control_ , vol. 53, no. 6, pp. 1370–1377, 2008.
  * [26] A. Pant, P. Seiler, and K. Hedrick, “Mesh stability of look-ahead interconnected systems,” _IEEE Transactions on Automatic Control_ , vol. 47, no. 2, pp. 403–407, 2002.
  * [27] M. Mirabilio, A. Iovine, E. De Santis, M. D. D. Benedetto, and G. Pola, “Scalable mesh stability of nonlinear interconnected systems,” _IEEE Control Systems Letters_ , vol. 6, pp. 968–973, 2022.
  * [28] B. Besselink and S. Knorn, “Scalable input-to-state stability for performance analysis of large-scale networks,” _IEEE Control Systems Letters_ , vol. 2, no. 3, pp. 507–512, 2018.


# Learning a Gaussian Mixture Model from Imperfect Training Data for Robust
Channel Estimation

Benedikt Fesl, Nurettin Turan, Michael Joham,
and Wolfgang Utschick This work was partly funded by Huawei Sweden
Technologies AB, Lund. The authors are with Professur für Methoden der
Signalverarbeitung, Technische Universität München, 80333 München, Germany
(e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>utschick@tum.de).

###### Abstract

In this letter, we propose a Gaussian mixture model (GMM)-based channel
estimator which is learned on imperfect training data, i.e., the training data
are solely comprised of noisy and sparsely allocated pilot observations. In a
practical application, recent pilot observations at the base station (BS) can
be utilized for training. This is in sharp contrast to state-of-the-art
machine learning (ML) techniques where a training dataset consisting of
perfect channel state information (CSI) samples is a prerequisite, which is
generally unaffordable. In particular, we propose an adapted training
procedure for fitting the GMM which is a generative model that represents the
distribution of all potential channels associated with a specific BS cell. To
this end, the necessary modifications of the underlying expectation-
maximization (EM) algorithm are derived. Numerical results show that the
proposed estimator performs close to the case where perfect CSI is available
for the training and exhibits a higher robustness against imperfections in the
training data as compared to state-of-the-art ML techniques.

###### Index Terms:

Robust channel estimation, imperfect data,
generative model, Gaussian mixture, OFDM system.

##  I Introduction

©This work has been submitted to the IEEE for possible publication. Copyright
may be transferred without notice, after which this version may no longer be
accessible.

Channel estimation plays a crucial role in enhancing wireless communications
systems. Recently, ML approaches were successfully leveraged for channel
estimation. The goal is to exploit a priori information about all possible
channels of mobile terminals (MTs) associated with a specific BS cell and its
radio propagation environment to improve the channel estimation quality. This
is generally intractable to model analytically but is represented in terms of
a training dataset that is available at the BS.

Thereby, different learning techniques can be distinguished. In end-to-end
learning, the channel estimation is not performed explicitly, but a network is
trained to directly perform signal detection [1]. A different approach is to
learn a nonlinear regression mapping from the pilot observation in the input
to a channel estimate at the output of a network by utilizing ground-truth CSI
samples [2, 3]. In contrast to that, it was recently proposed to train a
generative model that represents the channel distribution of the whole BS
cell, which is afterwards leveraged for channel estimation [4, 5].

A common prerequisite of ML-based approaches is the availability of a
representative training dataset consisting of perfect CSI samples. However,
the construction of such a dataset is a challenging task in practice. One
possibility is to perform costly measurement campaigns for each BS, which is
generally unaffordable. A different attempt is to use channel simulators,
e.g., [6]. The inherent problem of simulators is the mismatch between the real
and the artificially generated data, leading to performance losses. Another
problem is that the BS environment may change over time which is difficult to
track. This leads to the idea of utilizing pilot observations, which are
received in great numbers at the BS during regular operation and capture the
environmental information, as training data. Although the data aggregation is
then cheap and the dataset can be continuously updated to account for varying
conditions, the dataset contains imperfections, e.g., noise and sparse pilot
allocations. However, these imperfections can be mitigated by specific
training adaptations since they follow a known model.

Contributions: In this letter, we propose an adaptation of the EM algorithm to
fit a GMM with noisy, and possibly sparsely allocated, pilot observations. We
derive new update steps for the GMM parameters that take the system model,
i.e., the model of the imperfections, into account. Additionally, we show that
imposing structural features to the covariances of the GMM acts as a
regularization which enhances the channel estimation performance, especially
for sparse pilot allocations. It is discussed that the presented method is
particularly robust against imperfections in the training data, in contrast to
commonly used regression-based ML techniques. Finally, the discussed
properties are verified with simulations for both a spatial and a doubly-
selective fading model.

Notation: The identity matrix of dimension $N\times N$ is denoted by
$\operatorname{\mathbf{I}}_{N}$. The column-wise vectorization and the trace
of a matrix is denoted by $\operatorname{vec}(\cdot)$ and
$\operatorname{tr}(\cdot)$, respectively. We denote a positive semi-definite
(PSD) matrix ${\bm{C}}$ which fulfills
${\bm{x}}^{\operatorname{H}}{\bm{C}}{\bm{x}}\geq 0$ for all
${\bm{x}}\neq{\bm{0}}$ as ${\bm{C}}\succeq{\bm{0}}$.

##  II System and Channel Models

We consider the general system model of pilot observations

${\bm{y}}={\bm{A}}{\bm{h}}+{\bm{n}},$ (1)

where
${\bm{A}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\in\mathbb{C}^{M\times
N}}$ is the known observation matrix, ${\bm{h}}\in\mathbb{C}^{N}$ is the
wireless channel with an unknown probability density function (PDF)
$f_{{\bm{h}}}$, and additive white Gaussian noise (AWGN)
${\bm{n}}\sim\mathcal{N}_{\mathbb{C}}({\bm{0}},{\bm{C}}_{{\bm{n}}}=\sigma^{2}\operatorname{\mathbf{I}}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}M}})$.
In this letter, we consider the following two instances of (1).

###  II-A Spatial System Model

Consider a single-input multiple-output (SIMO) system where a BS equipped with
$N$ antennas serves single-antenna MTs in a single-user uplink transmission.
After correlating the pilot sequence for a single snapshot we obtain
${\bm{y}}={\bm{h}}+{\bm{n}}$, following the model in (1) with $M=N$ and
${\bm{A}}=\operatorname{\mathbf{I}}_{N}$.

We work with a spatial channel model [3] where channels are modeled
conditionally Gaussian:
${\bm{h}}|{\bm{\delta}}\sim\mathcal{N}_{\mathbb{C}}({\bm{0}},{\bm{C}}_{{\bm{\delta}}})$.
The random vector ${\bm{\delta}}$ collects the angles of arrival and path
gains of the main propagation clusters between a MT and the BS. A different
spatial channel covariance matrix ${\bm{C}}_{{\bm{\delta}}}$ is computed for
each sample by means of the steering vector of a uniform linear array (ULA)
antenna array at the BS, cf. [3].

###  II-B Doubly-Selective Fading System Model

In this case, we consider a single-input single-output (SISO) transmission in
the spatial domain over a doubly-selective fading channel
${\bm{h}}=\operatorname{vec}({\bm{H}})$, where
${\bm{H}}\in\mathbb{C}^{N_{c}\times N_{t}}$ represents the time-frequency
response of the channel for $N_{c}$ carriers and $N_{t}$ time slots. This is a
typical setup in orthogonal frequency-division multiplexing (OFDM) systems.
When only $N_{p}$ positions of the $N_{t}\times N_{c}$ time-frequency response
are occupied by pilot symbols, there is a selection matrix
${\bm{A}}\in\\{0,1\\}^{N_{p}\times N_{c}N_{t}}$ which represents the pilot
positions. This leads to observations as described in (1) where $N=N_{c}N_{t}$
and $M=N_{p}$. In this work, we consider a diamond-shaped pilot allocation
scheme which is known to be mean square error (MSE)-optimal [7].

For the construction of a scenario-specific channel dataset, we use the
QuaDRiGa channel simulator [6]. We consider an urban macrocell (UMa) scenario
following the 3GPP 38.901 specification, where the BS is placed at a height of
25m and covers a sector of 120°. Each MT is either placed indoors (80%) or
outdoors (20%) and moves with a certain velocity $v$ in a random direction,
which is captured by a drifting model.

##  III Learning a GMM from Imperfect Data

###  III-A Scenario-Specific Imperfect Training Dataset

Having access to training data that represent the environment of a BS cell,
i.e., the typically complex and intractable channel PDF $f_{{\bm{h}}}$, in
combination with ML techniques, has been shown to substantially improve the
channel estimation performance. The common prerequisite of these data-aided
techniques is the availability of a training dataset $\mathcal{H}$ consisting
of perfect CSI samples, i.e., $L$ channel realizations
$\mathcal{H}=\\{{\bm{h}}_{\ell}\\}_{\ell=1}^{L}$. However, this assumption is
impractical if the training data are collected at a BS during regular
operation. That is, the training data stem from certain pilot positions at a
finite signal-to-noise ratio (SNR), i.e.,
$\mathcal{Y}=\\{{\bm{y}}_{\ell}\\}_{\ell=1}^{L}$, cf. (1). Note that these
data can be pre-processed in different ways, e.g., by denoising or
interpolating the pilot positions. Although these training data are corrupted
by various imperfections, the statistical information about the system model
can be utilized to mitigate these effects, as shown in the following.

###  III-B GMM-based Channel Estimator with Adapted Training

In [4], a channel estimator was introduced that is based on a GMM whose
application consists of two phases. First, in an offline training phase a
$K$-components GMM of the form

$f_{{\bm{h}}}^{(K)}({\bm{h}})=\sum_{k=1}^{K}\pi_{k}\mathcal{N}_{\mathbb{C}}({\bm{h}};{\bm{\mu}}_{k},{\bm{C}}_{k}),$
(2)

where $\pi_{k}$, ${\bm{\mu}}_{k}$, and ${\bm{C}}_{k}$ are the mixing
coefficients, means, and covariances of the $k$th GMM component, respectively,
is fitted by maximum likelihood optimization with the well-known EM algorithm,
cf. [8, Sec. 9.2], via training samples from $\mathcal{H}$ in order to
approximate the underlying distribution $f_{{\bm{h}}}$ of the channels in the
whole BS cell. Due to the Gaussianity of the noise with the known covariance
matrix ${\bm{C}}_{{\bm{n}}}$ and the known observation matrix ${\bm{A}}$, the
PDF of the observations $f_{{\bm{y}}}$ can be approximated by means of the GMM
$f_{{\bm{y}}}^{(K)}({\bm{y}})=\sum_{k=1}^{K}\pi_{k}\mathcal{N}_{\mathbb{C}}({\bm{y}};{\bm{A}}{\bm{\mu}}_{k},{\bm{C}}_{{\bm{y}},k})$,
where
${\bm{C}}_{{\bm{y}},k}={\bm{A}}{\bm{C}}_{k}{\bm{A}}^{\operatorname{H}}+{\bm{C}}_{{\bm{n}}}$,
which can be straightforwardly computed once (2) is known, cf. (1). Second,
the trained GMM is leveraged to perform channel estimation by computing a
convex combination of linear minimum mean square error (MMSE) estimates for a
given pilot signal ${\bm{y}}$, i.e.,

$\hat{{\bm{h}}}_{\text{GMM}}=\sum_{k=1}^{K}\gamma_{k}({\bm{y}})\left({\bm{C}}_{k}{\bm{A}}^{\operatorname{H}}{\bm{C}}_{{\bm{y}},k}^{-1}({\bm{y}}-{\bm{A}}{\bm{\mu}}_{k})+{\bm{\mu}}_{k}\right)$
(3)

where $\gamma_{k}({\bm{y}})$ is the responsibility of the $k$th component,
i.e., the probability that the component $k$ is responsible for the
observation ${\bm{y}}$ [8, Sec. 9.2]. In [9], the estimator was extended to
impose structural features to the channel covariances ${\bm{C}}_{k}$ in the
training of the GMM for saving complexity and memory.

However, so far, the GMM was trained on a dataset comprised of perfect
channels $\mathcal{H}$, cf. Section III-A. If we naively replace $\mathcal{H}$
with the noisy, possibly pre-processed, training data $\mathcal{Y}$ from pilot
observations (1), a severe performance loss can be expected due to the
imperfections in the training data. Therefore, in the following, we propose
adapted training procedures for the GMM that drastically alleviates the
effects of the imperfect training data by utilizing the knowledge of the
system model (1) and possibly existing structural features of the covariances.

We first discuss the case of training data that is corrupted by AWGN following
the model in (1) with
${\bm{A}}=\operatorname{\mathbf{I}}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}N}}$,
cf. Section II-A. The main idea is to adapt the EM algorithm such that the PDF
$f_{{\bm{y}}}$ of the observations is approximated by the GMM
$f_{{\bm{y}}}^{(K)}({\bm{y}})$. Due to zero-mean AWGN, the updates of the
mixing coefficients $\pi_{k}$ and the means ${\bm{\mu}}_{k}$ are unchanged
with respect to the classical EM algorithm, cf. [8, Sec. 9.2]. However, we
would like to include the constraint that
${\bm{C}}_{{\bm{y}},k}={\bm{C}}_{k}+{\bm{C}}_{{\bm{n}}}$, following the model
(1). The GMM for the channel distribution (2) can then be directly obtained
since ${\bm{C}}_{{\bm{n}}}$ is known. The derivation for the update of
${\bm{C}}_{k}$ is given in the following.

######  Theorem 1.

Given noisy pilot observations $\mathcal{Y}$, the maximum likelihood solution
${\bm{C}}_{k}^{*}$ for the covariance in the EM algorithm is given by
computing the eigenvalue decomposition (EVD)

$\hat{{\bm{C}}}_{{\bm{y}},k}-{\bm{C}}_{{\bm{n}}}={\bm{V}}\operatorname{diag}({\bm{\xi}}){\bm{V}}^{\operatorname{H}}$
(4)

with
$\hat{{\bm{C}}}_{{\bm{y}},k}=\frac{1}{N_{k}}\sum_{\ell=1}^{L}\gamma_{k,\ell}({\bm{y}}_{\ell}-{\bm{\mu}}_{k})({\bm{y}}_{\ell}-{\bm{\mu}}_{k})^{\operatorname{H}}$,
$N_{k}=\sum_{\ell=1}^{L}\gamma_{k,\ell}$, and $\gamma_{k,\ell}$ being the
responsibility of component $k$ for data point ${\bm{y}}_{\ell}$. Afterwards,
an elementwise truncation of the negative eigenvalues via
${\bm{\xi}}^{\text{PSD}}=\max({\bm{0}},{\bm{\xi}})$ is performed such that

${\bm{C}}_{k}^{*}={\bm{V}}\operatorname{diag}({\bm{\xi}}^{\text{PSD}}){\bm{V}}^{\operatorname{H}}.$
(5)

###### Proof.

The maximization of the expected complete log-likelihood, cf. [8, Sec. 9.3],
of the pilot observations with respect to the $k$th channel covariance
${\bm{C}}_{k}$ is given by

$\displaystyle{\bm{C}}_{k}^{*}$
$\displaystyle=\operatorname*{arg\,max}_{{\bm{C}}_{k}\succeq{\bm{0}}}\sum_{j=1}^{K}\sum_{\ell=1}^{L}\gamma_{j,\ell}\log\mathcal{N}_{\mathbb{C}}({\bm{y}}_{\ell};{\bm{\mu}}_{j},{\bm{C}}_{j}+{\bm{C}}_{{\bm{n}}})$
(6)
$\displaystyle=\operatorname*{arg\,max}_{{\bm{C}}_{k}\succeq{\bm{0}}}\sum_{\ell=1}^{L}\gamma_{k,\ell}\log\mathcal{N}_{\mathbb{C}}({\bm{y}}_{\ell};{\bm{\mu}}_{k},{\bm{C}}_{k}+{\bm{C}}_{{\bm{n}}})$
(7)
$\displaystyle=\operatorname*{arg\,max}_{{\bm{C}}_{k}\succeq{\bm{0}}}-\log\det({\bm{C}}_{k}+{\bm{C}}_{{\bm{n}}})-\operatorname{tr}(\hat{{\bm{C}}}_{{\bm{y}},k}({\bm{C}}_{k}+{\bm{C}}_{{\bm{n}}})^{-1})$
(8)

where from (6) to (7) it is used that the maximization depends only on the
$k$th GMM component, and (8) follows from the definition of the Gaussian PDF
and $\hat{{\bm{C}}}_{{\bm{y}},k}$ from above. Thus, the optimization
simplifies to a maximization of the Gaussian likelihood of the $k$th GMM
component.  The optimal solution for the case of a Gaussian likelihood is
derived in [10, Appendix] and is obtained with the EVD and the truncation of
the negative eigenvalues. This can be applied to (8) through a substitution of
$\hat{{\bm{C}}}_{{\bm{y}},k}$ which yields the solution in (5). ∎

Next, we consider the case where the pilot observations stem from the OFDM
model in Section II-B, i.e., only a few channel entries are observed according
to the pilot pattern represented by ${\bm{A}}$. To adapt the training
procedure, we combine the insights from Theorem 1 together with the EM
algorithm for missing entries, cf. [11, Chapter 11]. We first define a
selection matrix
$\bar{{\bm{A}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\in\\{0,1\\}^{(N_{c}N_{t}-N_{p})\times
N_{c}N_{t}}}$ which represents the positions that are not allocated with
pilots, i.e., ${\bm{A}}\bar{{\bm{A}}}^{\operatorname{T}}={\bm{0}}$. After an
initial interpolation step, one EM iteration is performed in the usual way to
have initial estimates of the GMM parameters.  The updates of the mixing
coefficients $\pi_{k}$ and means ${\bm{\mu}}_{k}$ are again unchanged with
respect to the classical EM algorithm due to the zero-mean AWGN. The adapted
updates of the covariances ${\bm{C}}_{k}$ are derived in the following.

######  Theorem 2.

Given sparsely allocated and noisy pilot observations $\mathcal{Y}$, first, a
linear MMSE estimate of the unobserved and noisy channel entries
$\bar{{\bm{y}}}_{\ell,k}$ is computed via the current statistics of the $k$th
GMM component as

$\bar{{\bm{y}}}_{k,\ell}=\bar{{\bm{A}}}{\bm{\mu}}_{k}+\bar{{\bm{A}}}{\bm{C}}_{k}{\bm{A}}^{\operatorname{T}}({\bm{A}}{\bm{C}}_{k}{\bm{A}}^{\operatorname{T}}+{\bm{C}}_{{\bm{n}}})^{-1}({\bm{y}}_{\ell}-{\bm{A}}{\bm{\mu}}_{k}).$
(9)

Afterwards, a fully interpolated sample, given as
$\hat{{\bm{y}}}_{k,\ell}={\bm{A}}^{\operatorname{T}}{\bm{y}}_{\ell}+\bar{{\bm{A}}}^{\operatorname{T}}\bar{{\bm{y}}}_{k,\ell}$,
is used to update

$\hat{{\bm{C}}}_{{\bm{y}},k}=\frac{1}{N_{k}}\sum_{\ell=1}^{L}\gamma_{k,\ell}(\hat{{\bm{y}}}_{k,\ell}-{\bm{\mu}}_{k})(\hat{{\bm{y}}}_{k,\ell}-{\bm{\mu}}_{k})^{\operatorname{H}}+\bar{{\bm{A}}}^{\operatorname{T}}{\bm{\Sigma}}_{k}\bar{{\bm{A}}}$
(10)

where $N_{k}=\sum_{\ell=1}^{L}\gamma_{k,\ell}$ and

${\bm{\Sigma}}_{k}=\bar{{\bm{A}}}{\bm{C}}_{k}\bar{{\bm{A}}}^{\operatorname{T}}-\bar{{\bm{A}}}{\bm{C}}_{k}{\bm{A}}^{\operatorname{T}}({\bm{A}}{\bm{C}}_{k}{\bm{A}}^{\operatorname{T}}+{\bm{C}}_{{\bm{n}}})^{-1}{\bm{A}}{\bm{C}}_{k}\bar{{\bm{A}}}^{\operatorname{T}}.$
(11)

Finally, to account for the AWGN, the update of the channel covariance matrix
${\bm{C}}_{k}^{*}$ is computed via the EVD of
$\hat{{\bm{C}}}_{{\bm{y}},k}-{\bm{C}}_{{\bm{n}}}$ and the truncation of the
negative eigenvalues, cf. Theorem 1.

###### Proof.

The steps (9)–(11) are derived in the EM algorithm for missing data, cf. [11,
Chapter 11], where the additional covariance term (11) accounts for the
estimated covariance of the missing entries. Given the estimate of the
observation covariance $\hat{{\bm{C}}}_{{\bm{y}},k}$ in (10), the subsequent
projection via the EVD to get the maximum likelihood estimation of the channel
covariance ${\bm{C}}_{k}^{*}$ is a direct consequence of Theorem 1. ∎

Additionally, it is possible to enforce covariances with a specific structure
[9]. In this letter, we are focusing especially on the case of block-Toeplitz
covariances, constructed as
${\bm{C}}_{k}={\bm{Q}}^{\operatorname{H}}\operatorname{diag}({\bm{c}}_{k}){\bm{Q}}$,
where ${\bm{Q}}$ is a truncated 2D-discrete Fourier transform (DFT) matrix.
This structure fits particularly well in the OFDM case [9]. It is important to
mention that the imposed structural constraints on the covariances can be
understood as a regularization technique. Interestingly, as shown later in
Section IV, this regularization allows for performance gains, especially in
the case of sparse pilot allocations.

Algorithm 1 summarizes the proposed adapted EM algorithm for fitting the GMM
in the OFDM case from Section II-B with noisy and sparsely allocated pilot
observations with structured covariances. If we consider the spatial model
from Section II-A, the steps 8, 9, and 14 can be omitted. For unstructured
covariances, the steps 2, 20, 21, and 22 are dropped.

###  III-C Discussion about Robustness

The first stage of training the GMM is equivalent to learn a generative model
that represents the channel distribution of the whole BS cell. On the one
hand, this generative model allows to leverage prior information about the
distribution of the channels in the whole BS cell that enhances the channel
estimation performance [4]. On the other hand, it can be adapted to
imperfections in the training data as discussed above which is motivated by
model-based insights. This is a fundamental difference to common learning-
based estimators that are trained to learn a regression mapping between pilot
observations and channel estimates, which inherently rely on perfect CSI
samples, e.g., [3, 2, 1]. Thus, the proposed GMM approach allows for a more
robust solution with respect to various imperfections in the training data
because of its ability to adapt the training procedure accordingly and
introduce structural regularization, cf. Section III-B.

Algorithm 1 Adapted EM with Structured Covariances.

0: $\mathcal{Y}$, $K$, ${\bm{A}}$, ${\bm{C}}_{{\bm{n}}}$,
$\\{{\bm{\mu}}_{k}^{(1)},{\bm{C}}_{k}^{(1)},\pi_{k}^{(1)}\\}_{k=1}^{K}$,
${\bm{Q}}$, $i=1$, $i_{\text{max}}$

1: Get selection matrix $\bar{{\bm{A}}}$ such that
${\bm{A}}\bar{{\bm{A}}}^{\operatorname{T}}={\bm{0}}$

2: Initialize
$\\{{\bm{c}}_{k}^{(1)}\leftarrow{\bm{Q}}{\bm{C}}_{k}^{(1)}{\bm{Q}}^{\operatorname{H}}\\}_{k=1}^{K}$

3: while $i<i_{\text{max}}$ and convergence criterion not met do

4:
$\\{{\bm{C}}_{{\bm{y}},k}^{(i)}\leftarrow{\bm{A}}{\bm{C}}_{k}^{(i)}{\bm{A}}^{\operatorname{T}}+{\bm{C}}_{{\bm{n}}}\\}_{k=1}^{K}$

5: for $k=1$ to $K$ do

6: for $\ell=1$ to $L$ do

7:
$\gamma_{k,\ell}\leftarrow\frac{\pi_{k}^{(i)}\mathcal{N}_{\mathbb{C}}({\bm{y}}_{\ell};{\bm{A}}{\bm{\mu}}_{k}^{(i)},{\bm{C}}_{{\bm{y}},k}^{(i)})}{\sum_{j=1}^{K}\pi_{j}^{(i)}\mathcal{N}_{\mathbb{C}}({\bm{y}}_{\ell};{\bm{A}}{\bm{\mu}}_{j}^{(i)},{\bm{C}}_{{\bm{y}},j}^{(i)})}$
{E-step}

8:
$\bar{{\bm{y}}}_{k,\ell}\leftarrow\bar{{\bm{A}}}{\bm{\mu}}_{k}^{(i)}+\bar{{\bm{A}}}{\bm{C}}_{k}^{(i)}{\bm{A}}^{\operatorname{T}}{\bm{C}}_{{\bm{y}},k}^{(i),-1}({\bm{y}}_{\ell}-{\bm{A}}{\bm{\mu}}_{k}^{(i)})$

9:
$\hat{{\bm{y}}}_{k,\ell}\leftarrow{\bm{A}}^{\operatorname{T}}{\bm{y}}_{\ell}+\bar{{\bm{A}}}^{\operatorname{T}}\bar{{\bm{y}}}_{k,\ell}$
{interpolated sample}

10: end

11: $N_{k}\leftarrow\sum_{\ell=1}^{L}\gamma_{k,\ell}$

12: $\pi_{k}^{(i+1)}\leftarrow\frac{N_{k}}{L}$ {mixing coefficient update}

13:
${\bm{\mu}}_{k}^{(i+1)}\leftarrow\frac{1}{N_{k}}\sum_{\ell=1}^{L}\gamma_{k,\ell}\hat{{\bm{y}}}_{k,\ell}$
{mean update}

14:
${\bm{\Sigma}}_{k}\leftarrow\bar{{\bm{A}}}{\bm{C}}_{k}^{(i)}\bar{{\bm{A}}}^{\operatorname{T}}-\bar{{\bm{A}}}{\bm{C}}_{k}^{(i)}{\bm{A}}^{\operatorname{T}}{\bm{C}}_{{\bm{y}},k}^{(i),-1}{\bm{A}}{\bm{C}}_{k}^{(i)}\bar{{\bm{A}}}^{\operatorname{T}}$

15:
$\hat{{\bm{y}}}_{k,\ell}\leftarrow\hat{{\bm{y}}}_{k,\ell}-{\bm{\mu}}_{k}^{(i+1)}$

16:
${\bm{C}}_{{\bm{y}},k}^{(i+1)}\leftarrow\frac{1}{N_{k}}\sum\limits_{\ell=1}^{L}\gamma_{k,\ell}\hat{{\bm{y}}}_{k,\ell}\hat{{\bm{y}}}_{k,\ell}^{\operatorname{H}}+\bar{{\bm{A}}}^{\operatorname{T}}{\bm{\Sigma}}_{k}\bar{{\bm{A}}}$

17:
${\bm{V}}_{k},~{}{\bm{\xi}}_{k}\leftarrow\operatorname{EVD}({\bm{C}}_{{\bm{y}},k}^{(i+1)}-{\bm{C}}_{{\bm{n}}})$

18: ${\bm{\xi}}_{k}^{\text{PSD}}\leftarrow\max({\bm{0}},{\bm{\xi}}_{k})$
{elementwise max}

19:
${\bm{C}}_{k}^{(i+1)}\leftarrow{\bm{V}}_{k}\operatorname{diag}({\bm{\xi}}_{k}^{\text{PSD}}){\bm{V}}_{k}^{\operatorname{H}}$

20:
${\bm{\Theta}}_{k}\leftarrow{\bm{Q}}\left({\bm{C}}_{k}^{(i),-1}{\bm{C}}_{k}^{(i+1)}{\bm{C}}_{k}^{(i),-1}-{\bm{C}}_{k}^{(i),-1}\right){\bm{Q}}^{\operatorname{H}}$

21:
${\bm{c}}_{k}^{(i+1)}\leftarrow{\bm{c}}_{k}^{(i)}+\operatorname{diag}\left(\operatorname{diag}({\bm{c}}_{k}^{(i)}){\bm{\Theta}}_{k}\operatorname{diag}({\bm{c}}_{k}^{(i)})\right)$

22:
${\bm{C}}_{k}^{(i+1)}\leftarrow{\bm{Q}}^{\operatorname{H}}\operatorname{diag}({\bm{c}}_{k}^{(i+1)}){\bm{Q}}$
{covariance update}

23: end

24: $i\leftarrow i+1$

25: end

26: return
$\\{{\bm{\mu}}_{k}^{(i)},{\bm{C}}_{k}^{(i)},{\bm{\pi}}_{k}^{(i)}\\}_{k=1}^{K}$

An important advantage of the GMM estimator is that only the training phase
must be changed in contrast to [4] where perfect training CSI is assumed to be
available. Thus, the online complexity and the number of parameters of the
estimator do not change. Furthermore, the GMM is universal and requires
further adaptation only in case of changing parameters of the propagation
environment which can be conveniently tracked with the frequently received
pilots at the BS. Additionally, it has to be trained only once for a given SNR
and can then be applied to any other SNR value, which is in contrast to
learning-based estimators where the dependency on the SNR is crucial, e.g.,
[2, 3]. Fortunately, since the proposed covariance updates in Theorem 1 and
Theorem 2 are the maximum likelihood solutions, all favorable properties of
the EM algorithm, such as a monotonically increasing likelihood, are
preserved.

##  IV Simulation Results

We present numerical results to evaluate the proposed method for two different
channel models in comparison to state-of-the-art estimators. We set
${\operatorname{E}}[\|{\bm{h}}\|^{2}]=N$ such that the
$\text{SNR}=1/\sigma^{2}$. We choose $K=64$ GMM components which is
practically reasonable [4]. We utilize $L=10^{5}$ training samples for all
data-based approaches, and evaluate on $10^{4}$ test samples of different MTs
that are not part of the training data.

###  IV-A Spatial Channel Model

We first evaluate the channel estimation performance for the spatial model
from Section II-A with
${\bm{A}}=\operatorname{\mathbf{I}}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}N}}$.
The GMM estimator which is based on perfect CSI is denoted by “GMM
$\mathcal{H}$”, whereas the GMM estimator that is naively used without any
modifications on training data from $\mathcal{Y}$ is denoted by “GMM
mismatch”. The proposed approach with adapted training from Section III-B is
denoted by “GMM $\mathcal{Y}$”. We analyze the performance in comparison with
the following baseline estimators. The curve labeled “genie” represents the
utopian MMSE estimator that has full knowledge of ${\bm{C}}_{{\bm{\delta}}}$
for each sample. We further evaluate the least squares (LS) solution
$\hat{{\bm{h}}}_{\text{LS}}={\bm{y}}$. We include the convolutional neural
network (CNN) estimator from [3], that learns a regression mapping, labeled
“CNN”. The CNN consists of two layers with rectified linear unit (ReLU)
activation and we use the truncated DFT matrix for the input/output transform,
cf. [3]. To achieve a fair comparison, the CNN estimator is trained on data
from $\mathcal{Y}$, which clearly introduces an unavoidable mismatch in the
learning phase.

In Fig. 1, the channel estimation performance of the above discussed
estimators is shown for $N=128$ BS antennas and one propagation cluster,
consisting of multiple sub-paths. The depicted SNR is the same in both
training and evaluation, i.e., a different set of training data $\mathcal{Y}$
is given for each SNR value, which mimics a realistic situation. The GMM
estimator with perfect CSI $\mathcal{H}$ performs very close to the genie-MMSE
estimator which is in accordance with the findings in [4]. If the GMM
estimator is fitted naively with data from $\mathcal{Y}$ without modifications
to the training procedure, this leads to a severe performance loss of about
$5$dB for the whole SNR range. However, with the proposed modifications, a
performance close to the perfect CSI case, and thus to the genie-MMSE
estimator, is possible, purely based on training data from $\mathcal{Y}$. The
CNN estimator has a substantial performance loss as compared to the proposed
estimator, especially in the high SNR regime.

$-10$$-5$$0$$5$$10$$15$$20$$25$$30$$10^{0}$$10^{-1}$$10^{-2}$$10^{-3}$$10^{-4}$SNR
[dB]Normalized MSEgenieLSCNNGMM $\mathcal{H}$GMM mismatchGMM $\mathcal{Y}$
Figure 1: Spatial channel model from Section II-A with $N=128$ BS antennas and
one propagation cluster. The SNR is the same for training and evaluation.

###  IV-B Doubly-Selective Fading Channel Model

We consider a typical OFDM frame structure of $N_{c}=12$ carriers having a
spacing of $15$kHz, yielding a bandwidth of $180$kHz, and $N_{t}=14$ time
slots, cf. Section II-B. In both the training and the test dataset we consider
MTs which move with a fixed but random velocity between three and $130$km/h,
i.e., each MT in the datasets has a velocity drawn from the uniform
distribution $v\sim\mathcal{U}(3,130)$km/h. This should represent a practical
BS environment where different MTs are moving with a different velocity. We
evaluate once again the GMM estimator with perfect CSI from $\mathcal{H}$ and
a naive version where the GMM is fitted with linearly interpolated data from
$\mathcal{Y}$, labeled “GMM lin-int”. The proposed approach with adapted
training from Section III-B and linear interpolation as initialization is
denoted by “GMM $\mathcal{Y}$”. In the OFDM case, we can impose structural
features to the covariances as regularization on the adapted approach. In this
letter, we focus on block-Toeplitz matrices as discussed in [9], labeled as
“GMM $\mathcal{Y}$ toep”. We compare our proposed methods with the typically
used linear interpolator, labeled “lin-int”, and with a linear MMSE estimator
based on a cell-wide global sample covariance computed with linearly
interpolated data from $\mathcal{Y}$, labeled as “samp-cov lin-int”. In this
case the genie-MMSE estimator is not computable since no covariance statistics
are provided from the simulator [6]. We also compare with the ML-based
ChannelNet from [2] (we adopt the same architecture and hyperparameters) where
the training data are comprised of linearly interpolated samples from
$\mathcal{Y}$ for a fair comparison.

$-10$$-5$$0$$5$$10$$15$$20$$25$$30$$10^{0}$$10^{-1}$$10^{-2}$$10^{-3}$$10^{-4}$SNR
[dB]Normalized MSElin-intGMM $\mathcal{H}$samp-cov lin-intGMM lin-intGMM
$\mathcal{Y}$GMM $\mathcal{Y}$ toepChannelNet Figure 2: OFDM model from
Section II-B with $N_{c}=12$, $N_{t}=14$, and $N_{p}=18$. The SNR is the same
for training and evaluation. $5$$8$$13$$18$$25$$32$$10^{-2}$$10^{-1}$Pilots
$N_{p}$Normalized MSElin-intGMM $\mathcal{H}$samp-cov lin-intGMM lin-intGMM
$\mathcal{Y}$GMM $\mathcal{Y}$ toepChannelNet Figure 3: OFDM model from
Section II-B with $N_{c}=12$, $N_{t}=14$, and $\text{SNR}=15$dB. The SNR is
the same for training and evaluation.

In Fig. 2, the performances of the above discussed estimators are evaluated
for $N_{p}=18$ pilots over different SNR values. It can be seen that the
ChannelNet has no performance improvement over linear interpolation, which is
a consequence of the imperfect training data. The sample covariance and the
naive GMM approach based on linearly interpolated training data exhibit some
performance improvement over linear interpolation in the low and medium SNR
range, whereas for high SNRs, an error floor because of the imperfections in
the training data can be observed. The adapted GMM version without structural
constraints is performing close to the perfect CSI case in the low SNR range,
but also seems to saturate for high SNRs. In this SNR-region, the block-
Toeplitz structured adapted GMM performs especially well, showing performance
gains over all other approaches even for high SNR values. This leads to the
conclusion that the structural regularization, which is possible with the GMM,
is of great value when having noisy training data with missing entries.

In Fig. 3, we show a similar setup as above for a fixed SNR of $15$dB over
varying numbers of pilots. Thereby, linear interpolation and the ChannelNet
perform worst with a saturation at a high error floor, which was similarly
observed before. It can be seen that the differences between the sample
covariance, naive GMM, and adapted (unconstrained) GMM become more distinctive
for higher numbers of pilots where the adapted version is outperforming them
and has a decreasing gap to the perfect CSI based GMM. This can be reasoned
with a dominating systematic error for sparse pilot allocations. In contrast,
the block-Toeplitz based adapted GMM shows a substantial performance gain over
the baselines even for low numbers of pilots, which is an effect of the
structural regularization. As expected, the difference between the
structurally unconstrained and constrained adapted GMM decreases for higher
numbers of pilots.

##  V Conclusion

In this letter, a GMM-based robust channel estimator was proposed which can be
applied in spatial and OFDM systems. Thereby, the training data are purely
comprised of sparsely allocated and noisy pilot observations, without perfect
CSI samples. Simulation results demonstrated that the performance of the
proposed adapted GMM estimator is close to the version which utilizes perfect
training CSI with the same online complexity and memory overhead.
Additionally, state-of-the-art baselines are outperformed. Based on these
findings, the superior robustness properties against imperfect training data
of the generative model-aided estimator in contrast to regression-based ML
approaches were discussed. In future work, one may also account for
imperfections due to interference which is not considered in this letter.

## References

  * [1] H. Ye, G. Y. Li, and B.-H. Juang, “Power of Deep Learning for Channel Estimation and Signal Detection in OFDM Systems,” _IEEE Wireless Commun. Lett._ , vol. 7, no. 1, pp. 114–117, 2018.
  * [2] M. Soltani, V. Pourahmadi, A. Mirzaei, and H. Sheikhzadeh, “Deep Learning-Based Channel Estimation,” _IEEE Commun. Lett._ , vol. 23, no. 4, pp. 652–655, 2019.
  * [3] D. Neumann, T. Wiese, and W. Utschick, “Learning the MMSE Channel Estimator,” _IEEE Trans. Signal Process._ , vol. 66, no. 11, pp. 2905–2917, Jun. 2018.
  * [4] M. Koller, B. Fesl, N. Turan, and W. Utschick, “An Asymptotically MSE-Optimal Estimator Based on Gaussian Mixture Models,” _IEEE Trans. Signal Process._ , vol. 70, pp. 4109–4123, 2022.
  * [5] M. Baur, B. Fesl, M. Koller, and W. Utschick, “Variational Autoencoder Leveraged MMSE Channel Estimation,” in _56th Asilomar Conf. Signals, Syst., Comput._ , 2022.
  * [6] S. Jaeckel, L. Raschkowski, K. Börner, and L. Thiele, “QuaDRiGa: A 3-D Multi-Cell Channel Model With Time Evolution for Enabling Virtual Field Trials,” _IEEE Trans. Antennas Propag._ , vol. 62, no. 6, pp. 3242–3256, 2014.
  * [7] J.-W. Choi and Y.-H. Lee, “Optimum Pilot Pattern for Channel Estimation in OFDM Systems,” _IEEE Trans. Wireless Commun._ , vol. 4, no. 5, pp. 2083–2088, Sep. 2005.
  * [8] C. M. Bishop, _Pattern Recognition and Machine Learning (Information Science and Statistics)_. Berlin, Heidelberg: Springer-Verlag, 2006.
  * [9] B. Fesl, M. Joham, S. Hu, M. Koller, N. Turan, and W. Utschick, “Channel Estimation based on Gaussian Mixture Models with Structured Covariances,” in _56th Asilomar Conf. Signals, Syst., Comput._ , 2022.
  * [10] D. Neumann, M. Joham, L. Weiland, and W. Utschick, “Low-Complexity Computation of LMMSE Channel Estimates in Massive MIMO,” in _Proc. 19th Int. ITG Workshop Smart Antennas_ , 2015.
  * [11] R. J. A. Little and D. B. Rubin, _Statistical Analysis with Missing Data_. John Wiley & Sons, Incorporated, 2002.