# Blind Demixing of Diffused Graph Signals

Fernando J. Iglesias Garcia, , Santiago Segarra, ,
Antonio G. Marques Work in this paper is supported by Spanish Fed. Grants
KLINILYCS TEC2016-75361-R, SPGraph PID2019-105032GB-I00, by the Grants
F661-MAPPING-UCI and F663-AAGNCS funded by the Comunidad de Madrid (CAM) and
King Juan Carlos University (URJC), and by the USA NSF under award
CCF-2008555. Fernando J. Iglesias Garcia and Antonio G. Marques are with the
Dept. of Signal Theory and Comms., King Juan Carlos Univ. Santiago Segarra is
with the Dept. of Electrical and Computer Eng., Rice University. Emails:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Preliminary results were published in the conference paper [1].

###### Abstract

Using graphs to model irregular information domains is an effective approach
to deal with some of the intricacies of contemporary (network) data. A key
aspect is how the data, represented as graph signals, depend on the topology
of the graph. Widely-used approaches assume that the observed signals can be
viewed as outputs of graph filters (i.e., polynomials of a matrix
representation of the graph) whose inputs have a particular structure.
Diffused graph signals, which correspond to an originally sparse (node-
localized) signal percolated through the graph via filtering, fall into this
class. In that context, this paper deals with the problem of jointly
identifying graph filters and separating their (sparse) input signals from a
mixture of diffused graph signals, thus generalizing to the graph signal
processing framework the classical blind demixing (blind source separation) of
temporal and spatial signals. We first consider the scenario where the
supporting graphs are different across the signals, providing a theorem for
demixing feasibility along with probabilistic bounds on successful recovery.
Additionally, an analysis of the degenerate problem of demixing with a single
graph is also presented. Numerical experiments with synthetic and real-world
graphs empirically illustrating the main theoretical findings close the paper.

###### Index Terms:

Graph signal processing, blind demixing, blind source separation, blind system
identification, graph filter, network diffusion process.

##  I Introduction

Modern interconnected systems and networks are growing more pervasive and
complex by the day. While this is rendering data more abundant, its structure
is more intricate than the one exhibited by classical datasets. Relying on
graphs to model irregular information domains has emerged as a successful and
tractable approach to deal with some of the challenges present in contemporary
datasets. This approach, which has a long history in statistics [2], has been
recently gaining traction in the machine learning [3, 4] and signal processing
communities [5]. In the latter case, graph signal processing (GSP) puts forth
a systematic approach where variables are modeled as nodes, the information
domain is modeled as a graph (so that the focus is on sparse heterogeneous
relations among pairs of variables), and the data itself is modeled as graph
signals. Graph signals are abundant in nature, including the congestion levels
in a telecommunication network, the electrical neurological activity in a
brain network, and the information being diffused across the agents of a
social network, to name some diverse examples. The key in GSP is to postulate
meaningful and tractable: i) models relating the properties of the data with
the properties of the graph and ii) graph-signal operators that depend on the
topology of the graph [6, 7, 8]. With these two ingredients, a number of GSP-
related problems have been investigated, ranging from basic sampling and
reconstruction schemes [9, 10, 11] to advanced deep graph neural network
architectures [12, 13, 14]. An area of particular interest is that of linear
graph filters and their application to model linear network dynamics and solve
inverse GSP-related problems. Mathematically, graph filters are graph-signal
operators that can be represented by a matrix that is a polynomial of either
the adjacency or Laplacian matrix of the graph [5, 6, 15]. Practically, they
can be implemented distributedly, account for a multi-resolution (multi-hop)
description of the graph, and can capture spatio-temporal dynamics taking
place across the nodes (agents) of a network [16, 17, 18]. As a result, there
has been a growing interest in using graph filters to postulate generative
models of graph signals [5, 19, 20]. The most widely-used approach is to model
the signal at hand as the output of a graph filter excited with a parsimonious
(e.g., sparse, white, or bandlimited) input graph signal. One can then
formulate a number of inverse problems where (an imperfect version of) the
output graph signal is given and, depending on the particular setup and the
variables that are known, the goal is to recover the filter, the input, the
parts of the output that are unknown or any combination thereof [20]. When the
input graph signal is sparse (i.e., there are only a few nodes -sources- with
a non-zero value), the associated output is referred to as a diffused graph
signal, with sampling, reconstruction, blind deconvolution and filter
identification of diffused graph signals having been investigated in detail
[21, 22, 23, 24].

Contributions. In this paper, we focus on blind demixing of graph signals. In
particular, the observed graph signal $\mathbf{y}$ is assumed to be the
combination (sum) of the diffused graph signals
$\\{\mathbf{y}_{\text{i}}\\}_{\text{i}=1}^{\text{R}}$, each of which is
modeled as the output of a graph filter with coefficients
$\mathbf{h}_{\text{i}}$ excited with an input graph signal
$\mathbf{x}_{\text{i}}$. The signal $\mathbf{x}_{\text{i}}$ is assumed to be
localized in a subset of the nodes of the ith graph. In this setting we seek
to demix the observed signal $\mathbf{y}$ by _jointly_ estimating the input
graph signals _as well as_ the filters coefficients. This extends the
classical problem of blind source separation or demixing of signals in the
time and spatial domains. Apart from its theoretical interest, identification
of diffused graph signals also bears practical relevance, since they have been
shown to be germane to, e.g., explaining the structure of real-world datasets
[25, 22] as well as understanding the dynamics of linear network diffusion
processes [21, 22, 5].

Closely related work. We generalize the problem of blind demixing to _graph_
filters and signals, providing _both_ algorithms and theoretical analysis.
Relative to [21], which deals with the blind deconvolution of a single graph
signal (diffused by a single graph filter), the demixing setup considered here
deals with multiple graph signals (each associated with its own graph filter
and own graph) with the goal of separating them from a mixture. On top of
considering different operating conditions, demixing is intrinsically more
challenging since it is a more severely undetermined problem. Recently, a
data-driven deep learning approach was proposed for the demixing of non-linear
network processes [26]. Though relevant in practice, the approach in [26]
requires access to training data and, more importantly, lacks an in-depth
theoretical analysis as the one presented here. Finally, relative to the
classical (non-graph) demixing setup in [27], the approach in our paper not
only can be used in applications dealing with networks and irregular
information domains, but also entails different technical modeling
assumptions, namely: (1) sparsity of the input signals; (2) generalization of
random Fourier matrices (in the classical setup) to random unitary matrices
(in the GSP setup); and (3) consideration of multiple signal supports (one per
signal).

Paper outline. Section II summarizes general GSP concepts related to graph
filters and linear diffusion network processes. Section III formulates the
problem of blind demixing graph signals and proposes an efficient convex
relaxation. Section IV provides analytical results on the demixing
performance, along with graph-specific parameters that affect the recovery
guarantees. Section V deals with the special case of demixing with a single
graph _common_ to all the mixing components. Numerical experiments
illustrating the merits of our approach and the practical relevance of the
theoretical analysis are presented in Section VI.111Notational conventions:
The entries of a matrix $\mathbf{X}$ and a (column) vector $\mathbf{x}$ are
denoted as $X_{ij}$ and $x_{i}$ or $[\mathbf{x}]_{i}$. Operators
$(\cdot)^{\text{T}}$, $(\cdot)^{\text{H}}$, $\mathbb{E}{}$, $\circ$, $\otimes$
and $\odot$ stand for matrix transpose, conjugate transpose (Hermitian),
expectation, Hadamard (entry-wise), Kronecker, and Khatri-Rao (column-wise
Kronecker) products, respectively. The complex conjugate of $x$ is denoted by
$\bar{x}$; $\operatorname{diag}(\mathbf{x})$ is a diagonal matrix whose
$(i,i)$th entry is equal to $x_{i}$; $\mathcal{C}$ denotes a generic numerical
constant; and $|\cdot|$ is used for the cardinality of a set, and the
magnitude of a scalar. The $N\times N$ identity matrix is represented by
$\mathbf{I}_{N}$, while $\mathbf{0}_{N}$ stands for the $N\times 1$ vector of
all zeros, and $\mathbf{0}_{N\times P}=\mathbf{0}_{N}\mathbf{0}_{P}^{T}$. The
norms $\|\mathbf{X}\|_{\infty}$ and $\|\mathbf{X}\|$ stand for the entrywise
largest absolute value and the largest singular value of $\mathbf{X}$,
respectively. The sum of the absolute values of the elements of a vector
$\mathbf{x}$ is denoted as the norm $\|\mathbf{x}\|_{1}$. In the case of
matrices, the $\|\mathbf{X}\|_{2,1}$ norm is the sum of the Euclidean norms
(i.e. $\|\cdot\|_{2}$) of the rows of $\mathbf{X}$. The number of non-zero
elements of a vector $\mathbf{x}$ is denoted as $\|\mathbf{x}\|_{0}$.
$\|\mathbf{X}\|_{\ast}$ denotes the nuclear norm of matrix $\mathbf{X}$ and
$\|\mathbf{X}\|_{2,0}$ its number of non-zero rows. For a linear operator
$\mathcal{X}$, $\|\mathcal{X}\|:=\sup_{\mathbf{Z}\neq
0}\|\mathcal{X}(\mathbf{Z})\|_{\mathrm{F}}/\|\mathbf{Z}\|_{\mathrm{F}}=\|\mathbf{X}\|$,
where $\mathbf{X}$ is the matrix representation of $\mathcal{X}$. Otherwise,
standard vector and matrix norm notation is used.

##  II Graph Filtering Fundamentals

Let $\mathcal{G}$ denote a directed graph with a set of nodes $\mathcal{N}$
(with cardinality N) and a set of links $\mathcal{E}$ such that if $i$ is
connected to $j$ then $(i,j)\in\mathcal{E}$. Local connectivity is captured by
the set $\mathcal{N}_{i}:=\\{j\;|(j,i)\in\mathcal{E}\\}$, which stands for the
(incoming) neighborhood of $i$. For any given $\mathcal{G}$, we define the
adjacency matrix $\mathbf{A}\in\mathbb{R}^{\text{N}\times\text{N}}$ as a
sparse matrix with non-zero elements $A_{ji}$ if and only if
$(i,j)\in\mathcal{E}$. The value of $A_{ji}$ captures the strength of the
connection from $i$ to $j$. The focus of the paper is on analyzing and
modeling (graph) signals defined on $\mathcal{N}$. These signals can be
represented as vectors
$\mathbf{x}=[x_{1},...,x_{\text{N}}]^{\text{T}}\in\mathbb{R}^{\text{N}}$,
where $x_{i}$ represents the value of the signal at node $i$.

###  II-A Graph-shift operator

Since the vector representation
$\mathbf{x}=[x_{1},...,x_{\text{N}}]^{\text{T}}$ does not account explicitly
for the network structure, GSP endows $\mathcal{G}$ with the _graph-shift
operator_ (GSO) $\mathbf{S}$ [15, 5]. The shift
$\mathbf{S}\in\mathbb{R}^{\text{N}\times\text{N}}$ is a matrix whose entry
$S_{ji}$ can be non-zero only if $i=j$ or if $(i,j)\in\mathcal{E}$. The
sparsity pattern of the matrix $\mathbf{S}$ captures the local structure of
$\mathcal{G}$, with no specific assumptions on the values of its non-zero
entries. The intuition behind $\mathbf{S}$ is to represent a linear
transformation that can be computed locally at the nodes of the graph. More
rigorously, if $\mathbf{y}$ is defined as $\mathbf{y}=\mathbf{S}\mathbf{x}$,
then node $i$ can compute $y_{i}$ as linear combination of the signal values
$x_{j}$ where $j\in\mathcal{N}_{i}$. Typical choices for $\mathbf{S}$ are the
adjacency matrix $\mathbf{A}$ [15] and the graph Laplacian [6]. We assume
henceforth that $\mathbf{S}$ is diagonalizable, so that
$\mathbf{S}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1}$ with
$\boldsymbol{\Lambda}\in\mathbb{C}^{\text{N}\times\text{N}}$ being diagonal.
In particular, $\mathbf{S}$ is guaranteed to be diagonalizable when it is
normal, i.e., when it satisfies
$\mathbf{S}\mathbf{S}^{\text{H}}=\mathbf{S}^{\text{H}}\mathbf{S}$ (this is one
of the core assumptions in Section IV). In that case, we have that
$\mathbf{V}$ is unitary, which implies
$\mathbf{V}^{-1}=\mathbf{V}^{\text{H}}$, and leads to the decomposition
$\mathbf{S}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{\text{H}}$.

###  II-B Graph filters and diffusion processes

The shift $\mathbf{S}$ can be used to define _linear graph-signal operators_
of the form

$\mathbf{H}:=\sum_{\text{l}=0}^{\text{L}-1}h_{\text{l}}\mathbf{S}^{\text{l}},$
(1)

which are called _graph filters_ [15, 17]. For a given input $\mathbf{x}$, the
output of the filter is simply $\mathbf{y}=\mathbf{H}\mathbf{x}$. The
coefficients of the filter are collected into
$\mathbf{h}:=[h_{0},\ldots,h_{\text{L}-1}]^{\text{T}}$, with $\text{L}-1$
denoting the filter degree. Graph filters are of particular interest because
they model network diffusion processes and can be implemented in a distributed
fashion [17, 18]. Note also that the Cayley-Hamilton theorem guarantees that
even infinite-horizon processes (i.e., those depending on arbitrarily large
powers of $\mathbf{S}$) can be rewritten as a filter of at most degree
$\text{N}-1$.

###  II-C Spectral representation

Leveraging the spectral decomposition of $\mathbf{S}$, graph filters and
signals can be represented in the frequency domain. To be precise, let us use
the eigenvectors of $\mathbf{S}$ to define the $\text{N}\times\text{N}$ matrix
$\mathbf{U}:=\mathbf{V}^{-1}$, and the eigenvalues of $\mathbf{S}$ to define
the $\text{N}\times\text{L}$ Vandermonde matrix $\mathbf{\Psi}$, where
$\Psi_{ij}:=(\Lambda_{ii})^{j-1}$. Using these conventions, the frequency
representations of a _signal_ $\mathbf{x}$ and of a _filter_ $\mathbf{h}$ are
defined as $\widehat{\mathbf{x}}:=\mathbf{U}\mathbf{x}$ and
$\widehat{\mathbf{h}}:=\mathbf{\Psi}\mathbf{h}$, respectively [5, 17].
Exploiting such representations, the output
$\mathbf{y}\\!=\\!\mathbf{H}\mathbf{x}$ of a graph filter in the frequency
domain is given by

$\widehat{\mathbf{y}}=\operatorname{diag}\big{(}\mathbf{\Psi}\mathbf{h}\big{)}\mathbf{U}\mathbf{x}=\operatorname{diag}\big{(}\widehat{\mathbf{h}}\big{)}\widehat{\mathbf{x}}=\widehat{\mathbf{h}}\circ\widehat{\mathbf{x}},$
(2)

with $\circ$ being the Hadamart product. Identity (2) is the counterpart of
the celebrated convolution theorem for temporal signals, and follows from
$\mathbf{H}=\mathbf{V}\big{(}\sum_{\text{l}=0}^{\text{L}-1}h_{\text{l}}\boldsymbol{\Lambda}^{\text{l}}\big{)}\mathbf{U}$
[cf. (1)] and
$\sum_{\text{l}=0}^{\text{L}-1}h_{\text{l}}\boldsymbol{\Lambda}^{l}=\operatorname{diag}{(\mathbf{\Psi}\mathbf{h})}$.
Having $\widehat{\mathbf{y}}=\widehat{\mathbf{h}}\circ\widehat{\mathbf{x}}$
implies that the dependence of $\widehat{\mathbf{y}}$ on the frequency
representations of $\mathbf{h}$ and $\mathbf{x}$ is bilinear, with no mixing
among frequencies (i.e, $[\widehat{\mathbf{y}}]_{k}$ depends only on
$[\widehat{\mathbf{h}}]_{k}[\widehat{\mathbf{x}}]_{k}$). In contrast, if one
wants to write $\widehat{\mathbf{y}}$ as a function of the entries of
$\mathbf{h}$ and $\mathbf{x}$, it follows from (2) that the expression is

$\widehat{\mathbf{y}}=(\mathbf{\Psi}_{\text{i}}^{\text{T}}\odot\mathbf{U}_{\text{i}}^{\text{T}})^{\text{T}}\operatorname{vec}(\mathbf{x}\mathbf{h}^{\text{T}}),$
(3)

with $\odot$ being the Khatri-Rao product; see, e.g, [21]. That is, the
dependency is still bilinear, but any entry of $\widehat{\mathbf{y}}$ depends
on all the bilinear terms $[\mathbf{h}]_{\text{l}}[\mathbf{x}]_{n}$.
Interestingly, while in the time domain $\boldsymbol{\Psi}=\mathbf{U}$, in GSP
the matrices representing the graph Fourier transform (GFT) for filters and
signals are in general not the same.

##  III Blind Demixing with Multiple Diffusing Graphs

Given an observation
$\mathbf{y}=\sum_{\text{i}=1}^{\text{R}}\mathbf{y}_{\text{i}}$, we aim to
identify the R individual signals under the assumption that each
$\mathbf{y}_{\text{i}}$ can be modeled as a diffused graph signal, that is,
$\mathbf{y}_{\text{i}}$ is the output obtained by processing a _sparse_ input
signal $\mathbf{x}_{\text{i}}$ of _unknown_ support $\Omega_{\text{i}}$ (of
cardinality $|\Omega_{\text{i}}|=\text{S}_{\text{i}}$) with a graph filter of
unknown coefficients $\mathbf{h}_{\text{i}}$ and known GSO
$\mathbf{S}_{\text{i}}$. That is, by modeling each $\mathbf{y}_{\text{i}}$ as

$\mathbf{y}_{\text{i}}=\sum_{\text{l}=0}^{\text{L}-1}[\mathbf{h}_{\text{i}}]_{\text{l}}\mathbf{S}_{\text{i}}^{\text{l}}\mathbf{x}_{\text{i}},\;\;\text{where}\;\;[\mathbf{x}_{\text{i}}]_{n}=0\;\;n\notin\Omega_{\text{i}},$
(4)

we look at the problem of recovering
$\\{\mathbf{x}_{\text{i}},\mathbf{h}_{\text{i}},\mathbf{y}_{\text{i}}\\}_{\text{i}=1}^{\text{R}}$
from the observed $\mathbf{y}$.

With $\mathbf{\Psi}_{\text{i}}$, $\mathbf{V}_{\text{i}}$ and
$\mathbf{U}_{\text{i}}$ denoting the (inverse) GFT matrices associated with
the shift $\mathbf{S}_{\text{i}}$ (cf. Section II-A), and leveraging the
spectral representation of graph signals and filters (cf. Section II-C), the
joint demixing and blind deconvolution problem can be written as

$\displaystyle\mathop{\text{find}}$
$\displaystyle\quad\\{\mathbf{x}_{\text{i}},\mathbf{h}_{\text{i}},\mathbf{y}_{\text{i}}\\}_{\text{i}=1}^{\text{R}}$
(5a) s. to
$\displaystyle\quad\mathbf{y}_{\text{i}}=\mathbf{V}_{\text{i}}\left(\mathbf{\Psi}_{\text{i}}^{\text{T}}\odot\mathbf{U}_{\text{i}}^{\text{T}}\right)^{\text{T}}\operatorname{vec}(\mathbf{x}_{\text{i}}\mathbf{h}_{\text{i}}^{\text{T}})\;\;\;\forall\text{i}$
(5b)
$\displaystyle\quad\mathbf{y}={\textstyle\sum_{\text{i}=1}^{\text{R}}}\mathbf{y}_{\text{i}},\quad\quad\|\mathbf{x}_{\text{i}}\|_{0}\leq\text{S}_{\text{i}}\;\;\;\forall\text{i}.$
(5c)

When solving the above optimization, we assume that we know: i) the GSOs
$\\{\mathbf{S}_{\text{i}}\\}_{\text{i}=1}^{\text{R}}$, so that
$\\{\mathbf{\Psi}_{\text{i}},\mathbf{V}_{\text{i}},\mathbf{U}_{\text{i}}\\}_{\text{i}=1}^{\text{R}}$
are available; and ii) the full vector of observations $\mathbf{y}$ (see
Remark 1 for partially observable setups).

To handle the demixing problem in (5), note first that the constraint (5b)
implies that signals $\mathbf{y}_{\text{i}}$ are fully determined by
$\mathbf{x}_{\text{i}}$ and $\mathbf{h}_{\text{i}}$. Hence, (5b) can be
combined with the pertinent constraint in (5c), so that
$\\{\mathbf{y}_{\text{i}}\\}_{\text{i}=1}^{\text{R}}$ are no longer explicit
optimization variables, but obtained using (4) after the optimization over
$\mathbf{x}_{\text{i}}$ and $\mathbf{h}_{\text{i}}$ has been solved. Although
this reduces considerably the dimension of the feasible set, the optimization
is still challenging for a number of reasons. Firstly, even if the support of
the inputs were known, for the problem to be well-posed, the number of non-
trivial unknowns
$\sum_{\text{i}=1}^{\text{R}}(\text{S}_{\text{i}}+\text{L}_{\text{i}})$ needs
to be less than the number of observations N. Secondly, each of the terms
$\mathbf{x}_{\text{i}}\mathbf{h}_{\text{i}}^{\text{T}}$ is bilinear, which
introduces a source of non-convexity and gives rise to an inherent scaling
ambiguity. Thirdly, the presence of the $\ell_{0}$ norm renders the
optimization NP-hard.

The first step to design a tractable relaxation of (5), is to lift the problem
by defining the $\text{N}\\!\times\\!\text{L}_{\text{i}}$ rank-one matrices
$\mathbf{Z}_{\text{i}}\\!=\\!\mathbf{x}_{\text{i}}\mathbf{h}_{\text{i}}^{\text{T}}$,
together with the $\text{N}\\!\times\\!\text{N}\text{L}_{\text{i}}$ transfer
matrices
$\mathbf{M}_{\text{i}}=\mathbf{V}_{\text{i}}\left(\mathbf{\Psi}_{\text{i}}^{\text{T}}\odot\mathbf{U}_{\text{i}}^{\text{T}}\right)^{\text{T}}$.
With these notational conventions, and combining the linear constraints in
(5b)-(5c) into a single one, the demixing problem can be equivalently written
as

$\displaystyle\\{\mathbf{Z}_{\text{i}}^{(*)}\\}_{\text{i}=1}^{\text{R}}$
$\displaystyle=\text{find}\quad\\{\mathbf{Z}_{\text{i}}\\}_{\text{i}=1}^{\text{R}}$
(6) $\displaystyle\text{s. to}\,\,\,\,\mathbf{y}$
$\displaystyle\\!=\\!{\textstyle\sum_{\text{i}=1}^{\text{R}}}\mathbf{M}_{\text{i}}\text{vec}\\!\left(\mathbf{Z}_{\text{i}}\right),\;\,\text{rank}(\mathbf{Z}_{\text{i}})\\!=\\!1,\;\,\|\mathbf{Z}_{\text{i}}\|_{2,0}\leq\text{S}_{\text{i}}.$

Note that $\|\mathbf{Z}_{\text{i}}\|_{2,0}$, defined as the number of non-zero
rows of $\mathbf{Z}_{\text{i}}$, mimics the role of
$\|\mathbf{x}_{\text{i}}\|_{0}$ in (5). The solution vectors
$\mathbf{x}_{\text{i}}^{(*)}$ and $\mathbf{h}_{\text{i}}^{(*)}$ are given by
the scaled versions of the right and left principal singular vectors of the
rank-one matrix $\mathbf{Z}_{\text{i}}^{(*)}$. Although the equivalent
formulation (6) is still difficult to solve, it leads to a convex relaxation
which is described in the ensuing section, after the following remark.

######  Remark 1 (Robust demixing)

The demixing problem – in particular constraint (5b) – can be rewritten to
account for noisy and incomplete observations. If the values of $\mathbf{y}$
are observed only at a subset of M nodes, both sides of (5b) must be left
multiplied by the sampling matrix
$\mathbf{B}_{\mathcal{M}}\\!\in\\!\\{0,1\\}^{\text{M}\times\text{N}}$ whose
rows correspond to canonical basis vectors identifying the indexes
$\mathcal{M}=\\{n_{m}\\}_{m=1}^{\text{M}}$ of the observed entries. Presence
of noise in the observations would require replacing the equality in (5b) with
a term penalizing the difference between the two sides of (5b), with the
penalty depending on the specific type of noise.

###  III-A Convex relaxation and algorithmic approach

The _feasibility_ problem in (6) is reformulated as a rank-minimization
problem where, first, the objective is augmented with $\ell_{2,0}$ constraints
and, then, tractable surrogates of the (non-convex) rank and the $\ell_{2,0}$
norm are used. Specifically, consider replacing the rank with the nuclear
norm, and the $\ell_{2,0}$ norm with the $\ell_{2,1}$ norm [28, 29]. The
resultant optimization is

$\displaystyle\\{\mathbf{Z}_{\text{i}}^{(*)}\\}_{\text{i}=1}^{\text{R}}$
$\displaystyle=\operatorname*{arg\,min}_{\\{\mathbf{Z}_{\text{i}}\\}_{\text{i}=1}^{\text{R}}}\quad\textstyle\sum_{\text{i}=1}^{\text{R}}\eta_{\text{i}}\|\mathbf{Z}_{\text{i}}\|_{*}+\sum_{\text{i}=1}^{\text{R}}\beta_{\text{i}}\|\mathbf{Z}_{\text{i}}\|_{2,1}$
(7) $\displaystyle\text{s.
to}\quad\mathbf{y}={\textstyle\sum_{\text{i}=1}^{\text{R}}\mathbf{M}_{\text{i}}}\text{vec}\left(\mathbf{Z}_{\text{i}}\right),$

where $\\{\eta_{\text{i}},\beta_{\text{i}}\\}_{\text{i}=1}^{\text{R}}$ are
tuning constants that, if no prior information is available, can be set to
$\eta_{\text{i}}=1$ and $\beta_{\text{i}}=\beta$ for all i. The optimization
in (7) is convex, so that it can be handled with generic off-the-shelf
solvers. When computational complexity is a concern, one can alternatively
develop efficient algorithms tailored to the structure of (7), see, e.g., [30]
for a related problem.

Recoverability conditions can be found by ignoring the low-rank nature of all
the $\mathbf{Z}_{\text{i}}$ and focusing on the fact that the concatenation of
the vectorized forms of all these matrices
$\mathbf{z}=\mathrm{vec}([\mathbf{Z}_{1},\ldots,\mathbf{Z}_{\text{R}}])$ is
effectively a sparse vector. In this way, one can think of recovering a sparse
vector $\mathbf{z}$ satisfying
$\mathbf{y}=[\mathbf{M}_{1},\ldots,\mathbf{M}_{\text{R}}]\,\mathbf{z}$ and
derive recoverability conditions in this basic setting [31]. A more
sophisticated analysis accounting for the structure of the linear operators
$\mathbf{M}_{\text{i}}$ and the low-rank nature of $\mathbf{Z}_{\text{i}}$ is
our main theoretical contribution, which is presented in Section IV.

##  IV Exact Blind Demixing Via Nuclear Norm Minimization

This section presents our main theoretical result, namely, under a few
technical conditions (to be outlined next) all the filters
$\mathbf{h}_{\text{i}}$ and sparse input signals $\mathbf{x}_{\text{i}}$ in
(5) can be exactly recovered with high probability by solving a simplified
version of the convex problem (7) where we only keep the nuclear norm in the
objective, i.e.,

$\text{minimize}\sum_{\text{i}=1}^{\text{R}}\|\mathbf{Z}_{\text{i}}\|_{*}\quad\text{subject
to}\quad\mathbf{y}=\sum_{\text{i}=1}^{\text{R}}\mathbf{M}_{\text{i}}\text{vec}\left(\mathbf{Z}_{\text{i}}\right).$
(8)

To prove this result, we introduce _concentration functions_ to characterize
the model matrices. In [21], given an arbitrary matrix
$\mathbf{A}\in\mathbb{C}^{M\times N}$ and a positive integer $k\leq N$, the
function $\rho$ is defined as

$\rho_{\mathbf{A}}^{(k)}:=\max_{\ell\in\\{1,2,\dots,M\\}}\max_{\Omega\in\Omega_{k}^{N}}\|\mathbf{a}_{\ell_{\Omega}}\|_{2}^{2},$
(9)

where $\Omega_{k}^{N}$ represents the set of all $k$-subsets of
$\\{1,\dots,N\\}$ and $\mathbf{a}_{\ell_{\Omega}}$ the orthogonal projection
of the $\ell$th row of $\mathbf{A}$ onto the index set $\Omega$. Building on
this idea, let us now consider a second matrix
$\mathbf{B}\in\mathbb{C}^{M\times N}$ and, with positive integers $k_{1}$ and
$k_{2}$, define the function $\kappa$ as

$\kappa_{\mathbf{A},\mathbf{B}}^{(k_{1},k_{2})}:=\max_{\ell\in\\{1,2,\dots,M\\}}\max_{(\Omega_{1},\Omega_{2})\in\Omega_{k_{1}}^{N}\times\Omega_{k_{2}}^{N}}\|\mathbf{a}_{{\ell}_{\Omega_{1}}}\
\mathbf{b}_{{\ell}_{\Omega_{2}}}^{\text{H}}\|.$ (10)

In order to relate $\rho$ and $\kappa$ recall that the matrix norm
$\|\mathbf{A}\|$ is induced from the vector norm
$\sup_{\mathbf{x}}\|\mathbf{A}\mathbf{x}\|_{2}$ where $\|\mathbf{x}\|_{2}=1$.
If the matrix under consideration is of rank one, as in (10), it can be
written as $\mathbf{A}=\mathbf{a}\mathbf{b}^{\text{H}}$, whose norm is given
in terms of $\|\mathbf{a}\mathbf{b}^{\text{H}}\mathbf{x}\|_{2}$. Organizing
the terms,
$\|\mathbf{a}\mathbf{b}^{\text{H}}\mathbf{x}\|_{2}=|\mathbf{b}^{\text{H}}\mathbf{x}|\
\|\mathbf{a}\|_{2}$, which is maximum when $\mathbf{x}$ is parallel to
$\mathbf{b}$ and, therefore,
$\|\mathbf{a}\mathbf{b}^{\text{H}}\|=\|\mathbf{a}\|_{2}\|\mathbf{b}\|_{2}$.
Thus, we can relate $\rho$ with $\kappa$ as

$\kappa_{\mathbf{A},\mathbf{B}}^{(k_{1},k_{2})}\leq\sqrt{\rho_{\mathbf{A}}^{(k_{1})}}\sqrt{\rho_{\mathbf{B}}^{(k_{2})}}.$
(11)

The assumptions put in place in our analysis are outlined in Section IV-A,
after which the main result is formally stated in Section IV-B and proved in
Appendix -A.

###  IV-A Assumptions and scope of the analysis

We start by introducing two assumptions that render the analysis more
tractable.

Assumption 1 (AS1) The GSOs $\mathbf{S}_{\text{i}},\
\text{i}\\!\in\\!\\{1,2,\dots,\text{R}\\}$, are normal, i.e., they satisfy
$\mathbf{S}_{\text{i}}\mathbf{S}_{\text{i}}^{\text{H}}=\mathbf{S}_{\text{i}}^{\text{H}}\mathbf{S}_{\text{i}}$,
and for each $\mathbf{S}_{\text{i}}$ its eigenvalues are distinct.

Assumption 2 (AS2) The observation $\mathbf{y}$ is a mix of filtered graph
signals and its frequency representation adheres to the model

$\widehat{\boldsymbol{\nu}}=\sum_{\text{i}=1}^{\text{R}}\operatorname{diag}(\mathbf{\Psi}_{\text{i}}\mathbf{h}_{\text{i}})\widetilde{\mathbf{U}}_{\text{i}}\mathbf{x}_{\text{i}},$
(12)

where $\widetilde{\mathbf{U}}_{\text{i}}$ is a random $\text{N}\times\text{N}$
matrix obtained by concatenating N rows sampled independently and uniformly
with replacement from $\mathbf{U}_{\text{i}}$.

Under (AS1) the matrices $\mathbf{V}_{\text{i}}$ are unitary, which implies
$\mathbf{U}_{\text{i}}:=\mathbf{V}_{\text{i}}^{-1}=\mathbf{V}_{\text{i}}^{\text{H}}$
and leads to the decomposition
$\mathbf{S}_{\text{i}}=\mathbf{V}_{\text{i}}\boldsymbol{\Lambda}_{\text{i}}\mathbf{V}_{\text{i}}^{\text{H}}$.
Normality is for instance satisfied when $\mathcal{G}$ is undirected and the
GSO is chosen to be the adjacency matrix or the graph Laplacian. Furthermore,
that each $\mathbf{S}_{\text{i}}$ has distinct eigenvalues ensures that the
Vandermonde matrices $\mathbf{\Psi}_{\text{i}}$ are full-rank independently of
$\text{L}_{\text{i}}$, which is required for uniqueness.

Concerning the probabilistic model for the observations in (AS2), this class
of models is customary when establishing recovery guarantees of sparse signal
reconstruction and low-rank matrix completion in the context of, e.g.,
compressive sensing [31], or even classical demixing of temporal signals [27,
32].

###  IV-B Main theoretical result

Consider R $\text{S}_{\text{i}}$-sparse graph signals
$\mathbf{x}_{\text{i}}\\!\in\\!\mathbb{R}^{\text{N}},\text{i}\\!\in\\!\\{1,2,\dots,\text{R}\\}$,
each respectively diffused through its supporting graph using a filter with
coefficients $\mathbf{h}_{\text{i}}\in\mathbb{R}^{\text{L}_{\text{i}}}$, and
giving rise to an aggregated signal whose frequency representation adheres to
the model in (AS2). With the matrices
$\mathbf{U}_{\text{i}}\in\mathbb{C}^{\text{N}\times\text{N}}$ and
$\mathbf{\Psi}_{\text{i}}\in\mathbb{C}^{\text{N}\times\text{L}_{\text{i}}}$
denoting the GFT for signals and filters associated with the GSOs
$\mathbf{S}_{\text{i}}$, respectively, we assume without loss of generality
that the matrices $\mathbf{U}_{\text{i}}$ are normalized so that
$\mathbf{U}_{\text{i}}^{\text{H}}\mathbf{U}_{\text{i}}=\text{N}\mathbf{I}_{\text{N}}$.
In this setting, the following result holds.

###### Theorem 1

Under Assumptions (AS1) and (AS2), the unique solution to the convex program
(8) coincides with the original search (6) with probability at least
$1-\text{N}^{-\alpha+1}$ where
$1\leq\alpha\leq\min\\{\alpha_{1},\alpha_{2},\alpha_{3}\\}$,
$\rho_{\max}^{(k,j)}=\max\left\\{\rho_{\mathbf{U}_{k}}^{(S_{k})},\rho_{\mathbf{U}_{j}}^{(S_{j})}\right\\}$,
and

$\begin{split}&\alpha_{1}=\min_{\text{i}\in\\{1,2,\dots,\text{R}\\}}\frac{3}{128}\left(\rho_{\mathbf{\Psi}_{\text{i}}}^{(\text{L}_{\text{i}})}\rho_{\mathbf{U}_{\text{i}}}^{(\text{S}_{\text{i}})}\log\left(2\text{N}\text{L}_{\text{i}}\text{S}_{\text{i}}\right)\right)^{-1},\\\
&\alpha_{2}=\min_{j\neq
k}\frac{9}{32}\frac{\Big{(}\rho_{\mathbf{\Psi}_{k}}^{(L_{k})}\rho_{\mathbf{\Psi}_{j}}^{(L_{j})}\rho_{\max}^{(k,j)}+\frac{1}{\text{R}}\kappa_{\mathbf{U}_{k},\mathbf{U}_{j}}^{\left(\text{S}_{k},\text{S}_{j}\right)}\kappa_{\mathbf{\Psi}_{k},\mathbf{\Psi}_{j}}^{\left(\text{L}_{k},\text{L}_{j}\right)}\Big{)}^{-1}}{\text{R}^{2}\log(2\text{N})},\\\
&\alpha_{3}=\min_{\text{i}\in\\{1,2,\dots,\text{R}\\}}\frac{1}{\mathcal{C}}\left(\text{N}\text{R}\max\\{\mu_{\text{max}}^{2}\text{L}_{\text{i}},\mu_{h}^{2}\text{N}\\}\log^{2}{\text{N}}\right)^{-1}-\log{\text{R}}.\end{split}$
(13)

_Proof_ : See Appendix -A.

Note that the statement of the theorem leverages the definitions in (9)-(10)
and depends on the scalar quantities $\mu_{\text{h}}$ and $\mu_{\text{max}}$,
which can be interpreted as generalizations of (9)-(10) capturing the
interactions between $\mathbf{h}_{\text{i}}$ and $\mathbf{\Psi}_{\text{i}}$,
and are defined in detail in (23)-(24).

The main goal of the probabilistic recovery bounds in Theorem 1 is to
establish that a convex relaxation _can_ succeed in blind demixing of graph
signals and, in practice, extrapolate from the theoretical results a
_predictor_ of demixing performance based on the graph topologies. Although
admittedly loose, the bounds in the theorem motivate empirical predictors of
successful demixing as illustrated in Section VI.

Furthermore, we want to emphasize three differences between Theorem 1 and the
main result in [27]. First, Theorem 1 provides recovery guarantees for blind
demixing in arbitrary graphs [cf. (AS1)], whereas the results in [27] only
apply for the cases where $\widetilde{\mathbf{U}}_{\text{i}}$ are Gaussian
distributed. Second, [27] does not leverage sparsity in the mixing graph
signals and instead introduces auxiliary linear spaces to obtain probabilistic
bounds. Finally, we keep track of the constant terms throughout the proof,
which leads to bounds with finite values instead of order-wise asymptotic
results.

##  V Blind Demixing With A Single Diffusing Graph

So far, we have considered the case where each input $\mathbf{x}_{\text{i}}$
is diffused across its own graph characterized by its GSO
$\mathbf{S}_{\text{i}}$. The setting where all the signals are diffused in the
same graph and then combined leads to a more challenging blind demixing
problem. Notice that if $\mathbf{S}_{\text{i}}=\mathbf{S}$ for all i, then the
transfer matrix $\mathbf{M}_{\text{i}}=\mathbf{M}$ is equivalent for all the
inputs. As a result, the observation equation
$\mathbf{y}=\sum_{\text{i}=1}^{\text{R}}\mathbf{M}_{\text{i}}\text{vec}(\mathbf{Z}_{\text{i}})$
simplifies to
$\mathbf{y}=\mathbf{M}\text{vec}(\sum_{\text{i}=1}^{\text{R}}\\!\mathbf{Z}_{\text{i}})$,
which introduces a new source of ambiguity as explained in detail next.

###  V-A Single graph demixing ambiguity

First, consider the minimization form of the feasibility problem in (6), i.e.,
a non-relaxed predecessor of (7),

$\displaystyle\\{\mathbf{Z}_{\text{i}}^{(*)}\\}_{\text{i}=1}^{\text{R}}$
$\displaystyle=\operatorname*{arg\,min}_{\\{\mathbf{Z}_{\text{i}}\\}_{\text{i}=1}^{\text{R}}}\quad\textstyle\sum_{\text{i}=1}^{\text{R}}\eta_{\text{i}}\,\text{rank}(\mathbf{Z}_{\text{i}})+\sum_{\text{i}=1}^{\text{R}}\beta_{\text{i}}\|\mathbf{Z}_{\text{i}}\|_{2,0}$
s. to
$\displaystyle\mathbf{y}=\mathbf{M}\,\text{vec}\\!\left({\textstyle\sum_{\text{i}=1}^{\text{R}}\\!\mathbf{Z}_{\text{i}}}\right).$
(14)

To illustrate the ambiguity emerging in the single graph demixing problem,
consider the case where $\text{R}\\!=\\!2$ and, accordingly,
$\mathbf{Z}_{1}^{(*)}$ and $\mathbf{Z}_{2}^{(*)}$ denote the true rank-one
matrices. Next, use these matrices to define
$\mathbf{Z}_{1}^{\prime}=\mathbf{Z}_{1}^{(*)}+\mathbf{Z}_{2}^{(*)}$ and
$\mathbf{Z}_{2}^{\prime}=\mathbf{0}$. It readily follows that: 1) both
$\\{\mathbf{Z}_{1}^{(*)},\mathbf{Z}_{2}^{(*)}\\}$ and
$\\{\mathbf{Z}_{1}^{\prime},\mathbf{Z}_{2}^{\prime}\\}$ are in the feasible
set of (V-A); and 2) employing the triangle inequality:
$\text{rank}(\mathbf{Z}_{1}^{\prime})+\text{rank}(\mathbf{Z}_{2}^{\prime})=\text{rank}(\mathbf{Z}_{1}^{\prime})\leq\text{rank}(\mathbf{Z}_{1}^{(*)})+\text{rank}(\mathbf{Z}_{2}^{(*)})$.
In words, the true $\\{\mathbf{Z}_{1}^{(*)},\mathbf{Z}_{2}^{(*)}\\}$ involve a
(rank) cost larger than or equal to that of
$\\{\mathbf{Z}_{1}^{\prime},\mathbf{Z}_{2}^{\prime}\\}$. Furthermore, if
$\eta_{1}\neq\eta_{2}$ then a solution where either
$\mathbf{Z}_{1}^{\prime}=\mathbf{0}$ or $\mathbf{Z}_{2}^{\prime}=\mathbf{0}$
always achieves a smaller rank cost. A similar argument holds for the
$\ell_{2,0}$ norm that constitutes the objective in (V-A) together with the
rank just analyzed.

This demonstrates that the demixing problem with a single diffusing graph
presents additional challenges. With respect to the multiple-graph problem,
the individual $\mathbf{Z}_{\text{i}}$’s are in general non-identifiable
separately, albeit the combination
$\mathbf{Z}\\!=\\!\sum_{\text{i}=1}^{\text{R}}\\!\mathbf{Z}_{\text{i}}$ can be
recovered. Thus, we propose a two-step approach in which we first efficiently
obtain the true
$\mathbf{Z}^{(*)}\\!=\\!\sum_{\text{i}=1}^{\text{R}}\\!\mathbf{Z}^{(*)}_{\text{i}}$
and then outline conditions for which the constituent
$\mathbf{Z}^{(*)}_{\text{i}}$ matrices can be uniquely recovered.

###  V-B Convex relaxation and algorithmic approach

Following the strategy in Section III-A to design an efficient demixing
algorithm, the combinatorial norms in (V-A) are replaced with their convex
surrogates $\|\cdot\|_{*}$ and $\|\cdot\|_{2,1}$. This results in an instance
of (7) with all $\mathbf{M}_{\text{i}}=\mathbf{M}$ whose solutions
$\\{\mathbf{Z}_{\text{i}}^{(*)}\\}_{\text{i}=1}^{\text{R}}$ also suffer from
the linear ambiguity present in single graph demixing. In fact, given that the
nuclear as well as the $\ell_{2,1}$ norms are absolutely homogeneous (whereas
the rank and the $\ell_{2,0}$ penalties are not), the ambiguity is more severe
in the convex formulation. To see this, note that any solution of the form
$\mathbf{Z}_{\text{i}}\\!=\\!\theta_{\text{i}}\mathbf{Z}^{(*)}$ with
$\theta_{\text{i}}\\!\geq\\!0$ and
$\sum_{\text{i}=1}^{\text{R}}\theta_{\text{i}}\\!=\\!1$ is feasible and yields
the same value of $\sum_{\text{i}=1}^{\text{R}}\|\mathbf{Z}_{\text{i}}\|_{*}$.
Hence, the postulated optimization will have multiple solutions with the same
cost than that of the actual $\mathbf{Z}_{\text{i}}^{(*)}$ we aim to identify.

Another option consists in reformulating (7) based on
$\mathbf{Z}=\sum_{\text{i}=1}^{\text{R}}\mathbf{Z}_{\text{i}}$ as follows

$\mathbf{Z}^{(*)}=\operatorname*{arg\,min}_{\mathbf{Z}}\,\,\textstyle\|\mathbf{Z}\|_{*}+\beta\|\mathbf{Z}\|_{2,1}\quad\,\,\text{s.
to}\,\,\,\mathbf{y}=\textstyle\mathbf{M}\text{vec}\left(\mathbf{Z}\right).$
(15)

In terms of computational complexity, this reformulation presents a lower cost
because the number of optimization variables are effectively reduced by a
factor of R. Moreover, the formulation in (15) is similar to that of blind
deconvolution [21] (that is, a single source without demixing) with the
exception that the true $\mathbf{Z}^{(*)}$ to be recovered, albeit being low-
rank, is not necessarily rank-one. Still, the results in [21, Theorem 1] can
provide probabilistic guarantees for the recovery of the true
$\mathbf{Z}^{(*)}$ in terms of the matrices $\mathbf{\Psi}$ and $\mathbf{U}$
associated with the common diffusing graph.

(a)

(b)

(c)

Figure 1: Blind demixing performance in single-graph and multiple-graph
settings. (a) Demixing error as a function of the graph size N in a single
graph for four settings of S and R. (b) Rate of successful recovery between a
pair of graphs as a function of the _similarity_ parameter $\alpha$ between
them (when $\alpha=0$ the two graphs are statistically independent, when
$\alpha=1$ the graphs are the same). (c) Rate of successful recovery for
multiple brain graphs for several settings of L and R.

###  V-C Resolving the ambiguity

It is in general not feasible to recover
$\\{\mathbf{Z}_{\text{i}}^{(*)}\\}_{\text{i}=1}^{\text{R}}$ from
$\mathbf{Z}^{(*)}$ due to the single graph ambiguity already explained.
Consequently, we propose one way to ensure recoverability by means of imposing
orthogonality conditions on the filter coefficients and the input signals.
Specifically, we consider the following two separate assumptions:

(a1)

$\mathbf{x}_{\text{i}}^{\text{T}}\mathbf{x}_{\text{j}}=0$ or
$\widehat{\mathbf{x}}_{\text{i}}^{\text{T}}\widehat{\mathbf{x}}_{\text{j}}=0$
for all $\text{i}\neq\text{j}$.

(a2)

$\mathbf{h}_{\text{i}}^{\text{T}}\mathbf{h}_{\text{j}}=0$ or
$\widehat{\mathbf{h}}_{\text{i}}^{\text{T}}\widehat{\mathbf{h}}_{\text{j}}=0$
for all $\text{i}\neq\text{j}$.

These orthogonality conditions are satisfied, e.g., when the vectors are
sparse and their support does not overlap. In this way, two sparse inputs that
are non-zero in non-overlapping nodes satisfy the first requirement in (a1).
Similarly, two bandpass filters defined in non-overlapping bands satisfy the
second requirement in (a2).

Before formally stating recoverability in the single graph case, we define
auxiliary matrices to capture the orthogonality choices in (a1) and (a2).
Consider the matrix $\mathbf{T}_{\mathbf{x}}=\mathbf{I}$ if the first
requirement in (a1) holds and $\mathbf{T}_{\mathbf{x}}=\mathbf{U}$ if the
second requirement is true. Similarly, consider
$\mathbf{T}_{\mathbf{h}}=\mathbf{I}$ if the first requirement in (a2) is true
and $\mathbf{T}_{\mathbf{h}}=\mathbf{\Psi}$ if the second one holds. These
matrices allow us to state recoverability conditions and a final algorithmic
step for single graph demixing succinctly. That is, $\mathbf{T}_{\mathbf{x}}$
($\mathbf{T}_{\mathbf{h}}$) captures whether the input signals (filter taps)
are orthogonal in the node or spectral domains, rendering the remaining
analysis independent of the domain.

###### Proposition 1

If the input signals
$\\{\mathbf{x}_{\text{i}}^{(*)}\\}_{\text{i}=1}^{\text{R}}$ satisfy (a1), the
diffusing graph filters
$\\{\mathbf{h}_{\text{i}}^{(*)}\\}_{\text{i}=1}^{\text{R}}$ satisfy (a2) and
$\|\mathbf{T}_{\mathbf{x}}\mathbf{x}_{\text{i}}^{(*)}\|\|\mathbf{T}_{\mathbf{h}}\mathbf{h}_{\text{i}}^{(*)}\|\neq\|\mathbf{T}_{\mathbf{x}}\mathbf{x}_{\text{j}}^{(*)}\|\|\mathbf{T}_{\mathbf{h}}\mathbf{h}_{\text{j}}^{(*)}\|$
for all $\text{i}\neq\text{j}$; then
$\\{\mathbf{Z}^{(*)}_{\text{i}}\\}_{\text{i}=1}^{\text{R}}$ can be recovered
from $\mathbf{Z}^{(*)}$.

_Proof :_ First, form the matrix
$\mathbf{T}=\mathbf{T}_{\mathbf{x}}\mathbf{Z}^{(*)}\mathbf{T}_{\mathbf{h}}^{\text{T}}$
and factorize it using the singular value decomposition (SVD) as
$\mathbf{T}=\mathbf{L}\text{diag}(\boldsymbol{\sigma})\mathbf{R}^{\text{T}}$.
Then, the mixing components $\mathbf{Z}_{\text{i}}^{(*)}$ are estimated
individually by

$\breve{\mathbf{Z}}_{\text{i}}^{(*)}=\sigma_{\text{i}}(\mathbf{T}_{\mathbf{x}}^{\dagger})\mathbf{l}_{\text{i}}\mathbf{r}_{\text{i}}^{\text{T}}(\mathbf{T}_{\mathbf{h}}^{\dagger})^{\text{T}},$
(16)

where $\sigma_{\text{i}}$ is the ith largest singular value of $\mathbf{T}$,
$\mathbf{l}_{\text{i}}$ and $\mathbf{r}_{\text{i}}$ its associated left and
right singular vectors, and † denotes the pseudoinverse.

The remaining part is to show that
$\breve{\mathbf{Z}}_{\text{i}}^{(*)}=\mathbf{Z}_{\text{i}}^{(*)}$ for all i.
To this end, notice that we can write $\mathbf{T}$ as

$\mathbf{T}=\sum_{\text{i}}^{\text{R}}\|\mathbf{T}_{\mathbf{x}}\mathbf{x}_{\text{i}}^{(*)}\|\|\mathbf{T}_{\mathbf{h}}\mathbf{h}_{\text{i}}^{(*)}\|\frac{\mathbf{T}_{\mathbf{x}}\mathbf{x}_{\text{i}}^{(*)}}{\|\mathbf{T}_{\mathbf{x}}\mathbf{x}_{\text{i}}^{(*)}\|}\frac{(\mathbf{T}_{\mathbf{h}}\mathbf{h}_{\text{i}}^{(*)})^{\text{T}}}{\|\mathbf{T}_{\mathbf{h}}\mathbf{h}_{\text{i}}^{(*)}\|}.$
(17)

From (a1)-(a2), it follows that
${\mathbf{T}_{\mathbf{x}}\mathbf{x}_{\text{i}}^{(*)}}/{\|\mathbf{T}_{\mathbf{x}}\mathbf{x}_{\text{i}}^{(*)}\|}$
for all i forms an orthonormal basis and the same is true for
${\mathbf{T}_{\mathbf{h}}\mathbf{h}_{\text{i}}^{(*)}}/{\|\mathbf{T}_{\mathbf{h}}\mathbf{h}_{\text{i}}^{(*)}\|}$.
Moreover, it is assumed in the statement of the proposition that all
$\|\mathbf{T}_{\mathbf{x}}\mathbf{x}_{\text{i}}^{(*)}\|\|\mathbf{T}_{\mathbf{h}}\mathbf{h}_{\text{i}}^{(*)}\|$
are distinct. Thus, (17) must coincide with the SVD decomposition in that

$\mathbf{l}_{\text{i}}\\!=\\!\frac{\mathbf{T}_{\mathbf{x}}\mathbf{x}_{\text{i}}^{(*)}}{\|\mathbf{T}_{\mathbf{x}}\mathbf{x}_{\text{i}}^{(*)}\|},\,\,\,\,\,\mathbf{r}_{\text{i}}\\!=\\!\frac{\mathbf{T}_{\mathbf{h}}\mathbf{h}_{\text{i}}^{(*)}}{\|\mathbf{T}_{\mathbf{h}}\mathbf{h}_{\text{i}}^{(*)}\|},\,\,\,\,\,\sigma_{i}\\!=\\!\|\mathbf{T}_{\mathbf{x}}\mathbf{x}_{\text{i}}^{(*)}\|\|\mathbf{T}_{\mathbf{h}}\mathbf{h}_{\text{i}}^{(*)}\|.$
(18)

From (18) it readily follows that $\breve{\mathbf{Z}}_{\text{i}}^{(*)}$ in
(16) is equal to
$\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}$, as
sought. $\hfill\blacksquare$

Proposition 1 reveals that even for the highly challenging scenario of
multiple signals being diffused and aggregated in a single graph, demixing is
still possible under more stringent constraints. The analysis in this section
also hints at the fact that the similarity between graphs might be an
indicator of how hard the demixing problem is, with identical graphs being an
extreme setting. We revisit this point from an empirical viewpoint in the next
section.

Remark 2 (A priori information) Blind demixing benefits when prior knowledge
on
$\\{\mathbf{h}_{\text{i}}^{(*)},\mathbf{x}_{\text{i}}^{(*)}\\}_{\text{i}=1}^{\text{R}}$
is available. Even though this applies for both demixing cases with multiple
diffusing graphs or a single diffusing graph, a priori information is
particularly helpful in the single-graph scenario, where the inherent
ambiguity is larger. Two types of a priori information can be employed,
probabilistic and deterministic. Probabilistic information can be incorporated
by augmenting the objective with tractable (log) prior distributions.
Alternatively, deterministic knowledge about entries of the signals
$\mathbf{x}_{\text{i}}^{(*)}$ or the filter taps $\mathbf{h}_{\text{i}}^{(*)}$
can be incorporated in the form of linear equality constraints. In particular,
if $\\{\mathbf{h}_{\text{i}}^{(*)}\\}_{\text{i}=1}^{\text{R}}$ is known then
the structure of the optimization can be simplified. That is, the bilinearity
disappears and the resultant problem can be handled using classical sparse
signal reconstruction tools. One practical scenario occurs when
$\text{K}_{\text{i}}$ values of the signal $\mathbf{x}_{\text{i}}^{(*)}$ are
known. Since we optimize over $\mathbf{Z}_{\text{i}}$ instead of directly over
$\mathbf{x}_{\text{i}}$, exploiting this information requires considering the
constraints

$\mathbf{Z}_{\text{i},\ell}^{\text{T}}\mathbf{x}_{\text{i},\ell+1}^{(*)}=\mathbf{Z}_{\text{i},\ell+1}^{\text{T}}\mathbf{x}_{\text{i},\ell}^{(*)}\quad\forall\
{\ell\in\mathcal{K}_{\text{i}}},$ (19)

where $\mathcal{K}_{\text{i}}=\\{n_{\ell}\\}_{\ell=1}^{\text{K}_{\text{i}}}$
collects the indexes of the entries of $\mathbf{x}_{\text{i}}^{(*)}$ that are
known. Note that the introduction of this type of constraints still renders a
convex program.

##  VI Numerical Results

Five experiments are performed to illustrate and validate our proposed blind
demixing approach. In the first two, we assess and compare the effectiveness
of single-graph and multiple-graph demixing using synthetic random graphs. In
the third experiment, we analyze multiple-graph demixing using real-world
brain graphs (where nodes represent regions of the human brain and edges
anatomical connections) associated with different individuals. In the fourth
experiment, we analyze the recovery performance of single-graph demixing using
a real-world social network and, in the last experiment, we rely on some of
the graph-dependent parameters identified in Theorem 1 to predict (assess)
demixing performance in real-world graphs.

Throughout this section, we solve a collection of blind demixing problems,
where we change the type of graphs, the number of mixing components R, the
number of filter coefficients $\text{L}_{\text{i}}$, or the number of seeding
(non-zero signal) nodes $\text{S}_{\text{i}}$. The signal and filter-
coefficient estimates will be denoted by
$\\{\breve{\mathbf{x}}_{\text{i}},\breve{\mathbf{h}}_{\text{i}}\\}$. Unless
otherwise stated, we: set the GSO as
$\mathbf{S}_{\text{i}}=\mathbf{A}_{\text{i}}$, the adjacency matrix of the
graph $\mathcal{G}_{\text{i}}$; use the same values of
$\text{L}_{\text{i}}=\text{L}$ and $\text{S}_{\text{i}}=\text{S}$ across i;
and draw the ‘true’ vectors $\mathbf{x}^{(*)}_{\text{i}}$ and
$\mathbf{h}^{(*)}_{\text{i}}$ from standard multivariate Gaussian
distributions, normalizing them to have unit norm. In the single-graph
experiments, $\mathbf{x}^{(*)}_{\text{i}}$ and $\mathbf{h}^{(*)}_{\text{i}}$
are selected according to Proposition 1. That is, in addition to satisfying
(a1) and (a2) in Section V-C, the vectors are normalized to have the same
$\ell_{1}$ norm but different $\ell_{2}$ norm, so that it becomes possible to
resolve the ambiguity intrinsic to single-graph demixing. Figures of merit
derived from the demixing error (DE)

$\text{DE}:=\frac{1}{\text{R}}\sum_{\text{i}=1}^{\text{R}}\left\|\breve{\mathbf{x}}_{\text{i}}\breve{\mathbf{h}}_{\text{i}}^{\text{T}}-\mathbf{x}^{(*)}_{\text{i}}{\mathbf{h}^{(*)}_{\text{i}}}^{\text{T}}\right\|_{\text{F}},$
(20)

are adopted to assess demixing. These are the median DE across realizations
and the _rate of successful recovery_ defined to be the probability that
$\text{DE}<10^{-3}$. The results reported correspond to $1000$ demixing
realizations each of them solved using the log-det relaxation (see Section
III-A and [33, 23] for additional details).

Blind demixing in a single random graph. Figure 1(a) shows the results of
demixing with a single diffusing graph of varying N and $\text{L}=3$. Random
graphs are drawn according to the Erdős-Rényi model [34] with edge probability
$0.1$. Four settings with different values of R and S are presented. As
expected, the performance worsens for larger values of R and S. Interestingly,
demixing for $(\text{R}\\!=\\!3,\text{S}\\!=\\!6)$ is less successful than for
$(\text{R}\\!=\\!2,\text{S}\\!=\\!9)$. Since the number of signal coefficients
to estimate in these two settings is the same (RS), this suggests that the
recovery is more sensitive to R.

Effect of graph similarity on recovery. The second experiment explores
demixing with two Barabási-Albert graphs [35] of varying N (built following
the scale-free model from a seeding Erdős-Rényi graph of $\text{N}/10$ nodes
and edge probability $0.1$) and $\text{L}=3$. The two graphs are generated so
that $\alpha$ percent of their edges are the same. Note that $\alpha=1$ leads
to demixing with a single diffusing graph. Figure 1(b) shows that, as
expected, the multi-graph demixing strategy (cf. Section III) fails for
$\alpha=1$, although the performance does not considerably drop until
$\alpha>0.9$.

Blind demixing with multiple brain graphs. Figure 1(c) depicts the rate of
successful recovery using up to five different _brain graphs_ corresponding to
different people. The graphs belong to a dataset of six undirected graphs of
the human brain, consisting of $\text{N}=66$ regions of interest (ROI) and
edges given by the density of anatomical connections between regions [36]. The
level of activity of each ROI in a brain graph can be represented by a graph
signal $\mathbf{x}_{\text{i}}$, thus, successive applications of
$\mathbf{S}_{\text{i}}$ model a linear evolution of the brain activity
pattern. Under the hypothesis that we observe a sum of linear combinations
(filtered signals) resulting from diffusing originally sparse brain signals,
then blind demixing jointly estimates which regions of the brains were
originally active, the activity in these regions, and the diffusing
coefficients $\mathbf{h}_{\text{i}}$ of the linear combinations. In this
experiment, L and R vary, $\text{S}=1$, and the GSO $\mathbf{S}_{\text{i}}$ is
chosen to be the Laplacian matrix. The results show that demixing is indeed
feasible, although the performance decreases quickly as R increase. This is
not surprising since the brain graphs at hand exhibit strong similarities. On
the other hand, for a fixed value of R it is interesting that the performance
increases with L growing from 2 to 4, and it starts to drop afterwards. This
reveals a trade-off between the number of parameters to estimate (increasing
with L) and the amount of information of $\mathbf{x}_{\text{i}}$ in
$\mathbf{y}$ (which increases with the successive diffusions).

Effect of noise on recovery. In this experiment, we study the performance of
single-graph demixing (15) with $\text{L}=3$, varying R and S, and in the
presence of additive white Gaussian noise $\mathcal{N}(0,\sigma_{n}^{2})$ (see
Remark 1 for a discussion on handling noisy observations). The graph at hand
corresponds to the adjacency matrix of Zachary’s karate club network [37],
which contains $\text{N}=34$ nodes (the members of the club) and $78$
undirected edges (symbolizing friendships among members). The results reported
in Figure 2 reveal that, in all cases, the demixing error decreases as
$\sigma_{n}$ decreases and it increases as either R or S increase.
Furthermore, if the value $\text{R}+\text{S}$ is fixed, successful demixing is
easier for smaller values of R. Once again, this suggests that the sensitivity
of the demixing performance with respect to R is larger than that for S.

Figure 2:  Blind demixing in the karate network. Demixing Error (DE) as a
function of the noise level $\sigma_{n}$ for different number of filters (R)
and seeding nodes (S). The DE decreases as the sum $\text{R}+\text{S}$ also
decreases. When $\text{R}+\text{S}$ is fixed, smaller R results in higher
demixing accuracy. In the three cases with $\text{R}+\text{S}\geq 5$ where the
error does not converge to zero as the noise power goes to zero, demixing is
too hard given the size of the karate network ($\text{N}=34$) to achieve
successful recovery.

Predicting demixing performance for brain graphs. In Section IV, we presented
a theoretical result characterizing demixing performance as a function of a
number of parameters. When taking a close look at the proof (Section -A2), we
observe that the values of $\rho_{\mathbf{U}_{k}}^{(\text{S}_{k})}$,
$\rho_{\mathbf{U}_{j}}^{(\text{S}_{j})}$,
$\rho_{\mathbf{\Psi}_{k}}^{(\text{L}_{k})}$ and
$\rho_{\mathbf{\Psi}_{j}}^{(\text{L}_{j})}$ are intrinsic to demixing, account
for the interaction between mixing components, and have a larger impact on
recovery (see Section -A4 for details).

This motivates us to combine the values of
$\rho_{\mathbf{U}_{k}}^{(\text{S}_{k})}$,
$\rho_{\mathbf{U}_{j}}^{(\text{S}_{j})}$,
$\rho_{\mathbf{\Psi}_{k}}^{(\text{L}_{k})}$ and
$\rho_{\mathbf{\Psi}_{j}}^{(\text{L}_{j})}$ using the simple expression

$\bar{\rho}_{k,j}:=\rho_{\mathbf{U}_{k}}^{(\text{S}_{k})}\
\rho_{\mathbf{U}_{j}}^{(\text{S}_{j})}\
\rho_{\mathbf{\Psi}_{k}}^{(\text{L}_{k})}\
\rho_{\mathbf{\Psi}_{j}}^{(\text{L}_{j})},$ (21)

and, then, employ $\bar{\rho}_{k,j}$ as a predictor of the demixing
performance between a pair of graphs. Specifically, these quantities stem from
the center parameter of the concentration inequality (i.e. the bound in (37))
and the relation between $\kappa$ and $\rho$ in (11).

Hence, the results in this experiment try to asses the dependence of the
demixing performance with respect to $\bar{\rho}_{k,j}$. To that end, we
consider the six brain graphs used in the second experiment (cf. [36]), and
solve 15 (six choose two) different demixing problems, each of them involving
two out of the six graphs with $\text{L}=2$ and $\text{S}=1$. For each pair of
graphs, we use the Laplacian as GSO, compute $\bar{\rho}_{k,j}$ according to
(21) and quantify the rate of successful recovery.

(a)

(b)

Figure 3: (a) Recovery rate of the convex blind demixing minimization [cf.
(7)] as a function of the predictor (21). Each point corresponds to a pair of
brain graphs, where nodes of the same color share one fixed brain graph (of
the six in total). There are three outliers on the right side of the plot. The
straight line, result of a robust linear regression, shows the correlation
between the predictor and the recovery rate. (b) Explanation of the outliers
in (a) based on graph similarity for all pairs of brain graphs. The three
peaks correspond to the outliers in (a).

The obtained results, which are reported in Figure 3(a), confirm that pairs of
graphs with a larger value of $\bar{\rho}_{k,j}$ give rise to more challenging
demixing scenarios (i.e., smaller rate of recovery). In fact, if we fit the
obtained data using a linear model robust against outliers, we obtain the
dashed-straight line shown in the figure. This simple analysis reveals that
most of the points fall in the close vicinity of the linear fit, pointing out
that the predictor in (21) is a good _graph-dependent predictor_ for demixing
performance and, implicitly, validating the analysis and assumptions in
Theorem 1. Interestingly, the three pairs of graphs that are deemed as
outliers are the ones with most overlapping edges, as illustrated in Figure
3(b), which reports the number of overlapping edges for each pair of graphs
and highlights the three outliers identified in Figure 3(a). This could imply
that some of the assumptions adopted during the proof of the theorem could be
too pessimistic when the graphs at hand are very similar, opening the door to
future work focused on those types of graphs.

##  VII Conclusions

We formulated and analyzed the problem of blind demixing of filtered diffused
graph signals, an extension of blind demixing of time (or spatial) domain
signals to graphs. Envisioned application domains encompass a particular class
of bilinear inverse problems, where the observed graph signal is modeled as a
linear combination of diffusion processes driven by a number of sources. While
the graph signal observations are bilinear functions of the filter
coefficients and the sparse input signal, we rely on the frequency
interpretation of graph signals and filters to show that they are also
linearly related to the entries of a lifted rank-one, row-sparse matrix.
Accordingly, the blind demixing problem of graph signals was tackled via rank
and sparsity minimization subject to linear constraints, an inverse problem
amenable to convex relaxations offering provable recovery guarantees under
simplifying assumptions. The theoretical results reveal key features driving
the demixing performance including the number of sources and graph nodes, the
sparsity levels, the filter orders, as well as spectral properties of the
underlying graphs and their pairwise coherences. These spectral properties are
captured via conveniently defined concentration parameters. Numerical tests
validated the theoretical claims and demonstrated that the proposed approach
offers satisfactory recovery performance, even for settings well beyond the
scope of the analysis. Building on the insights gained here, it is of interest
to develop spectral graph theoretical tools to facilitate identifying the key
relationships between the (e.g., topological) properties of graphs and the
resulting recovery performance. Future research includes also estimating the
shift operator employing network topology inference [38, 39, 40].

###  -A Proof of Theorem 1

The proof is based on a set of conditions outlined in Lemma 1. In Section -A1
we show that, under these conditions, the nuclear norm minimization (8)
produces the same result as the original constrained search (6). For the
remaining of the proof, Sections -A2 and -A3 establish the conditions and,
finally, Section -A4 concludes the proof by implementing a relevant union
bound. We begin by introducing notation that will be used throughout the
proof.

Without loss of generality, it is assumed throughout that the true input
signals $\mathbf{x}_{\text{i}}^{(*)}$ and filter taps
$\mathbf{h}_{\text{i}}^{(*)}$ are real vectors with unit norm. Furthermore,
let $\mathbf{u}_{\text{i},n}^{\text{T}}$ and
$\boldsymbol{\psi}_{\text{i},n}^{\text{T}}$ denote the $n$th row of the GFT
matrices $\mathbf{U}_{\text{i}}$ and $\mathbf{\Psi}_{\text{i}}$, respectively
(cf. Section II-C). Then, consider the following linear operators
$\mathcal{M}_{\text{i}}:\mathbb{C}^{\text{N}\times\text{L}_{\text{i}}}\to\mathbb{C}^{\text{N}}$
and their adjoints
$\mathcal{M}^{*}_{\text{i}}:\mathbb{C}^{\text{N}}\to\mathbb{C}^{\text{N}\times\text{L}_{\text{i}}}$,

$\displaystyle\mathcal{M}_{\text{i}}(\mathbf{Z})$
$\displaystyle=\\{\mathbf{u}_{\text{i},n}^{\text{T}}\mathbf{Z}\boldsymbol{\psi}_{\text{i},n}\\}_{n=1}^{\text{N}},\quad\mathbf{Z}\in\mathbb{C}^{\text{N}\times\text{L}_{\text{i}}},$
$\displaystyle\mathcal{M}^{*}_{\text{i}}(\mathbf{z})$
$\displaystyle=\sum_{n=1}^{\text{N}}z_{n}\overline{\mathbf{u}_{\text{i},n}}\boldsymbol{\psi}_{\text{i},n}^{\text{H}},\quad\mathbf{z}\in\mathbb{C}^{\text{N}}.$

Let us also introduce a _common partition_
$\\{\Gamma_{\text{p}}\\}_{\text{p}=1}^{\text{Q}}$ of
$\\{1,2,\dots,\text{N}\\}$ with $|\Gamma_{\text{p}}|=\text{P}\ \forall\
\text{p}\in\\{1,2,\dots,\text{Q}\\}$ and $\text{N}=\text{Q}\text{P}$. See [27,
section II-C] and references therein for further discussion on the
partitioning. Associated with the partition, define the matrices

$\mathbf{D}_{\text{i},\text{p}}=\sum_{\ell\in\Gamma_{\text{p}}}\boldsymbol{\psi}_{\text{i},\ell}\boldsymbol{\psi}_{\text{i},\ell}^{\text{H}},\qquad\mathbf{F}_{\text{i},\text{p}}=\mathbf{D}_{\text{i},\text{p}}^{-1},$
(22)

the sampled linear operators

$\displaystyle\mathcal{M}_{\text{i},\text{p}}(\mathbf{Z})$
$\displaystyle=\\{\mathbf{u}_{\text{i},\ell}^{\text{T}}\mathbf{Z}\boldsymbol{\psi}_{\text{i},\ell}\\}_{\ell\in\Gamma_{\text{p}}},$
$\displaystyle\mathcal{M}^{*}_{\text{i},\text{p}}(\mathbf{z})$
$\displaystyle=\sum_{\ell\in\Gamma_{\text{p}}}z_{\ell}\overline{\mathbf{u}_{\text{i},\ell}}\boldsymbol{\psi}_{\text{i},\ell}^{\text{H}},$

and the following scalar quantities

$\mu_{\text{h}}^{2}:=\max\Bigg{\\{}\frac{\text{Q}^{2}}{\text{N}}\max_{\begin{subarray}{c}\ell\in\Gamma_{\text{p}},\\\
1\leq\text{p}\leq\text{Q},\\\
1\leq\text{i}\leq\text{R}\end{subarray}}\frac{|\boldsymbol{\psi}_{\text{i},n}^{\text{H}}\mathbf{F}_{\text{i},\text{p}}\mathbf{h}_{\text{i}}|^{2}}{\|\mathbf{h}_{\text{i}}\|^{2}},\text{N}\max_{\begin{subarray}{c}1\leq
n\leq\text{N},\\\
1\leq\text{i}\leq\text{R}\end{subarray}}\frac{|\boldsymbol{\psi}_{\text{i},n}^{\text{H}}\mathbf{h}_{\text{i}}|^{2}}{\|\mathbf{h}_{\text{i}}\|^{2}}\Bigg{\\}}$
(23)

capturing the interaction between $\mathbf{h}_{\text{i}}$ and
$\mathbf{\Psi}_{\text{i}}$, and

$\mu_{\text{max}}^{2}:=\max_{1\leq
n\leq\text{N},1\leq\text{i}\leq\text{R}}\frac{\text{N}}{\text{L}_{\text{i}}}\|\boldsymbol{\psi}_{\text{i},n}\|^{2},$
(24)

providing a notion of (in)coherence of the matrices
$\mathbf{\Psi}_{\text{i}}$.

Additionally, we will handle projections on spaces such as those constrained
to the signals supports $\Omega_{\text{i}}$,

$\displaystyle{\mathcal{M}_{\text{i}}}_{\Omega_{\text{i}}}(\mathbf{Z})$
$\displaystyle=\\{{\mathbf{u}_{\text{i}}^{\text{T}}}_{\Omega_{\text{i}}}\mathbf{Z}\boldsymbol{\psi}_{\text{i}}\\}_{\text{i}=1}^{\text{N}},$
$\displaystyle{\mathcal{M}^{*}_{\text{i}}}_{\Omega_{\text{i}}}(\mathbf{z})$
$\displaystyle=\sum_{\text{i}=1}^{\text{R}}\mathbf{z}{\overline{\mathbf{u}_{\text{i}}}}_{\Omega_{\text{i}}}\boldsymbol{\psi}_{\text{i}}^{\text{H}}.$

Naturally, projected _and_ sampled operators will be denoted by

$\displaystyle{\mathcal{M}_{\text{i},\text{p}}}_{\Omega_{\text{i}}}(\mathbf{Z})$
$\displaystyle=\\{{\mathbf{u}_{\text{i},\ell}^{\text{T}}}_{\Omega_{\text{i}}}\mathbf{Z}\boldsymbol{\psi}_{\text{i},\ell}\\}_{\ell\in\Gamma_{\text{p}}},$
$\displaystyle{\mathcal{M}^{*}_{\text{i},\text{p}}}_{\Omega_{\text{i}}}(\mathbf{z})$
$\displaystyle=\sum_{\ell\in\Gamma_{\text{p}}}\mathbf{z}_{\ell}{\overline{\mathbf{u}_{\text{i},\ell}}}_{\Omega_{\text{i}}}.\boldsymbol{\psi}_{\text{i},\ell}^{\text{H}}.$

Although making slight abuse overloading the symbols, the subindex notation
denoting projection onto a space will be used indistinctly for linear
operators, matrices, and vectors.

####  -A1 Sufficient conditions

In [27, Lemma 5.1] it is shown that for any real
$\\{\mathbf{H}_{\text{i}}\\}_{\text{i}=1}^{\text{R}}$ satisfying that
$\sum_{\text{i}=1}^{\text{R}}\mathcal{M}_{\text{i}}(\mathbf{H}_{\text{i}})=0$
with at least one $\mathbf{H}_{\text{i}}$ nonzero, if

$\sum_{\text{i}=1}^{\text{R}}\left\langle\mathbf{H}_{\text{i}},\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}\right\rangle+\|\mathbf{H}_{\text{i},\Omega_{\text{i}}^{\perp}}\|_{*}>0$
(25)

then the solution to the convex program,

$\operatorname*{minimize}_{\\{\mathbf{x}_{\text{i}},\mathbf{h}_{\text{i}}\\}_{\text{i}=1}^{\text{R}}}\sum_{\text{i}=1}^{\text{R}}\|\mathbf{x}_{\text{i}}\mathbf{h}_{\text{i}}^{\text{T}}\|_{*}\quad\text{s.
to}\quad\widehat{\boldsymbol{\nu}}=\sum_{\text{i}=1}^{\text{R}}\mathcal{M}_{\text{i}}\left(\mathbf{x}_{\text{i}}\mathbf{h}_{\text{i}}^{\text{T}}\right),$
(26)

is the set formed by the true $\mathbf{x}_{\text{i}}^{(*)}$ and
$\mathbf{h}_{\text{i}}^{(*)}$. Nevertheless, this lemma does not provide with
conditions on the linear operators, which capture the blind demixing model. To
this end, inspired by [27, Lemma 5.2], we state the following lemma.

###### Lemma 1

Assume that, for all $i$,

$\|{\mathcal{M}_{\text{i}}^{*}}_{\Omega_{\text{i}}}{\mathcal{M}_{\text{i}}}_{\Omega_{\text{i}}}-\mathbf{I}_{\Omega_{\text{i}}}\|\leq\frac{1}{4},$
(27) $\mu:=\max_{j\neq
k}\|{\mathcal{M}_{j}^{*}}_{\Omega_{j}}{\mathcal{M}_{k}}_{\Omega_{k}}\|\leq\frac{1}{4\text{R}},$
(28) $\|\mathcal{M}_{\text{i}}\|\leq\gamma,$ (29)

for some $\gamma>0$, and that there exists a
$\boldsymbol{\lambda}\in\mathbb{C}^{\text{N}}$ (denominated _dual
certificate_) such that, for all $i$,

${\|\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}-{\mathcal{M}_{\text{i}}^{*}}_{\Omega_{\text{i}}}(\boldsymbol{\lambda})\|_{\text{F}}\leq\frac{1}{5\text{R}\gamma}},$
(30)
$\|{\mathcal{M}_{\text{i}}^{*}}_{\Omega_{\text{i}}^{\perp}}(\boldsymbol{\lambda})\|\leq\frac{1}{2}.$
(31)

Then, the solution to (26) is given by the true signals
$\mathbf{x}_{\text{i}}^{(*)}$ and filter taps $\mathbf{h}_{\text{i}}^{(*)}$.

_Proof :_ We begin by defining $I_{1}$ as

$I_{1}:=\sum_{\text{i}=1}^{\text{R}}\left\langle\mathbf{H}_{\text{i}},\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}-\mathcal{M}_{\text{i}}^{*}(\boldsymbol{\lambda})\right\rangle+\|\mathbf{H}_{\text{i},\Omega_{\text{i}}^{\perp}}\|_{*}.$

Our goal is to show that for any nonzero
$\\{\mathbf{H}_{\text{i}}\\}_{\text{i}=1}^{\text{R}}$ such that
$\sum_{\text{i}=1}^{\text{R}}\mathcal{M}_{\text{i}}(\mathbf{H}_{\text{i}})=0$,
it holds that $I_{1}>0$. We begin by decomposing the inner product into
$\Omega_{\text{i}}$ and $\Omega_{\text{i}}^{\perp}$ to obtain

$\begin{split}I_{1}=\sum_{\text{i}=1}^{\text{R}}&\
\Big{(}\langle\mathbf{H}_{\text{i},\Omega_{\text{i}}},\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}-{\mathcal{M}_{\text{i}}^{*}}_{\Omega_{\text{i}}}(\boldsymbol{\lambda})\rangle\\\
&-\langle\mathbf{H}_{\text{i},\Omega_{\text{i}}^{\perp}},{\mathcal{M}_{\text{i}}^{*}}_{\Omega_{\text{i}}^{\perp}}(\boldsymbol{\lambda})\rangle+\|\mathbf{H}_{\text{i},\Omega_{\text{i}}^{\perp}}\|_{*}\Big{)}.\end{split}$

Let us additionally define $I_{2}$ as follows

$\begin{split}I_{2}:=\sum_{\text{i}=1}^{\text{R}}\Big{(}-&\|\mathbf{H}_{\text{i},\Omega_{\text{i}}}\|_{\text{F}}\
\|\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}-{\mathcal{M}_{\text{i}}^{*}}_{\Omega_{\text{i}}}(\boldsymbol{\lambda})\|_{\text{F}}\\\
+&\|\mathbf{H}_{\text{i},\Omega_{\text{i}}^{\perp}}\|_{*}(1-\|{\mathcal{M}_{\text{i}}^{*}}_{\Omega_{\text{i}}^{\perp}}(\boldsymbol{\lambda})\|)\Big{)}.\end{split}$

We are now going to rely on the following two matrix inequalities. First,

$\langle\mathbf{U},\mathbf{V}\rangle\geq-|\langle\mathbf{U},\mathbf{V}\rangle|\geq-\|\mathbf{U}\|_{\text{F}}\|\mathbf{V}\|_{\text{F}},$
(32)

which follows from the Cauchy-Schwarz inequality. Second,

$-\langle\mathbf{U},\mathbf{V}\rangle\geq-|\langle\mathbf{U},\mathbf{V}\rangle|\geq-\|\mathbf{U}\|_{*}\|\mathbf{V}\|,$
(33)

which can be shown by applying a SVD of $\mathbf{U}$ followed by the fact that
the matrix norm $\|\cdot\|$ is invariant to orthogonal transformations.
Leveraging (32) and (33), it follows that $I_{1}\geq I_{2}$. Thus, we now seek
to show that $I_{2}>0$.

Making use of $\|\cdot\|_{*}\geq\|\cdot\|_{\text{F}}$ we have that

$\begin{split}I_{2}\geq\sum_{\text{i}=1}^{\text{R}}\Big{(}&-\|\mathbf{H}_{\text{i},\Omega_{\text{i}}}\|_{\text{F}}\|\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}-{\mathcal{M}_{\text{i}}^{*}}_{\Omega_{\text{i}}}(\boldsymbol{\lambda}))\|_{\text{F}}\\\
&+\|\mathbf{H}_{\text{i},\Omega_{\text{i}}^{\perp}}\|_{\text{F}}(1-\|{\mathcal{M}_{\text{i}}^{*}}_{\Omega_{\text{i}}^{\perp}}(\boldsymbol{\lambda}))\|)\Big{)}.\end{split}$

Next, from (30) and (31) it follows that

$I_{2}\geq-\frac{1}{5\text{R}\gamma}\sum_{\text{i}=1}^{\text{R}}\|\mathbf{H}_{\text{i},\Omega_{\text{i}}}\|_{\text{F}}+\frac{1}{2}\sum_{\text{i}=1}^{\text{R}}\|\mathbf{H}_{\text{i},\Omega_{\text{i}}^{\perp}}\|_{\text{F}}.$

Moreover, from (27)-(29), it can be shown that

$\frac{1}{2\text{R}}\sum_{\text{i}=1}^{\text{R}}\|\mathbf{H}_{\text{i},\Omega_{\text{i}}}\|_{\text{F}}\leq\gamma\sum_{\text{i}=1}^{\text{R}}\|\mathbf{H}_{\text{i},\Omega_{\text{i}}^{\perp}}\|_{\text{F}}.$
(34)

The derivation of (34) follows the same steps as those in the proof of [27,
Lemma 5.2]. Hence, we avoid it here to reduce redundancy. Leveraging (34) we
obtain that

$I_{2}\geq-\frac{2}{5}\sum_{\text{i}=1}^{\text{R}}\|\mathbf{H}_{\text{i},\Omega_{\text{i}}^{\perp}}\|_{\text{F}}+\frac{1}{2}\sum_{\text{i}=1}^{\text{R}}\|\mathbf{H}_{\text{i},\Omega_{\text{i}}^{\perp}}\|_{\text{F}}=\frac{1}{10}\sum_{\text{i}=1}^{\text{R}}\|\mathbf{H}_{\text{i},\Omega_{\text{i}}^{\perp}}\|_{\text{F}},$

which is indeed larger than zero since at least one of the terms
$\|\mathbf{H}_{\text{i},\Omega_{\text{i}}^{\perp}}\|_{\text{F}}$ in the sum
must be nonzero owing to the fact that via (34)
$\mathbf{H}_{\text{i},\Omega_{\text{i}}^{\perp}}\\!=\\!0\ \forall\
1\\!\leq\\!\text{i}\\!\leq\\!\text{R}$ would imply that all
$\mathbf{H}_{\text{i}}\\!=\\!0$. Finally, knowing that $I_{2}>0$ and therefore
$I_{1}>0$ too, from the definition of an adjoint operator it follows that

$\displaystyle I_{1}$
$\displaystyle=\sum_{\text{i}=1}^{\text{R}}\left\langle\mathbf{H}_{\text{i}},\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}\right\rangle-\sum_{\text{i}=1}^{\text{R}}\left\langle\mathcal{M}_{\text{i}}(\mathbf{H}_{\text{i}}),\boldsymbol{\lambda}\right\rangle+\|\mathbf{H}_{\text{i},\Omega_{\text{i}}^{\perp}}\|_{*}$
$\displaystyle=\sum_{\text{i}=1}^{\text{R}}\left\langle\mathbf{H}_{\text{i}},\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}\right\rangle+\|\mathbf{H}_{\text{i},\Omega_{\text{i}}^{\perp}}\|_{*}>0.$

Having shown (25), the statement of the lemma holds. $\hfill\blacksquare$

At this point we are left occupied with showing (27)-(31). Condition (27)
follows from [21, Lemma 2] and (29) from [21, Lemma 3]. The probabilities with
which these two conditions hold are revisited in Section -A4. In the following
two sections we develop proofs for the remaining three conditions: the
coherence between pairs of linear operators (28) is shown in Section -A2 and
the two dual certificate conditions (30) and (31) are shown in Section -A3.
Finally, Section -A4 synthesizes the probabilities to fulfill each of the five
conditions.

####  -A2 Coherence between pairs of linear operators

A measure of the _coherence_ between two linear operators $\mathcal{M}_{j}$
and $\mathcal{M}_{k}$ is $\mu$ as defined in (28). The following lemma
provides a lower bound on the probability that (28) holds.

###### Lemma 2

Condition (28) holds with probability at least

$\displaystyle 1-2\text{N}$
$\displaystyle\exp\Bigg{(}-\frac{1}{32\text{R}^{2}}\Big{[}\frac{16\text{N}^{2}}{9Q^{2}}\rho_{\mathbf{\Psi}_{k}}^{(L_{k})}\rho_{\mathbf{\Psi}_{j}}^{(L_{j})}\max\left\\{\rho_{\mathbf{U}_{k}}^{(S_{k})},\rho_{\mathbf{U}_{j}}^{(S_{j})}\right\\}$
$\displaystyle+\frac{\text{N}}{9Q\text{R}}\kappa_{\mathbf{U}_{k},\mathbf{U}_{j}}^{\left(\text{S}_{k},\text{S}_{j}\right)}\kappa_{\mathbf{\Psi}_{k},\mathbf{\Psi}_{j}}^{\left(\text{L}_{k},\text{L}_{j}\right)}\
\Big{]}^{-1}\Bigg{)}.$ (35)

_Proof :_ Consider the following expression, more general than $\mu$ [cf.
(28)]

$\|{\mathcal{M}_{j,\text{p}}^{*}}_{\Omega_{j}}{\mathcal{M}_{k,\text{p}}}_{\Omega_{k}}(\mathbf{Z}\mathbf{F}_{k,\text{p}})\|.$
(36)

To see that (36) indeed generalizes $\mu$, consider a partition given by a
single set $\Gamma_{1}$. In this case, from
$\mathbf{\Psi}_{\text{i}}^{\text{H}}\mathbf{\Psi}_{\text{i}}=\mathbf{I}_{\text{L}_{\text{i}}}$
222Under the assumption that the GSOs have distinct eigenvalues [cf. (AS1) in
Section IV-A], it can be assumed without loss of generality that
$\mathbf{\Psi}_{\text{i}}^{\text{H}}\mathbf{\Psi}_{\text{i}}=\mathbf{I}_{\text{L}_{\text{i}}}$.
First of all, since each $\mathbf{S}_{\text{i}}$ has distinct eigenvalues,
then the matrices $\mathbf{\Psi}_{\text{i}}$ are full-rank. Second, consider
the SVD of
$\mathbf{\Psi}_{\text{i}}=\mathbf{L}_{\text{i}}\mathbf{\Sigma}_{\text{i}}\mathbf{R}_{\text{i}}^{\text{H}}$;
accordingly the frequency response of the filter
$\mathbf{\Psi}_{\text{i}}\mathbf{h}_{\text{i}}=\mathbf{L}_{\text{i}}\mathbf{\Sigma}_{\text{i}}\mathbf{R}_{\text{i}}^{\text{H}}\mathbf{h}_{\text{i}}$,
and define
${\mathbf{h}^{\prime}}_{\text{i}}:=\mathbf{\Sigma}_{\text{i}}\mathbf{R}_{\text{i}}\mathbf{h}_{\text{i}}$.
Since $\mathbf{\Psi}_{\text{i}}$ is full-rank, $\mathbf{h}_{\text{i}}$ can be
recovered from ${\mathbf{h}^{\prime}}_{\text{i}}$, and finally since the
matrix of left-singular vectors matrix holds that
$\mathbf{L}_{\text{i}}^{\text{H}}\mathbf{L}_{\text{i}}=\mathbf{I}_{\text{L}_{\text{i}}}$,
we can assume that
$\mathbf{\Psi}_{\text{i}}^{\text{H}}\mathbf{\Psi}_{\text{i}}=\mathbf{I}_{\text{L}_{\text{i}}}$.
it follows that
$\mathbf{D}_{\text{i},1}=\mathbf{F}_{\text{i},1}=\mathbf{I}_{\text{L}_{\text{i}}}$.
Also, since $\Gamma_{1}=\\{1,2,\dots,\text{N}\\}$ then
$\mathcal{M}_{\text{i},1}=\mathcal{M}_{\text{i}}$ in the single-set partition
case.

We shall proceed setting $j=2$ and $k=1$ to simplify the notation and work
with the expression divided into two parts,

${\mathcal{M}_{2,\text{p}}^{*}}_{\Omega_{2}}(\mathbf{z})=\sum_{\ell\in\Gamma_{\text{p}}}z_{\ell}{\overline{\mathbf{u}_{2,\ell}}}_{\Omega_{2}}\boldsymbol{\psi}_{2,\ell}^{\text{H}},$
${\mathcal{M}_{1,\text{p}}}_{\Omega_{1}}(\mathbf{Z}\mathbf{F}_{1,\text{p}})=\left\\{{\mathbf{u}_{1,\ell}^{\text{T}}}_{\Omega_{1}}\mathbf{Z}\mathbf{F}_{1,\text{p}}\boldsymbol{\psi}_{1,\ell}\right\\}_{\ell\in\Gamma_{\text{p}}}.$

Concatenating the two operators above,

$\sum_{\ell\in\Gamma_{\text{p}}}({\mathbf{u}_{1,\ell}^{\text{T}}}_{\Omega_{1}}\mathbf{Z}\mathbf{F}_{1,\text{p}}\boldsymbol{\psi}_{1,\ell})\overline{\mathbf{u}_{2,\ell}}_{\Omega_{2}}\boldsymbol{\psi}_{2,\ell}^{\text{H}}$

and examining each summand we have that

$\left({\mathbf{u}_{1,\ell}^{\text{T}}}_{\Omega_{1}}\mathbf{Z}\mathbf{F}_{1,\text{p}}\boldsymbol{\psi}_{1,\ell}\right)\overline{\mathbf{u}_{2,\ell}}_{\Omega_{2}}\boldsymbol{\psi}_{2,\ell}^{\text{H}}=\overline{\mathbf{u}_{2,\ell}}_{\Omega_{2}}{\mathbf{u}_{1,\ell}^{\text{T}}}_{\Omega_{1}}\mathbf{Z}\mathbf{F}_{1,\text{p}}\boldsymbol{\psi}_{1,\ell}\boldsymbol{\psi}_{2,\ell}^{\text{H}}.$

We seek to bound (36) by applying the non-commutative Bernstein inequality for
matrices [21, Theorem 2] to a sum of matrices given by
$\overline{\mathbf{u}_{2,\ell}}_{\Omega_{2}}{\mathbf{u}_{1,\ell}^{\text{T}}}_{\Omega_{1}}\mathbf{Z}\
\mathbf{F}_{1,\text{p}}\boldsymbol{\psi}_{1,\ell}\boldsymbol{\psi}_{2,\ell}^{\text{H}}$.
Following the notation in [21, Theorem 2], the operators to bound in the
concentration inequality are a sequence of

$\mathcal{Z}_{\ell}(\mathbf{Z})=\overline{\mathbf{u}_{2,\ell}}_{\Omega_{2}}{\mathbf{u}_{1,\ell}^{\text{T}}}_{\Omega_{1}}\mathbf{Z}\
\mathbf{F}_{1,\text{p}}\boldsymbol{\psi}_{1,\ell}\boldsymbol{\psi}_{2,\ell}^{\text{H}}.$

Note that, given its definition in (22), $\mathbf{F}_{1,\text{p}}$ is
Hermitian. Then, the adjoint operators become

$\mathcal{Z}_{\ell}^{*}(\mathbf{Z})=\overline{\mathbf{u}_{1,\ell}}_{\Omega_{1}}{\mathbf{u}_{2,\ell}^{\text{T}}}_{\Omega_{2}}\mathbf{Z}\
\boldsymbol{\psi}_{2,\ell}\boldsymbol{\psi}_{1,\ell}^{\text{H}}\mathbf{F}_{1,\text{p}}.$

Recalling that $\text{S}_{\text{i}}=|\Omega_{\text{i}}|$ and leveraging the
definition of $\kappa$ in (10), we have that

$\|\mathcal{Z}_{\ell}\|\leq\|\mathbf{F}_{1,\text{p}}\|\kappa_{\mathbf{U}_{1},\mathbf{U}_{2}}^{\left(\text{S}_{1},\text{S}_{2}\right)}\kappa_{\mathbf{\Psi}_{1},\mathbf{\Psi}_{2}}^{\left(\text{L}_{1},\text{L}_{2}\right)}\leq\frac{4\text{N}}{3\text{Q}}\kappa_{\mathbf{U}_{1},\mathbf{U}_{2}}^{\left(\text{S}_{1},\text{S}_{2}\right)}\kappa_{\mathbf{\Psi}_{1},\mathbf{\Psi}_{2}}^{\left(\text{L}_{1},\text{L}_{2}\right)}$
(37)

where $\|\mathbf{F}_{\text{i},\text{p}}\|\\!\leq\\!4\text{N}/(3\text{Q})\
\forall\
(i,\text{p})\\!\in\\!\\{1,2,\dots,\text{R}\\}\times\\{1,2,\dots,\text{P}\\}$
from [27, (12)]. This is a direct consequence of the definitions of the
partition $\Gamma_{\text{p}}$ and the matrices
$\mathbf{D}_{\text{i},\text{p}}$ lying close to
$\mathbf{I}_{\text{L}_{\text{i}}}$, refer to [27, section II-C] for further
discussion on this.

In order to implement the concentration inequality, we compute the expected
value of the concatenation operators
$\mathcal{Z}_{\ell}\mathcal{Z}_{\ell}^{*}(\mathbf{Z})$ and
$\mathcal{Z}_{\ell}^{*}\mathcal{Z}_{\ell}(\mathbf{Z})$ to obtain

$\displaystyle\mathbb{E}\left(\mathcal{Z}_{\ell}\mathcal{Z}_{\ell}^{*}(\mathbf{Z})\right)=$
$\displaystyle\left(\boldsymbol{\psi}_{1,\ell}^{\text{H}}\mathbf{F}_{1,\text{p}}\mathbf{F}_{1,\text{p}}\boldsymbol{\psi}_{1,\ell}\right)\mathbb{E}\left({\mathbf{u}_{1,\ell}^{\text{T}}}_{\Omega_{1}}\overline{\mathbf{u}_{1,\ell}}_{\Omega_{1}}\right)$
$\displaystyle\mathbb{E}\left(\overline{\mathbf{u}_{2,\ell}}_{\Omega_{2}}{\mathbf{u}_{2,\ell}^{\text{T}}}_{\Omega_{2}}\right)\mathbf{Z}\
\boldsymbol{\psi}_{2,\ell}\boldsymbol{\psi}_{2,\ell}^{\text{H}},$
$\displaystyle\mathbb{E}\left(\mathcal{Z}_{\ell}^{*}\mathcal{Z}_{\ell}(\mathbf{Z})\right)=$
$\displaystyle\left(\boldsymbol{\psi}_{2,\ell}^{\text{H}}\boldsymbol{\psi}_{2,\ell}\right)\mathbb{E}\left({\mathbf{u}_{2,\ell}^{\text{T}}}_{\Omega_{2}}\overline{\mathbf{u}_{2,\ell}}_{\Omega_{2}}\right)$
$\displaystyle\mathbb{E}\left(\overline{\mathbf{u}_{1,\ell}}_{\Omega_{1}}{\mathbf{u}_{1,\ell}^{\text{T}}}_{\Omega_{1}}\right)\mathbf{Z}\
\mathbf{F}_{1,\text{p}}\boldsymbol{\psi}_{1,\ell}\boldsymbol{\psi}_{1,\ell}^{\text{H}}\mathbf{F}_{1,\text{p}}.$

Due to the normalization
$\mathbf{U}_{\text{i}}^{\text{H}}\mathbf{U}_{\text{i}}=\text{N}\mathbf{I}_{\text{N}}$,
it follows that
$\mathbb{E}(\overline{\mathbf{u}_{\text{i},\ell}}_{\Omega_{\text{i}}}{\mathbf{u}_{\text{i},\ell}^{\text{T}}}_{\Omega_{\text{i}}})=\mathbf{I}_{\text{N}_{\Omega_{\text{i}}}}$.
Moreover, from the definition in (9), we have that the norms of
$\mathbb{E}({\mathbf{u}_{\text{i},\ell}^{\text{T}}}_{\Omega_{\text{i}}}\overline{\mathbf{u}_{\text{i},\ell}}_{\Omega_{\text{i}}})$
are bounded by $\rho_{\mathbf{U}_{\text{i}}}^{(\text{S}_{\text{i}})}$.
Therefore,

$\displaystyle\left\|\sum_{\ell\in\Gamma_{\text{p}}}\mathbb{E}\left(\mathcal{Z}_{\ell}\mathcal{Z}_{\ell}^{*}\right)\right\|$
$\displaystyle\leq\rho_{\mathbf{U}_{1}}^{(\text{S}_{1})}\sum_{\ell\in\Gamma_{\text{p}}}\|\boldsymbol{\psi}_{1,\ell}^{\text{H}}\mathbf{F}_{1,\text{p}}\mathbf{F}_{1,\text{p}}\boldsymbol{\psi}_{1,\ell}\|\
\|\boldsymbol{\psi}_{2,\ell}\boldsymbol{\psi}_{2,\ell}^{\text{H}}\|$
$\displaystyle\leq\frac{16\text{N}^{2}}{9\text{Q}^{2}}\rho_{\mathbf{U}_{1}}^{(\text{S}_{1})}\rho_{\mathbf{\Psi}_{1}}^{(\text{L}_{1})}\rho_{\mathbf{\Psi}_{2}}^{(\text{L}_{2})},$
$\displaystyle\left\|\sum_{\ell\in\Gamma_{\text{p}}}\mathbb{E}\left(\mathcal{Z}_{\ell}^{*}\mathcal{Z}_{\ell}\right)\right\|$
$\displaystyle\leq\rho_{\mathbf{U}_{2}}^{(\text{S}_{2})}\sum_{\ell\in\Gamma_{\text{p}}}\|\boldsymbol{\psi}_{2,\ell}^{\text{H}}\boldsymbol{\psi}_{2,\ell}\|\
\|\mathbf{F}_{1,\text{p}}\boldsymbol{\psi}_{1,\ell}\boldsymbol{\psi}_{1,\ell}^{\text{H}}\mathbf{F}_{1,\text{p}}\|$
$\displaystyle\leq\frac{16\text{N}^{2}}{9\text{Q}^{2}}\rho_{\mathbf{U}_{2}}^{(\text{S}_{2})}\rho_{\mathbf{\Psi}_{1}}^{(\text{L}_{1})}\rho_{\mathbf{\Psi}_{2}}^{(\text{L}_{2})}$

where we have used, once again, that $\|\mathbf{F}_{\text{i},\text{p}}\|\leq
4\text{N}/(3\text{Q})$. Finally, the variance parameter in the Bernstein
inequality is given by

$\sigma^{2}=\frac{16\text{N}^{2}}{9\text{Q}^{2}}\rho_{\mathbf{\Psi}_{1}}^{(\text{L}_{1})}\rho_{\mathbf{\Psi}_{2}}^{(\text{L}_{2})}\max\left\\{\rho_{\mathbf{U}_{1}}^{(\text{S}_{1})},\rho_{\mathbf{U}_{2}}^{(\text{S}_{2})}\right\\}.$
(38)

Employing the concentration inequality as stated in [21, Theorem 2] with the
parameters obtained in (37) and (38) we have that

$\displaystyle\mathbf{P}\\!\left(\mu\\!\geq\\!t\right)\\!\leq\\!2\text{N}\exp\Bigg{(}\\!\\!$
$\displaystyle-\frac{t^{2}}{2}\Big{[}\frac{16\text{N}^{2}}{9Q^{2}}\rho_{\mathbf{\Psi}_{k}}^{(\text{L}_{k})}\rho_{\mathbf{\Psi}_{j}}^{(\text{L}_{j})}\max\left\\{\rho_{\mathbf{U}_{k}}^{(\text{S}_{k})},\rho_{\mathbf{U}_{j}}^{(\text{S}_{j})}\right\\}$
$\displaystyle+\frac{t}{3}\frac{4\text{N}}{3\text{Q}}\kappa_{\mathbf{U}_{k},\mathbf{U}_{j}}^{\left(\text{S}_{k},\text{S}_{j}\right)}\kappa_{\mathbf{\Psi}_{k},\mathbf{\Psi}_{j}}^{\left(\text{L}_{k},\text{L}_{j}\right)}\Big{]}^{-1}\Bigg{)}.$

The final result of the lemma is obtained by setting $t=1/(4\text{R})$ and
switching from a lower to an upper bound. $\hfill\blacksquare$

####  -A3 Dual certificate

The approximate dual certificate
$\\{\mathbf{Y}_{\text{i}}:=\mathcal{M}^{*}_{\text{i}}(\boldsymbol{\lambda})\\}_{\text{i}=1}^{\text{R}}$
is constructed via the golfing scheme, producing a sequence of random matrices
[27]. The purpose behind the method is that the sequence
$\\{\mathbf{Y}_{\text{i},\text{p}}\\}_{\text{p}=1}^{\text{P}}$ approaches
$\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}$ on
$\Omega_{\text{i}}$ exponentially fast while keeping
$\mathbf{Y}_{\text{i},\text{p}}$ small on $\Omega^{\perp}_{\text{i}}$. The
method is now outlined.

  1. 1.

Initialization:

$\displaystyle\boldsymbol{\lambda}_{0}$
$\displaystyle:=\sum_{\text{i}=1}^{\text{R}}\mathcal{M}_{\text{i},\text{p}}(\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}\mathbf{F}_{\text{i},\text{p}})\
\in\mathbb{C}^{\text{N}},$ (39) $\displaystyle\mathbf{Y}_{\text{i},0}$
$\displaystyle:=\mathbf{0}_{\text{N}\times\text{L}_{\text{i}}},\quad\forall\
\text{i}\in\\{1,2,\dots,\text{R}\\}.$ (40)

  2. 2.

Recursion:

$\displaystyle\boldsymbol{\lambda}_{\text{p}-1}$
$\displaystyle:=\sum_{\text{i}=1}^{\text{R}}\mathcal{M}_{\text{i},\text{p}}\left((\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}-{\mathbf{Y}_{\text{i},\text{p}-1}}_{\Omega_{\text{i}}})\mathbf{F}_{\text{i},\text{p}}\right),$
(41) $\displaystyle\mathbf{Y}_{\text{i},\text{p}}$
$\displaystyle:=\mathbf{Y}_{\text{i},\text{p}-1}+\mathcal{M}^{*}_{\text{i},\text{p}}(\boldsymbol{\lambda}_{\text{p}-1}),\quad\forall\
\text{i}\in\\{1,2,\dots,\text{R}\\}.$ (42)

Let $\mathbf{Y}_{\text{i}}:=\mathbf{Y}_{\text{i},\text{P}}$ for all i and
$\boldsymbol{\lambda}:=\boldsymbol{\lambda}_{\text{P}}$. Furthermore, define
the difference

$\mathbf{W}_{\text{i},\text{p}}=\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}-{\mathbf{Y}_{\text{i},\text{p}}}_{\Omega_{\text{i}}},\quad\mathbf{W}_{\text{i},0}=\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}$
(43)

which leads to rewriting

$\boldsymbol{\lambda}_{\text{p}-1}=\sum_{\text{i}=1}^{\text{R}}\mathcal{M}_{\text{i},\text{p}}(\mathbf{W}_{\text{i},\text{p}-1}\mathbf{F}_{\text{i},\text{p}})$
(44)

and the recursion

$\mathbf{W}_{\text{i},\text{p}}=\mathbf{W}_{\text{i},\text{p}-1}-\sum_{\text{j}=1}^{\text{R}}{\mathcal{M}^{*}_{\text{j},\text{p}}}_{\Omega_{\text{i}}}{\mathcal{M}_{\text{j},\text{p}}}_{\Omega_{\text{i}}}(\mathbf{W}_{\text{j},\text{p}-1}\mathbf{F}_{\text{j},\text{p}}).$
(45)

Note that by definition the residuals hold that
${\mathbf{W}_{\text{j},\text{p}-1}}_{\Omega^{\perp}_{\text{i}}}=0$.

The dual certificate comprises two conditions: the exponential decay of
$\|\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}-{\mathbf{Y}_{\text{i}}}_{\Omega_{\text{i}}}\|_{\text{F}}$,
given in Lemma 3, and
$\|{\mathbf{Y}_{\text{i}}}_{\Omega^{\perp}_{\text{i}}}\|\leq 1/2$ in Lemma 4.

###### Lemma 3

Under the condition that $\text{P}\\!\geq\\!\log_{2}(5\text{R}\gamma)$, then
$\|\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}-{\mathbf{Y}_{\text{i}}}_{\Omega_{\text{i}}}\|_{\text{F}}\leq(5\text{R}\gamma)^{-1}$
with probability at least $1-\text{N}^{-\alpha+1}$.

_Proof :_ The following development follows the steps from [27, Lemma 5.12].
Expanding the sum in (45),

$\begin{split}\mathbf{W}_{\text{i},\text{p}}=\mathbf{W}_{\text{i},\text{p}-1}-&\mathcal{M}^{*}_{\text{i},\text{p}}\mathcal{M}_{\text{i},\text{p}}(\mathbf{W}_{\text{i},\text{p}-1}\mathbf{F}_{\text{i},\text{p}})\\\
-\sum_{\text{j}=1,\
\text{j}\neq\text{i}}^{\text{R}}&{\mathcal{M}^{*}_{\text{j},\text{p}}}_{\Omega_{\text{i}}}{\mathcal{M}_{\text{j},\text{p}}}_{\Omega_{\text{i}}}(\mathbf{W}_{\text{j},\text{p}-1}\mathbf{F}_{\text{j},\text{p}}).\end{split}$
(46)

Applying [21, Lemma 2] ($\delta=1/4$, $\text{P}=1$ therein), $\mu\leq
1/(4\text{R})$ (cf. Lemma 2), and noting that
$\mathbf{W}_{\text{j},\text{p}-1}={\mathbf{W}_{\text{j},\text{p}-1}}_{\Omega_{\text{j}}}$:

$\|\mathbf{W}_{\text{i},\text{p}}\|_{\text{F}}\leq\frac{1}{4}\|\mathbf{W}_{\text{i},\text{p}-1}\|_{\text{F}}+\frac{1}{4\text{R}}\sum_{\text{j}\neq\text{i}}\|\mathbf{W}_{\text{j},\text{p}-1}\|_{\text{F}},\quad
1\leq\text{i}\leq\text{R}$ (47)

and simplifying with the max

$\|\mathbf{W}_{\text{i},\text{p}}\|_{\text{F}}\leq\frac{1}{2}\max_{1\leq\text{i}\leq\text{R}}\|\mathbf{W}_{\text{i},\text{p}-1}\|_{\text{F}},\quad
1\leq\text{i}\leq\text{R}.$ (48)

Finally, recalling that
$\|\mathbf{W}_{\text{i},0}\|_{\text{F}}=\|\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}\|_{\text{F}}=1$,
it is proven that
$\|\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}-{\mathbf{Y}_{\text{i}}}_{\Omega_{\text{i}}}\|_{\text{F}}$
decays exponentially because by induction

$\|\mathbf{W}_{\text{i},\text{p}}\|_{\text{F}}\leq 2^{-\text{p}},\quad
1\leq\text{p}\leq\text{P},\quad 1\leq\text{i}\leq\text{R}.$ (49)

Equating the most constraining upper bound from (49) $2^{-\text{P}}$ to the
upper bound sought $(5\text{R}\gamma)^{-1}$ (30), a lower bound on the number
of partition elements is obtained,
$\text{P}\\!\geq\\!\log_{2}(5\text{R}\\!\gamma)$. Such partitioning guarantees
the fulfillment of
$\|\mathbf{x}_{\text{i}}^{(*)}\mathbf{h}_{\text{i}}^{(*)}{{}^{\text{T}}}-{\mathbf{Y}_{\text{i}}}_{\Omega_{\text{i}}}\|_{\text{F}}\leq(5\text{R}\gamma)^{-1}$
(cf. [27, section II-C] for further discussion on the partitioning).
$\hfill\blacksquare$

###### Lemma 4

Under the condition that
$\text{Q}\geq\mathcal{C}_{\alpha+\log\text{R}}\text{N}^{2}\text{R}\max\\{\mu_{\text{max}}^{2}\text{L}_{\text{i}},\mu_{h}^{2}\text{N}\\}\log^{2}\text{N}$,
then $\|{\mathbf{Y}_{\text{i}}}_{\Omega^{\perp}_{\text{i}}}\|\leq 1/2$ with
probability at least $1-\text{N}^{-\alpha+1}$.

_Proof :_ Expanding the recursion in (42),

$\mathbf{Y}_{\text{i}}=\sum_{\text{p}=1}^{\text{P}}\mathcal{M}_{\text{i},\text{p}}^{*}\left(\boldsymbol{\lambda}_{\text{p}-1}\right).$
(50)

Since ${\mathbf{W}_{\text{j},\text{p}-1}}_{\Omega^{\perp}_{\text{i}}}=0$ (cf.
comment after (45)), it holds that

$\begin{split}\|{\mathbf{Y}_{\text{i}}}_{\Omega^{\perp}_{\text{i}}}\|&=\left\|\sum_{\text{p}=1}^{\text{P}}{\mathcal{M}_{\text{i},\text{p}}^{*}}_{\Omega^{\perp}_{\text{i}}}(\boldsymbol{\lambda}_{\text{p}-1})-{\mathbf{W}_{\text{j},\text{p}-1}}_{\Omega^{\perp}_{\text{i}}}\right\|\\\
&\leq\sum_{\text{p}=1}^{\text{P}}\|\mathcal{M}_{\text{i},\text{p}}^{*}\left(\boldsymbol{\lambda}_{\text{p}-1}\right)-\mathbf{W}_{\text{j},\text{p}-1}\|.\end{split}$
(51)

From here, if we show that
$\|\mathcal{M}_{\text{i},\text{p}}^{*}\left(\boldsymbol{\lambda}_{\text{p}-1}\right)-\mathbf{W}_{\text{j},\text{p}-1}\|\leq
2^{-\text{p}-1}$, then
$\|{\mathbf{Y}_{\text{i}}}_{\Omega^{\perp}_{\text{i}}}\|\leq\sum_{\text{p}=1}^{\text{P}}2^{-\text{p}-1}\leq
1/2$ follows directly.

The upper bound
$\|\mathcal{M}_{\text{i},\text{p}}^{*}(\boldsymbol{\lambda}_{\text{p}-1})-\mathbf{W}_{\text{j},\text{p}-1}\|\leq
2^{-\text{p}-1}$ is based on [27, Lemma 5.14]. Nevertheless, since our model
is based on random unitary $\mathbf{U}_{\text{i}}$ whereas [27] was derived
for random Gaussian, the condition on Q is slightly different. Let us now
illustrate this difference, that is, the results in [27, Lemma 5.14] adapted
to random unitary matrices.

Firstly, the random matrices in the concentration inequality [27, (56)] are
bounded by

$\begin{split}\|\mathcal{Z}_{\ell}\|&\leq\frac{\text{N}\text{R}\sqrt{\text{L}_{\text{i}}\text{N}}}{\text{Q}}2^{-\text{p}+1}\mu_{\text{max}}\
\mu_{h}\\\
&\leq\frac{\text{N}\text{R}}{\text{Q}}2^{-\text{p}+1}\max\\{\mu_{\text{max}}^{2}\text{L}_{\text{i}},\
\mu_{h}^{2}\text{N}\\}.\end{split}$ (52)

The remaining steps follow the development in [27, Lemma 5.14], we only need
to keep track of the extra N factors stemming from the random unitary model.
The center and standard deviation parameters become, respectively,

$R=\max_{\ell\in\Gamma_{\text{p}}}\|\mathcal{Z}_{\ell}\|\leq\frac{\text{N}\text{R}}{\text{Q}}2^{-\text{p}+1}\max\\{\mu_{\text{max}}^{2}\text{L}_{\text{i}},\
\mu_{h}^{2}\text{N}\\}$ (53) $\sigma^{2}\leq\mathcal{C}\
4^{-\text{p}+1}\frac{\text{R}\text{N}^{2}}{Q}\max\\{\mu_{\text{max}}^{2}\text{L}_{\text{i}},\mu_{h}^{2}\\}.$
(54)

Having the vectors $\mathbf{u}_{\text{i},n}$ random unitary instead of
Gaussian, to obtain (53), [27, (90) and (91)] become for deterministic
$\mathbf{q}\in\mathbb{C}^{\text{N}}$, respectively,

$\|(\mathbf{u}_{\text{i},\ell}\mathbf{u}_{\text{i},\ell}^{\text{H}}-\mathbf{I}_{\text{N}})\mathbf{q}\|\leq\|\mathbf{u}_{\text{i},\ell}\mathbf{u}_{\text{i},\ell}^{\text{H}}\mathbf{q}\|+\|\mathbf{q}\|\leq(\text{N}+1)\|\mathbf{q}\|$
(55)
$\mathbb{E}\left[\left(\mathbf{u}_{\text{i},\ell}\mathbf{u}_{\text{i},\ell}^{\text{H}}-\mathbf{I}_{\text{N}}\right)\mathbf{q}\mathbf{q}^{\text{H}}\left(\mathbf{u}_{\text{i},\ell}\mathbf{u}_{\text{i},\ell}^{\text{H}}-\mathbf{I}_{\text{N}}\right)\right]=\text{N}\|\mathbf{q}\|^{2}\mathbf{I}_{\text{N}}-\mathbf{q}\mathbf{q}^{\text{H}}.$
(56)

In addition, to obtain (54), [27, (86)] becomes

$\mathbb{E}\left[\left(\mathbf{u}_{\text{i},\ell}\mathbf{u}_{\text{i},\ell}^{\text{H}}-\mathbf{I}_{\text{N}}\right)^{2}\right]=\left(\text{N}-1\right)\mathbf{I}_{\text{N}}$
(57)

since
$\mathbf{u}_{\text{i},\ell}^{\text{H}}\mathbf{u}_{\text{i},\ell}\\!=\\!\text{N}$
and
$\mathbb{E}\left[\mathbf{u}_{\text{i},\ell}\mathbf{u}_{\text{i},\ell}^{\text{H}}\right]\\!=\\!\mathbf{I}_{\text{N}}$.
Finally, we let
$\text{Q}\geq\mathcal{C}_{\alpha+\log\text{R}}\text{N}^{2}\text{R}\max\\{\mu_{\text{max}}^{2}\text{L}_{\text{i}},\mu_{h}^{2}\text{N}\\}\log^{2}\text{N}$
so that
$\|\mathcal{M}_{\text{i},\text{p}}^{*}\left(\boldsymbol{\lambda}_{\text{p}-1}\right)-\mathbf{W}_{\text{j},\text{p}-1}\|\leq
2^{-\text{p}-1}$ with probability at least $1-\text{N}^{-\alpha+1}$
simultaneously $\forall\ \text{i}\in\\{1,2,\dots,\text{R}\\}$ and
$\text{p}\in\\{1,2,\dots,\text{P}\\}$. $\hfill\blacksquare$

####  -A4 Union bound: Fusing the sufficient conditions

The local isometry property (27),
$\|{\mathcal{M}_{\text{i}}^{*}}_{\Omega_{\text{i}}}{\mathcal{M}_{\text{i}}}_{\Omega_{\text{i}}}-\mathbf{I}_{\Omega_{\text{i}}}\|\\!\leq\\!1/4$,
is proven in [21, Lemma 2]. In particular, it is obtained setting
$\text{Q}=\text{N}$ and $\delta=1/4$ in [21, (30)] and the upper bound holds
with probability at least $1-\text{N}^{-\alpha+1}$ if $\alpha\geq 1$ and

$\begin{split}\alpha\leq\alpha_{1}:=&\min_{\text{i}\in\\{1,2,\dots,\text{R}\\}}\frac{1/16\
(5/2+1/6)^{-1}}{\rho_{\mathbf{\Psi}_{\text{i}}}^{(\text{L}_{\text{i}})}\rho_{\mathbf{U}_{\text{i}}}^{(\text{S}_{\text{i}})}\log\left(2\text{N}\text{L}_{\text{i}}\text{S}_{\text{i}}\right)}=\\\
&\min_{\text{i}\in\\{1,2,\dots,\text{R}\\}}\frac{3}{128}\left(\rho_{\mathbf{\Psi}_{\text{i}}}^{(\text{L}_{\text{i}})}\rho_{\mathbf{U}_{\text{i}}}^{(\text{S}_{\text{i}})}\log\left(2\text{N}\text{L}_{\text{i}}\text{S}_{\text{i}}\right)\right)^{-1}.\end{split}$
(58)

(cf. [21, (31)]).

The coherence between pairs (28), $\mu\leq 1/(4\text{R})$, is proven in
Section -A2. Setting $\text{Q}=\text{N}$ in (2), it holds with probability at
least $1-\text{N}^{-\alpha+1}$ if $\alpha\geq 1$ and

$\begin{split}\alpha\leq\alpha_{2}&:=\min_{j\neq
k}\frac{9}{32\text{R}^{2}\log(2\text{N})}\\\
\Bigg{(}\rho_{\mathbf{\Psi}_{k}}^{(\text{L}_{k})}\rho_{\mathbf{\Psi}_{j}}^{(\text{L}_{j})}&\max\left\\{\rho_{\mathbf{U}_{k}}^{(\text{S}_{k})},\rho_{\mathbf{U}_{j}}^{(\text{S}_{j})}\right\\}+\frac{1}{\text{R}}\kappa_{\mathbf{U}_{k},\mathbf{U}_{j}}^{\left(\text{S}_{k},\text{S}_{j}\right)}\kappa_{\mathbf{\Psi}_{k},\mathbf{\Psi}_{j}}^{\left(\text{L}_{k},\text{L}_{j}\right)}\Bigg{)}^{-1}.\end{split}$
(59)

The third condition (29), $\|\mathcal{M}_{\text{i}}\|\leq\gamma$, is proven in
[21, Lemma 3] and holds with probability at least $1-\text{N}^{-\alpha+1}$,
$\gamma=\sqrt{2\text{N}(\log(2\text{L}_{\text{i}}\text{N})+1)+1}$, $\alpha\geq
1$, and

$\alpha\leq\min_{\text{i}\in\\{1,2,\dots,\text{R}\\}}\left(\rho_{\mathbf{\Psi}_{\text{i}}}^{(\text{L}_{\text{i}})}\log(\text{N})\right)^{-1},$
(60)

which is already fulfilled via $\alpha\leq\alpha_{1}$ in (58).

The first dual certificate condition (30) is based on two of the previous
conditions as shown in Section -A3, the local isometry property (27) and the
coherence between pairs (28).

The second dual certificate and last condition (31),
$\|{\mathcal{M}_{\text{i}}^{*}}_{\Omega_{\text{i}}^{\perp}}(\boldsymbol{\lambda})\|\leq
1/2$, holds with probability at least $1-\text{N}^{-\alpha+1}$ and
$\text{Q}\geq\mathcal{C}_{\alpha+\log\text{R}}\text{N}^{2}\text{R}\max\\{\mu_{\text{max}}^{2}\text{L}_{\text{i}},\mu_{h}^{2}\text{N}\\}\log^{2}\text{N}$
(cf. Lemma 4). Translating the latter condition into an upper bound on
$\alpha$ and simplifying to $\text{Q}=\text{N}$,

$\alpha\\!\leq\\!\alpha_{3}\\!:=\\!\\!\min_{\text{i}\in\\{1,2,\dots,\text{R}\\}}\frac{1}{\mathcal{C}}\\!\left(\text{N}\text{R}\max\\{\mu_{\text{max}}^{2}\text{L}_{\text{i}},\mu_{h}^{2}\text{N}\\}\\!\log^{2}{\text{N}}\right)^{-1}\\!\\!-\log{\text{R}}.$
(61)

By setting $\alpha\leq\min\\{\alpha_{1},\alpha_{2},\alpha_{3}\\}$ as in the
statement of the theorem, we guarantee that all conditions (27)-(31) are
satisfied with high probability.

## References

  * [1] F. Iglesias, S. Segarra, S. Rey-Escudero, A. G. Marques, and D. Ramírez, “Demixing and blind deconvolution of graph-diffused sparse signals,” in _IEEE Intl. Conf. Acoust., Speech and Signal Process. (ICASSP)_. IEEE, 2018, pp. 4189–4193.
  * [2] E. D. Kolaczyk, _Statistical Analysis of Network Data: Methods and Models_. New York, NY: Springer, 2009.
  * [3] W. L. Hamilton, R. Ying, and J. Leskovec, “Representation learning on graphs: Methods and applications,” _arXiv preprint arXiv:1709.05584_ , 2017.
  * [4] W. Hu, M. Fey, M. Zitnik, Y. Dong, H. Ren, B. Liu, M. Catasta, and J. Leskovec, “Open graph benchmark: Datasets for machine learning on graphs,” _arXiv preprint arXiv:2005.00687_ , 2020.
  * [5] P. Djuric and C. Richard, _Cooperative and Graph Signal Processing: Principles and Applications_. Academic Press, 2018.
  * [6] D. Shuman, S. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,” _IEEE Signal Process. Mag._ , vol. 30, no. 3, pp. 83–98, Mar. 2013.
  * [7] A. Ortega, P. Frossard, J. Kovačević, J. M. Moura, and P. Vandergheynst, “Graph signal processing: Overview, challenges, and applications,” _Proc. IEEE_ , vol. 106, no. 5, pp. 808–828, 2018.
  * [8] A. G. Marques, N. Kiyavash, J. M. Moura, D. Van De Ville, and R. Willett, “Graph signal processing: Foundations and emerging directions [from the guest editors],” _IEEE Signal Process. Mag._ , vol. 37, no. 6, pp. 11–13, 2020.
  * [9] A. Anis, A. Gadde, and A. Ortega, “Towards a sampling theorem for signals on arbitrary graphs,” in _IEEE Intl. Conf. Acoust., Speech and Signal Process. (ICASSP)_ , 2014, pp. 3864–3868.
  * [10] S. Chen, R. Varma, A. Sandryhaila, and J. Kovačević, “Discrete signal processing on graphs: Sampling theory,” _IEEE Trans. Signal Process._ , vol. 63, no. 24, pp. 6510–6523, 2015.
  * [11] A. G. Marques, S. Segarra, G. Leus, and A. Ribeiro, “Sampling of graph signals with successive local aggregations,” _IEEE Trans. Signal Process._ , vol. 64, no. 7, pp. 1832–1843, 2016.
  * [12] M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst, “Geometric deep learning: going beyond euclidean data,” _IEEE Signal Process. Mag._ , vol. 34, no. 4, pp. 18–42, 2017.
  * [13] M. Defferrard, X. Bresson, and P. Vandergheynst, “Convolutional neural networks on graphs with fast localized spectral filtering,” in _Advances Neural Info. Process. Systems_ , 2016, pp. 3844–3852.
  * [14] F. Gama, A. G. Marques, G. Leus, and A. Ribeiro, “Convolutional neural network architectures for signals supported on graphs,” _IEEE Trans. Signal Process._ , vol. 67, no. 4, pp. 1034–1049, 2019.
  * [15] A. Sandryhaila and J. Moura, “Discrete signal processing on graphs,” _IEEE Trans. Signal Process._ , vol. 61, no. 7, pp. 1644–1656, Apr. 2013\.
  * [16] A. Sandryhaila and J. M. F. Moura, “Finite-time distributed consensus through graph filters,” in _IEEE Intl. Conf. Acoust., Speech and Signal Process. (ICASSP)_ , 2014, pp. 1080–1084.
  * [17] S. Segarra, A. G. Marques, and A. Ribeiro, “Optimal graph-filter design and applications to distributed linear network operators,” _IEEE Trans. Signal Process._ , vol. 65, no. 15, pp. 4117–4131, 2017.
  * [18] M. Coutino, E. Isufi, and G. Leus, “Advances in distributed graph filtering,” _IEEE Trans. Signal Process._ , vol. 67, no. 9, pp. 2320–2333, 2019.
  * [19] A. G. Marques, S. Segarra, G. Leus, and A. Ribeiro, “Stationary graph processes and spectral estimation,” _IEEE Trans. Signal Process._ , vol. 65, no. 22, pp. 5911–5926, 2017.
  * [20] A. G. Marques, S. Segarra, and G. Mateos, “Signal processing on directed graphs: The role of edge directionality when processing and learning from network data,” _IEEE Signal Processing Magazine_ , vol. 37, no. 6, pp. 99–116, 2020.
  * [21] S. Segarra, G. Mateos, A. G. Marques, and A. Ribeiro, “Blind identification of graph filters,” _IEEE Trans. Signal Process._ , vol. 65, no. 5, pp. 1146–1159, 2016.
  * [22] Y. Zhu, F. Iglesias, A. G. Marques, and S. Segarra, “Estimating network processes via blind identification of multiple graph filters,” _IEEE Trans. Signal Process._ , vol. 68, pp. 3049–3063, 2020.
  * [23] D. Ramírez, A. G. Marques, and S. Segarra, “Graph-signal reconstruction and blind deconvolution for diffused sparse inputs,” in _IEEE Intl. Conf. Acoust., Speech and Signal Process. (ICASSP)_ , Mar. 2017, pp. 4104–4108.
  * [24] D. Thanou, X. Dong, D. Kressner, and P. Frossard, “Learning heat diffusion graphs,” _IEEE Trans. Signal Inf. Process. Netw._ , vol. 3, no. 3, pp. 484–499, 2017.
  * [25] S. Rey-Escudero, F. Iglesias, C. Cabrera, and A. G. Marques, “Sampling and reconstruction of diffused sparse graph signals from successive local aggregations,” _IEEE Signal Process. Lett._ , vol. 26, no. 8, pp. 1142–1146, 2019.
  * [26] G. Cutura, B. Li, A. Swami, and S. Segarra, “Deep demixing: Reconstructing the evolution of epidemics using graph neural networks,” _arXiv preprint arXiv:2011.09583_ , 2020.
  * [27] S. Ling and T. Strohmer, “Blind deconvolution meets blind demixing: Algorithms and performance bounds,” _IEEE Transactions on Information Theory_ , vol. 63, no. 7, pp. 4497–4520, 2017.
  * [28] M. Fazel, H. Hindi, and S. P. Boyd, “A rank minimization heuristic with application to minimum order system approximation,” in _American Control Conf._ , 2001, pp. 4734–4739.
  * [29] J. A. Tropp, “Just relax: Convex programming methods for identifying sparse signals in noise,” _IEEE Trans. Inf. Theory_ , vol. 52, no. 3, pp. 1030–1051, 2006.
  * [30] B. Recht, M. Fazel, and P. A. Parrilo, “Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization,” _Society for Industrial and Applied Mathematics Review_ , vol. 52, no. 3, pp. 471–501, Jan. 2010.
  * [31] S. Foucart and H. Rauhut, _A Mathematical Introduction to Compressive Sensing_. New York, NY, USA: Springer, 2013\.
  * [32] A. Flinth, “Sparse blind deconvolution and demixing through $\ell_{1,2}$-minimization,” _Advances in Computational Mathematics_ , vol. 44, no. 1, pp. 1–21, 2018.
  * [33] M. Fazel, H. Hindi, and S. P. Boyd, “Log-det heuristic for matrix rank minimization with applications to Hankel and Euclidean distance matrices,” in _American Control Conf._ , 2003.
  * [34] B. Bollobás, _Random Graphs_. Springer, 1998.
  * [35] R. Albert and A.-L. Barabási, “Statistical mechanics of complex networks,” _Reviews of modern physics_ , vol. 74, no. 1, p. 47, 2002.
  * [36] P. Hagmann, L. Cammoun, X. Gigandet, R. Meuli, C. J. Honey, V. J. Wedeen, and O. Sporns, “Mapping the structural core of human cerebral cortex,” _PLoS Biol_ , vol. 6, no. 7, p. e159, 2008.
  * [37] W. W. Zachary, “An information flow model for conflict and fission in small groups,” _Journal of anthropological research_ , vol. 33, no. 4, pp. 452–473, 1977.
  * [38] J. Mei and J. M. Moura, “Signal processing on graphs: Estimating the structure of a graph,” in _IEEE Intl. Conf. Acoust., Speech and Signal Process. (ICASSP)_. IEEE, 2015, pp. 5495–5499.
  * [39] S. Segarra, A. G. Marques, G. Mateos, and A. Ribeiro, “Network topology inference from spectral templates,” _IEEE Trans. Signal Inf. Process. Netw._ , vol. 3, no. 3, pp. 467–483, 2017.
  * [40] Y. Zhu, M. T. Schaub, A. Jadbabaie, and S. Segarra, “Network inference from consensus dynamics with unknown parameters,” _IEEE Trans. Signal Inf. Process. Netw._ , vol. 6, pp. 300–315, 2020.