# Robust Graph Filter Identification and Graph Denoising from Signal
Observations

Samuel Rey, , Victor M. Tenorio, and Antonio G. Marques S. Rey, V. M. Tenorio
and A. G. Marques are with the Department of Signal Theory and Comms., King
Juan Carlos University, Madrid, Spain,
<EMAIL_ADDRESS>conference version of this manuscript with some preliminary results was
published in [1]. Work in this paper was partially supported by the Spanish
Grants SPGRAPH (PID2019-105032GB-I00), FPU17/04520 and FPU20/05554. The
authors would like to thank Pablo Espana-Gutierrez for his help reviewing the
literature and developing the early versions of the convergence results.

###### Abstract

When facing graph signal processing tasks, the workhorse assumption is that
the graph describing the support of the signals is known. However, in many
relevant applications _the available graph suffers from observation errors and
perturbations_. As a result, any method relying on the graph topology may
yield suboptimal results if those imperfections are ignored. Motivated by
this, we propose a novel approach for handling perturbations on the links of
the graph and apply it to the problem of robust graph filter (GF)
identification from input-output observations. Different from existing works,
we formulate a non-convex optimization problem that operates in the vertex
domain and jointly performs GF identification and graph denoising. As a
result, on top of learning the desired GF, an estimate of the graph is
obtained as a byproduct. To handle the resulting bi-convex problem, we design
an algorithm that blends techniques from alternating optimization and
majorization minimization, showing its convergence to a stationary point. The
second part of the paper i) generalizes the design to a robust setup where
several GFs are jointly estimated, and ii) introduces an alternative
algorithmic implementation that reduces the computational complexity. Finally,
the detrimental influence of the perturbations and the benefits resulting from
the robust approach are numerically analyzed over synthetic and real-world
datasets, comparing them with other state-of-the-art alternatives.

###### Index Terms:

Graph Filter Identification, Graph Denoising, Robust Graph Signal Processing,
Graph Perturbations.

##  I Introduction

Nowadays, a significant number of datasets are defined over an irregular
(heterogeneous) support that can be conveniently represented by a graph. As a
result, the data at hand can be readily understood as graph signals
(alternatively, network processes) whose structure and properties depend on
the topology of the generating graph. Illustrative examples include
measurements from power, communications, social, biological, or financial
networks, to name a few [2, 3, 4, 5, 6]. Characterizing and modeling graph and
network processes entails a prevalent and relevant task that not only enhances
our understanding of the data at hand, but also opens the door to more
sophisticated processing and knowledge extraction schemes. A popular approach
within the graph signal processing (GSP) framework [5, 7, 8, 9, 10] is to
represent the generative process by a _graph filter_ (GF) and, then, model the
data as graph signals generated by applying a (low-pass, smooth, bandlimited…)
GF to a simple (sparse, white, constant…) input. GFs are topology-aware linear
operators whose output can be interpreted as the outcome of a network
diffusion or spreading process [11, 12, 13]. GFs can be expressed as
polynomials of the graph-shift operator (GSO), a matrix encoding the topology
of the graph, and on top of its theoretical interest, the task of GF
identification is practically relevant to, e.g., understanding the dynamics of
network diffusion processes [12, 14, 7], as well as explaining the structure
of real-world datasets [15, 16, 17].

Since GSP is a relatively recent area of research, it is not surprising that
most GSP works focus on how to harness the graph structure while assuming that
the topology of the graph is perfectly known. Nevertheless, this assumption is
unlikely to hold in many practical setups, where graphs suffer from
_imperfections and perturbations_. When networks are given explicitly,
perturbations may be due to observational noise and errors (e.g., link
failures in power or wireless networks [18]). When in lieu of physical
entities, the graphs model (statistical) pairwise relationships among the
observed variables, they need to be inferred from the data [19, 20, 21, 22].
Since this is a challenging (oftentimes ill-posed) task, the estimated graphs
inherit the imperfections (limitations) of the graph learning scheme adopted
(e.g., the thresholding operation implemented in correlation networks [3]).
Intuitively, the presence of perturbations hinders any GSP scheme or GSP tool
applied to the data. While this is true regardless of the task and hand, it is
even more relevant for those involving spectral transforms and GFs, since
eigenvectors and high-order matrix polynomials are more sensitive to errors in
the matrix codifying the graph.

To address these challenges, this paper investigates the problem of estimating
a GF from input-output signal pairs assuming that both the signals and the
supporting graph have errors. The proposed approach is formulated in the
vertex domain, avoiding the numerical instability of computing large
polynomials and, at the same time, bypassing the challenges associated with
robust _spectral_ graph theory. We recast the robust estimation as a joint
optimization problem where the GF identification objective is augmented with a
graph-denoising regularizer, so that, on top of the desired GF, we also obtain
an enhanced estimate of the supporting graph. The joint formulation leads to a
non-convex bi-convex optimization problem, for which a provably-convergent
efficient (alternating minimization) algorithm able to find an approximate
solution is developed. Furthermore, to address scenarios where multiple GFs
are present (e.g., when dealing with vector autoregressive (AR) spatio-
temporal processes or in setups where nodes collect multi-feature vector
measurements), we generalize our framework so that multiple GFs, all defined
over the same graph, are jointly identified.

Despite their theoretical and practical relevance, the number of robust GSP
works is limited, due in part to the challenges emanating from the presence of
graph perturbations [23, 24, 25, 26]. In the spectral domain,[23] employs a
small perturbation analysis to study the impact of perturbations in the
spectrum of the graph Laplacian. In the vertex domain, [27, 25] postulates a
graphon-based perturbation model and analyzes how perturbations affect GFs of
order one. More recently, [26] combines structural equation models (SEMs) with
total least squares (TLS) to jointly infer the GF and the perturbations when
the observed data is explained by a SEM. A different robust alternative is
presented in [28], where the support of the graph is assumed to be known and
the goal is to estimate the weights of the network topology and the
coefficients of the GF. The resultant problem is non-convex and the authors
adopt a sequential convex programming (SCP) approach to solve it. Finally, the
presence of perturbations has also been considered in non-linear GSP tasks. An
alternative definition of GFs robust to perturbations is proposed in [29], and
the transferability of GFs when employed in graph neural networks is studied
in [30, 31, 32].

Contributions and outline. We close the section by summarizing the
organization and contributions of the manuscript. After reviewing preliminary
GSP concepts in Sec. II, we analyze the influence of edge perturbations in
polynomial GFs and state the robust GF identification problem in Sec. III.
After that, our main contributions are:

  1. 1.

We formulate a non-convex optimization problem to jointly estimate the graph
and the GF, develop an alternating optimization algorithm to solve it and
prove its convergence to a stationary point (Sec. IV).

  2. 2.

We consider a generalization where several GFs are jointly estimated by
exploiting the fact that they are polynomials of the same GSO (Sec. V).

  3. 3.

We propose an efficient implementation of the GF identification algorithm to
handle graphs with a large number of nodes (Sec. VI).

The effectiveness of the proposed algorithms is evaluated numerically in Sec.
VII, and some concluding remarks are provided in Sec. VIII. Last but not
least, while we focus on GF identification from input-output pairs, the
approach put forth in this paper can be generalized to other GSP tasks, which
is a research path we plan to pursue in the near future.

##  II GSP preliminaries

Graphs and graph signals. Consider a directed graph
${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}})$ formed by the set of $N$ nodes
(vertices) collected in ${\mathcal{V}}$ and the set of edges
${\mathcal{E}}\subset{\mathcal{V}}\times{\mathcal{V}}$, such that
$(i,j)\in{\mathcal{E}}$ if node $i$ is connected to node $j$. The topology of
${\mathcal{G}}$ can be represented by the adjacency matrix
${\mathbf{A}}\in{\mathbb{R}}^{N\times N}$, a sparse matrix whose elements
$A_{ij}$ are non-zero if and only if $(i,j)\in{\mathcal{E}}$. When
${\mathcal{G}}$ is weighted, the entries $A_{ij}$ capture the strength of the
link between nodes $i$ and $j$. Otherwise, the elements of ${\mathbf{A}}$ are
binary. Together with the graph ${\mathcal{G}}$, we focus on modeling the data
as signals defined on the nodes in ${\mathcal{V}}$. Formally, a (nodal) _graph
signal_ is a function from the vertex set to the real field
$x:{\mathcal{V}}\to{\mathbb{R}}$, which can be alternatively represented as an
$N$-dimensional vector ${\mathbf{x}}\in{\mathbb{R}}^{N}$ whose $i$-th entry
$x_{i}$ denotes the value of the signal at node $i$. The foundational
assumption of GSP is that the properties of the ${\mathbf{x}}$ and the
topology of ${\mathcal{G}}$ are related. As a simple example, if
${\mathcal{G}}$ captures the similarity between nodes and $A_{ij}$ is high,
then the values of $x_{i}$ and $x_{j}$ are expected to be akin to each other.

Graph-shift operator (GSO). The GSO associated with a graph ${\mathcal{G}}$ of
$N$ nodes is as a generic matrix ${\mathbf{S}}\in{\mathbb{R}}^{N\times N}$
that: i) captures the topology of the underlying graph and ii) represents a
_local_ and _linear_ transformation that can be applied to graph signals
defined over ${\mathcal{G}}$. The entries of ${\mathbf{S}}$ satisfy
$S_{ij}\neq 0$ only if $(i,j)\in{\mathcal{E}}$ or $i=j$ and, as a result, the
application of the GSO to a graph signal involves mixing values among one-hop
neighborhoods. Two typical choices for the GSO are the adjacency matrix
${\mathbf{A}}$ and the graph Laplacian
${\mathbf{L}}:=\\!\text{diag}({\mathbf{A}}\textbf{1})\\!-\\!{\mathbf{A}}$ [5,
7], where $\text{diag}(\cdot)$ transforms a vector into a diagonal matrix and
1 is the vector of all ones. We assume that ${\mathbf{S}}$ is diagonalized as
${\mathbf{S}}={\mathbf{V}}\text{diag}(\bm{\lambda}){\mathbf{V}}^{-1}$, where
${\mathbf{V}}$ is an $N\times N$ matrix collecting the eigenvectors of
${\mathbf{S}}$, and vector $\bm{\lambda}$ collects its eigenvalues. The matrix
${\mathbf{V}}^{-1}$ is commonly adopted as the Graph Fourier Transform (GFT)
for graph signals with ${\tilde{\mathbf{x}}}={\mathbf{V}}^{-1}{\mathbf{x}}$
denoting the graph frequency representation of ${\mathbf{x}}$ [6].

Graph filtering. GFs are topology-aware operators whose inputs and outputs are
graph signals. More specifically, GFs implement a linear transformation that
can be written as a polynomial of ${\mathbf{S}}$

${\mathbf{H}}=\sum_{r=0}^{N-1}h_{r}{\mathbf{S}}^{r}={\mathbf{V}}\text{diag}(\bm{\Psi}{\mathbf{h}}){\mathbf{V}}^{-1}={\mathbf{V}}\text{diag}({\tilde{\mathbf{h}}}){\mathbf{V}}^{-1},$
(1)

where ${\mathbf{h}}=[h_{0},...,h_{N-1}]$ is the vector collecting the GF
coefficients. The $N\times N$ Vandermonde matrix $\bm{\Psi}$ defined as
$\Psi_{ij}:=\lambda_{i}^{j-1}$ represents the GFT for GFs, and thus,
${\tilde{\mathbf{h}}}=\bm{\Psi}{\mathbf{h}}$ is the vector of size $N$
representing the frequency response of ${\mathbf{H}}$ [6, 12]. Since
${\mathbf{S}}^{r}$ encodes the $r$-hop neighborhood of the graph, a graph
signal given by
${\mathbf{y}}=\sum_{r=0}^{N-1}h_{r}{\mathbf{S}}^{r}{\mathbf{x}}={\mathbf{H}}{\mathbf{x}}$
can be interpreted as a version of the input signal ${\mathbf{x}}$ diffused
across $N\\!-\\!1$ neighborhoods with $h_{r}$ being the coefficients of the
linear combination [12]. Moreover, there are scenarios where $h_{r}=0$ for
$r\geq R$. In those cases, the order of the filter is $R$ and, if more
convenient, ${\mathbf{h}}$ and $\bm{\Psi}$ can be redefined so that the number
of elements of ${\mathbf{h}}$ (columns of $\bm{\Psi}$) is $R$ in lieu of $N$.

Graph stationarity. A zero-mean _random_ graph signal ${\mathbf{x}}$ is said
to be stationary on ${\mathcal{G}}$ if its covariance matrix
${\mathbf{C}}_{\mathbf{x}}=\mathbb{E}[{\mathbf{x}}{\mathbf{x}}^{\top}]$ is a
positive-semidefinite polynomial of the GSO [33]. A common example of
stationary graph signals arises when ${\mathbf{x}}$ is the output of a linear
graph diffusion process whose input (initial condition) is a white signal
${\mathbf{w}}$ and whose diffusion dynamics can be accurately represented by a
GF. Mathematically, if we have that ${\mathbf{x}}={\mathbf{H}}{\mathbf{w}}$
with ${\mathbf{H}}$ being a GF [cf. (1)] and
${\mathbf{C}}_{\mathbf{w}}=\mathbb{E}[{\mathbf{w}}{\mathbf{w}}^{\top}]={\mathbf{I}}$,
it follows that the covariance of ${\mathbf{x}}$ is
${\mathbf{C}}_{\mathbf{x}}={\mathbf{H}}{\mathbf{H}}^{\top}=\sum_{r=0,r^{\prime}=0}^{N-1,N-1}h_{r}h_{r^{\prime}}{\mathbf{S}}^{r+r^{\prime}}$.
Since the latter is a polynomial of the GSO, it follows that ${\mathbf{x}}$ is
stationary on ${\mathcal{G}}$.

###  II-A GF identification from input-output pairs

In the context of linear operators, let us consider that we observe $M$ input
and output pairs ${\mathbf{X}}:=[{\mathbf{x}}_{1},...,{\mathbf{x}}_{M}]$ and
${\mathbf{Y}}:=[{\mathbf{y}}_{1},...,{\mathbf{y}}_{M}]$ whose relation is
given by

${\mathbf{Y}}={\mathbf{H}}{\mathbf{X}}+{\mathbf{W}},$ (2)

with ${\mathbf{W}}$ being a zero-mean random matrix (typically assumed to have
i.i.d. entries) that accounts for noisy measurements and model inaccuracies.
Leveraging (2), the GF identification task amounts to using the input-output
pairs to estimate ${\mathbf{H}}$ under the model in (1), which, if the GSO
${\mathbf{S}}$ is known, boils down to estimating the GF coefficients
collected in ${\mathbf{h}}\in{\mathbb{R}}^{N}$.

Hence, we can approach the GF identification task in the node domain by
solving the convex problem

$\min_{\mathbf{h}}\big{\|}{\mathbf{Y}}-\sum_{r=0}^{N-1}h_{r}{\mathbf{S}}^{r}{\mathbf{X}}\big{\|}_{F}^{2}.$
(3)

Leveraging the frequency definition of GFs in (1), we use the GFT matrices
${\mathbf{V}}^{-1}$ and $\bm{\Psi}$ to rewrite the least-squares (LS) cost in
(3) and obtain its (closed-form) solution as

$\\!\\!\\!{\hat{\mathbf{h}}}\\!=\\!\operatornamewithlimits{argmin}_{\mathbf{h}}\|{\mathrm{vec}}({\mathbf{Y}})\\!-\\!(\\!({\mathbf{V}}^{-1}{\mathbf{X}})^{\top}\\!\\!\odot\\!{\mathbf{V}})\bm{\Psi}{\mathbf{h}}\|_{2}^{2}\\!=\\!\bm{\Theta}^{\dagger}{\mathrm{vec}}({\mathbf{Y}}),$
(4)

where ${\mathrm{vec}}(\cdot)$ denotes the vectorization operation,
${\mathbf{V}}^{-1}{\mathbf{X}}$ is the frequency representation of the input
signals, $\odot$ denotes the Khatri–Rao product, $\bm{\Psi}$ is the GFT
Vandermonde matrix,
$\bm{\Theta}:=(({\mathbf{V}}^{-1}{\mathbf{X}})^{\top}\odot{\mathbf{V}})\bm{\Psi}$,
and † is the pseudoinverse operator.

From (4) we observe that estimating ${\mathbf{H}}$ is straightforward under
the assumptions of: i) $\bm{\Theta}$ being full rank (i.e., the inputs are
sufficiently rich) and ii) ${\mathbf{S}}$ being perfectly known. However, the
(critical) assumption in ii) does not hold in most practical settings. The
remainder of the paper approaches the GF identification problem assuming
imperfect knowledge of the GSO.

##  III GF identification with imperfect graph knowledge

This section introduces and discusses the problem of estimating a GF
${\mathbf{H}}=\sum_{r=0}^{N-1}h_{r}{\mathbf{S}}^{r}$ from noisy input-output
signal pairs $({\mathbf{X}}\in{\mathbb{R}}^{N\times
M},{\mathbf{Y}}\in{\mathbb{R}}^{N\times M})$ assuming that we have access to
an _imperfect_ GSO ${\bar{\mathbf{S}}}\in{\mathbb{R}}^{N\times N}$, which can
be modeled as

${\bar{\mathbf{S}}}={\mathbf{S}}+\bm{\Delta},$ (5)

where ${\mathbf{S}}\in{\mathbb{R}}^{N\times N}$ represents the true GSO and
$\bm{\Delta}\in{\mathbb{R}}^{N\times N}$ is a _perturbation matrix_. Before
discussing models for the perturbation matrix, we find illustrative to
demonstrate the impact of $\bm{\Delta}$ on the GSP problem at hand.

Figure 1: Absolute error for different powers of the matrix ${\mathbf{S}}$ and
its perturbed version ${\bar{\mathbf{S}}}$. The true GSO is the adjacency
matrix of an Erdős-Rényi graph with link probability of 0.15, and
${\bar{\mathbf{S}}}$ is perturbed by creating and destroying links
independently with a probability of $0.05$.

As pointed out in the introduction, the presence of uncertainties in the
topology of ${\mathcal{G}}$ is particularly relevant when dealing with GFs.
Indeed, due to the polynomial definition of ${\mathbf{H}}$, even small
perturbations can lead to significant errors when ${\bar{\mathbf{S}}}$ (and
not ${\mathbf{S}}$) is used as the true GSO. To see this more clearly, Fig. 1
provides an example that illustrates how the errors encoded in $\bm{\Delta}$
propagate for different matrix powers, demonstrating that the discrepancies
between ${\bar{\mathbf{S}}}^{r}$ and ${\mathbf{S}}^{r}$ increase swiftly as
the power $r$ grows. More rigorously, let $C$ be a positive constant such that
$\|{\mathbf{S}}\|\leq C$ and $\|{\bar{\mathbf{S}}}\|\leq C$, and define
${\bar{\mathbf{H}}}:=\sum_{r=0}^{N-1}h_{r}{\bar{\mathbf{S}}}^{r}$. Then, the
error generated by the perturbations is upper-bounded by

$\|{\bar{\mathbf{H}}}-{\mathbf{H}}\|\leq\sum_{r=1}^{N-1}|h_{r}|\|{\bar{\mathbf{S}}}^{r}-{\mathbf{S}}^{r}\|\leq\sum_{r=1}^{N-1}|h_{r}|rC^{r-1}\|\bm{\Delta}\|,$
(6)

where the last inequality follows from [30, Lemma 3]. In words, the maximum
difference between the true ${\mathbf{H}}$ and the perturbed
${\bar{\mathbf{H}}}$ increases _exponentially_ with the degree of the GF.

From the previous discussion, it is not surprising that the imperfect
knowledge of the graph topology is also relevant when estimating the filter
coefficients. In fact, ignoring the errors in $\bm{\Delta}$ and attempting to
estimate ${\mathbf{h}}$ solving (4) when ${\bar{\mathbf{S}}}$ is used in lieu
of the true (unknown) ${\mathbf{S}}$ leads to a poor solution, as we
illustrate numerically in Sec. VII. Motivated by this, we approach the GF
identification problem from a robust perspective by taking into account the
imperfect knowledge of the GSO. The resultant robust estimation task is
formally stated next.

######  Problem 1.

Let ${\mathcal{G}}$ be a graph with $N$ nodes, let
${\mathbf{S}}\in{\mathbb{R}}^{N\times N}$ be the true (unknown) GSO associated
with ${\mathcal{G}}$, and let ${\bar{\mathbf{S}}}\in{\mathbb{R}}^{N\times N}$
be the perturbed (observed) GSO. Moreover, let
${\mathbf{X}}\in{\mathbb{R}}^{N\times M}$ and
${\mathbf{Y}}\in{\mathbb{R}}^{N\times M}$ be a pair of matrices collecting $M$
observed input and output signals defined over ${\mathcal{G}}$ and related by
the model in (2). Our goal is to use the triplet
$({\mathbf{X}},{\mathbf{Y}},{\bar{\mathbf{S}}})$ to: i) learn the GF
${\mathbf{H}}$ that best fits the model in (2) and ii) recover an enhanced
estimation of ${\mathbf{S}}$. To that end, we make the following assumptions:
(AS1) ${\mathbf{H}}$ is a polynomial of ${\mathbf{S}}$ [cf. (1)].
(AS2) ${\mathbf{S}}$ and ${\bar{\mathbf{S}}}$ are close according to some
metric $d({\mathbf{S}},{\bar{\mathbf{S}}})$, i.e., the observed perturbations
are “small” in some sense.

On top of the previous two assumptions, we also consider that the norm of the
noise observation matrix ${\mathbf{W}}$ in (2) is small, which is a workhorse
assumption in this type of problems. Similar to standard GF identification
approaches, (AS1) limits the degrees of freedom of the linear operator in (2).
However, the fact of the true ${\mathbf{S}}$ being unknown adds uncertainty to
the problem and, as a result, additional signal observations are required to
achieve an identification performance comparable to the one obtained when
${\bar{\mathbf{S}}}={\mathbf{S}}$. Regarding the recovery of the true GSO,
(AS2) accounts for the hypothesis that ${\bar{\mathbf{S}}}$ is a perturbed
observation of ${\mathbf{S}}$ and, hence, matrices ${\mathbf{S}}$ and
${\bar{\mathbf{S}}}$ are not extremely different. Note that this guarantees
that “some” information about the true GSO is available, so that (AS1) can be
effectively leveraged. While not exploited in our formulation, additional
assumptions constraining the GSO could also be incorporated into the problem.
Finally, the metric $d(\cdot,\cdot)$ employed to quantify the similarity
between ${\mathbf{S}}$ and ${\bar{\mathbf{S}}}$ should depend on the model for
the perturbation $\bm{\Delta}$, a subject that is briefly discussed next.

###  III-A Modeling graph perturbations

The development and analysis of graph perturbation models that combine
practical relevance and analytical tractability constitutes an interesting yet
challenging open line of research [27, 25]. Due to its flexibility and
tractability, this paper considers an additive perturbation model [cf. (5)],
so that the focus is constrained to understanding the structural (statistical)
properties of matrix $\bm{\Delta}={\bar{\mathbf{S}}}-{\mathbf{S}}$.

Consider first the case where perturbations only _create or destroy links_
independently. If ${\mathcal{G}}$ is an _unweighted graph_ , a simple approach
is to consider perturbations modeled as independent Bernoulli variables with
possibly different creation/destruction probabilities. In this case, the
entries of $\bm{\Delta}$ would be

$\Delta_{ij}=\left\\{\begin{array}[]{rl}1&\mathrm{if}\;\mathrm{link}\;(i,j)\;\mathrm{is}\;\mathrm{created},\\\
-1&\mathrm{if}\;\mathrm{link}\;(i,j)\;\mathrm{is}\;\mathrm{destroyed},\\\
0&\mathrm{otherwise}.\end{array}\right.$ (7)

Since $\bm{\Delta}$ models the creation and destruction of links, it is worth
noting that $\Delta_{ij}=1$ only if $S_{ij}=0$ and $\Delta_{ij}=-1$ only if
$S_{ij}=1$. In the more general case of ${\mathcal{G}}$ being a _weighted
graph_ , $\Delta_{ij}=-S_{ij}$ destroys an existing link while $\Delta_{ij}=z$
creates a new link. Here, $z$ is a random variable sampled from a particular
distribution (typically mimicking the weight distribution of the true
${\mathbf{S}}$). When facing this type of perturbations, a suitable distance
function is the $\ell_{0}$ norm

$d({\mathbf{S}},{\bar{\mathbf{S}}})=\|{\mathbf{S}}-{\bar{\mathbf{S}}}\|_{0},$
(8)

with the $\ell_{1}$ norm $\|{\mathbf{S}}-{\bar{\mathbf{S}}}\|_{1}$ being a
prudent convex relaxation.

Alternatively, rather than creating or destroying links, perturbations may
represent uncertainty over the edge weights. This entails the support of
matrix $\bm{\Delta}$ matching that of ${\mathbf{S}}$ and ${\bar{\mathbf{S}}}$,
and the non-zero entries of $\bm{\Delta}$ being sampled from a distribution
that models the observation noise. For example, if the noise is zero-mean,
Gaussian and white, it holds that $\Delta_{ij}\sim{\mathcal{N}}(0,\sigma^{2})$
when ${\mathbf{S}}_{ij}\neq 0$ and $\Delta_{ij}=0$ when ${\mathbf{S}}_{ij}=0$.
Under this setting, an appropriate distance metric is given by

$d({\mathbf{S}},{\bar{\mathbf{S}}})=\|{\mathbf{S}}_{\mathcal{E}}-{\bar{\mathbf{S}}}_{\mathcal{E}}\|_{2}^{2},$
(9)

where ${\mathbf{S}}_{\mathcal{E}}$ and ${\bar{\mathbf{S}}}_{\mathcal{E}}$ only
select the non-zero entries (edges) in ${\mathbf{S}}$ and
${\bar{\mathbf{S}}}$. Additionally, one can have setups where the two types of
perturbations are simultaneously present. That is, perturbations may create
and destroy links while the actual value of the existing links is also
uncertain. In such a case, a combination of $\ell_{1}$ and $\ell_{2}$ norms
like in elastic nets [34] is adequate.

The models previously described only consider the perturbation of edges in an
independent fashion. However, there may be scenarios where the perturbations
are correlated. Consider for example a communication network. If the power
supply of a node stalls, the signal-to-noise ratio of all its links will be
poor, and hence, links involving that node will be more likely to fail.
Perturbations dependent across links can be modeled by means of a multivariate
correlated Bernoulli distribution, an Ising model, or more sophisticated
random graph models [35]. When prior information about the dependence of the
perturbations is available, it can be incorporated into the function
$d({\mathbf{S}},{\bar{\mathbf{S}}})$ to better extract the information encoded
in ${\bar{\mathbf{S}}}$.

##  IV Robust GF identification

This section presents the optimization problem and the proposed algorithm to
estimate ${\mathbf{H}}$ and ${\mathbf{S}}$ under the setting described in 1.
Given the matrices ${\mathbf{X}}$, ${\mathbf{Y}}$, and ${\bar{\mathbf{S}}}$,
we approach the robust GF identification task by means of the following non-
convex optimization

$\displaystyle\\!{\hat{\mathbf{H}}},{\hat{\mathbf{S}}}=\operatornamewithlimits{argmin}_{{\mathbf{H}},{\mathbf{S}}}$
$\displaystyle\;\;\|{\mathbf{Y}}-{\mathbf{H}}{\mathbf{X}}\|_{F}^{2}+\lambda
d({\mathbf{S}},{\bar{\mathbf{S}}})+\beta\|{\mathbf{S}}\|_{0}$
$\displaystyle\\!\hskip 34.14322pt\mathrm{\;\;s.\;to:}$
$\displaystyle\;\;{\mathbf{S}}\in{\mathcal{S}},\;\;{\mathbf{S}}{\mathbf{H}}={\mathbf{H}}{\mathbf{S}},$
(10)

where $\mathrm{s.\;to}$ stands for $\mathrm{subject\;to}$. The first term in
the objective promotes the linear input-output relation in (2), encouraging
the norm of ${\mathbf{W}}={\mathbf{Y}}-{\mathbf{H}}{\mathbf{X}}$ to be small.
The use of the Frobenius norm is well-justified when the observation noise is
Gaussian and white, but other types of noise could be accommodated by using a
different norm. The second term incorporates the assumption (AS2) as a
regularizer to obtain an estimate ${\hat{\mathbf{S}}}$ that is related to the
given GSO ${\bar{\mathbf{S}}}$. The $\ell_{0}$ norm in the third term accounts
for the fact of ${\mathbf{S}}$ being sparse. Clearly, if additional
information about ${\mathbf{S}}$ is available, it can be incorporated into
(IV), either as a regularizer (e.g., a statistical prior quantifying the log-
likelihood of a class of GSOs) or as a constraint that _must_ be satisfied
(e.g., the GSO being symmetric). The latter is indeed the role of
${\mathbf{S}}\in{\mathcal{S}}$ in (2), with ${\mathcal{S}}$ representing a
(desired) family of GSOs such as the set of adjacency matrices with no self-
loops (${\mathcal{S}}$ is the set of matrices with non-negative entries whose
diagonal entries are zero) or the set of combinatorial graph Laplacians
(matrices with non-positive off-diagonal entries and zero row-sum). Finally,
the (key) constraint ${\mathbf{S}}{\mathbf{H}}={\mathbf{H}}{\mathbf{S}}$
captures the fact of ${\mathbf{H}}$ being a polynomial of ${\mathbf{S}}$ and
not of ${\bar{\mathbf{S}}}$ (AS1). Note first that the constraint is
pertinent, if ${\mathbf{H}}$ is a polynomial of ${\mathbf{S}}$, then
${\mathbf{H}}$ and ${\mathbf{S}}$ have the same eigenvectors and, as a result,
their product commutes [12]. More importantly for the GF-identification at
hand, when the GSO is perfectly known the model
${\mathbf{H}}=h_{0}{\mathbf{I}}+h_{1}{\mathbf{S}}+...+h_{N-1}{\mathbf{S}}^{N-1}$
is linear in the unknown ${\mathbf{h}}$. As a result, a formulation that
estimates ${\mathbf{h}}$ directly (as carried out in classical non-robust
approaches) is well-motivated. However, when both ${\mathbf{h}}$ and
${\mathbf{S}}$ are unknown, the model
${\mathbf{H}}=h_{0}{\mathbf{I}}+h_{1}{\mathbf{S}}+...+h_{N-1}{\mathbf{S}}^{N-1}$
is highly non-linear in ${\mathbf{S}}$, challenging the development of a
tractable solution that jointly estimates ${\mathbf{h}}$ and ${\mathbf{S}}$.
Our formulation bypasses this problem by recasting the optimization variables
as ${\mathbf{H}}$ and ${\mathbf{S}}$, leading to the (more tractable) bilinear
constraint in (IV). Nonetheless, if estimating ${\mathbf{h}}$ is the ultimate
goal, this can be readily achieved from ${\hat{\mathbf{H}}}$ and
${\hat{\mathbf{S}}}$ as

${\hat{\mathbf{h}}}=\Big{(}{\mathrm{vec}}({\mathbf{I}}),{\mathrm{vec}}({\hat{\mathbf{S}}}),...,{\mathrm{vec}}({\hat{\mathbf{S}}}^{N-1})\Big{)}^{\dagger}{\mathrm{vec}}({\hat{\mathbf{H}}}).$
(11)

The approach put forth in (IV) has two main advantages. First, while most
works formulate the recovery of the GF in the spectral domain, our formulation
operates in the vertex domain. Working on the spectral domain would imply
finding the Vandermonde GFT matrix $\bm{\Psi}$. Since this matrix involves
high-order polynomials of the eigenvalues of the GSO, it is also prone to
numerical instability and error accumulation [7]. Even if approaches that
bypass this issue by estimating the graph-frequency response
${\tilde{\mathbf{h}}}=\bm{\Psi}{\mathbf{h}}$ in lieu of ${\mathbf{h}}$ are
adopted, the estimation would still be challenging since they require
computing the eigenvectors ${\mathbf{V}}$, which are known to be highly
sensitive to errors in the GSO (especially those associated with small
eigenvalues)[24, 23]. On top of this, characterizing the spectral errors and
incorporating those to the optimization is not a trivial task. The second
advantage emanates from casting the true GSO ${\mathbf{S}}$ as an explicit
optimization variable. As already explained, this approach is robust to error
accumulation and facilitates the incorporation of the (additive) effect of the
perturbations into the optimization. An additional benefit is that we obtain a
denoised version (enhanced estimation) of the true GSO, which can be
practically relevant in most real-world applications.

In a nutshell, in the context of robust GF identification, choosing a
formulation that: i) works entirely in the vertex domain, ii) considers
${\mathbf{S}}$ as an explicit optimization variable, and iii) codifies the GF
structure via the constraint
${\mathbf{H}}{\mathbf{S}}={\mathbf{S}}{\mathbf{H}}$, exhibits multiple
advantages. However, it must be noted that the number of optimization
variables is larger than in classical approaches (adding computational
complexity) and that the bilinear filtering constraint
${\mathbf{H}}{\mathbf{S}}={\mathbf{S}}{\mathbf{H}}$, while more tractable than
its polynomial counterpart, is still non-convex. Alternatives to deal with
these issues are discussed in later sections.

###  IV-A Alternating minimization for robust GF identification

This section presents a systematic efficient approach to find an approximate
solution to (IV). Since the goal is to design specific algorithms, from this
section onwards, we particularize the GSO distance to
$d({\mathbf{S}},{\bar{\mathbf{S}}})=\|{\mathbf{S}}-{\bar{\mathbf{S}}}\|_{0}$,
so that, according to the discussion in Sec. III-A, the focus is on graph
perturbations that create and destroy links. Apart from its practical
relevance, the reason for choosing the $\ell_{0}$ norm as a distance is also
motivated by its more intricate (challenging) structure. Indeed, the
algorithms presented next can be easily adapted to (more tractable) distances
associated with alternative perturbation models. Having clarified this, the
main obstacle to solving (IV) is its lack of convexity, which emanates from
two different _sources_ : _(s1)_ the $\ell_{0}$ norms in the objective, and
_(s2)_ the bilinear constraint involving ${\mathbf{S}}$ and ${\mathbf{H}}$.
Next, we explain the strategy adopted to deal with them and find a solution to
(IV) by solving a succession of convex problems.

  * •

Regarding the $\ell_{0}$ norm in _(s1)_ , a workhorse approach is to replace
it with its convex surrogate, the $\ell_{1}$ norm. However, it is possible to
exploit more sophisticated (non-convex) alternatives that typically lead to
sparser solutions. The one chosen in this paper is to approximate the
$\ell_{0}$ norm of a generic matrix ${\mathbf{Z}}\in{\mathbb{R}}^{I\times J}$
using the logarithmic penalty

$\|{\mathbf{Z}}\|_{0}\approx
r_{\delta}({\mathbf{Z}}):=\sum_{i=1}^{I}\sum_{j=1}^{J}\log(|Z_{ij}|+\delta),$
(12)

where $\delta$ is a small positive constant [36]. The non-convexity of the
logarithm can be handled efficiently by relying on a majorization-minimization
approach (MM) [37], which considers an iterative linear approximation leading
to an iterative re-weighted $\ell_{1}$ norm. It is worth noting that, since we
will consider an iterative algorithm to deal with the bilinearity of (IV), the
iterative nature of the re-weighted $\ell_{1}$ norm will not impose a
significant computational burden. Details on the exact form of this sparse
regularizer will be provided soon, when describing the estimation of
${\mathbf{S}}$.

  * •

To deal with the bilinear terms in _(s2)_ , we adopt an alternating
optimization approach [38] resulting in an iterative algorithm where the
optimization variables ${\mathbf{H}}$ and ${\mathbf{S}}$ are updated in _two
separate iterative steps_. At each step, we optimize over one of the
optimization variables with the other remaining fixed, resulting in two
simpler problems that can be solved efficiently. The details about the
specific steps will be provided shortly.

Taking into account these considerations, the first task to implement our
approach is to rewrite the problem in (IV) as

$\min_{{\mathbf{S}}\in{\mathcal{S}},{\mathbf{H}}}\|{\mathbf{Y}}\\!-\\!{\mathbf{H}}{\mathbf{X}}\|_{F}^{2}\\!+\\!\lambda
r_{\delta_{1}}\\!({\mathbf{S}}\\!-\\!{\bar{\mathbf{S}}})\\!+\\!\beta
r_{\delta_{2}}\\!({\mathbf{S}})\\!+\\!\gamma\|{\mathbf{S}}{\mathbf{H}}\\!-\\!{\mathbf{H}}{\mathbf{S}}\|_{F}^{2},$
(13)

where we recall that $r_{\delta}(\cdot)$ was introduced in (12). Note that: i)
the logarithmic penalty has also been used to promote sparsity in the term
${\mathbf{S}}-{\bar{\mathbf{S}}}$ since we selected the $\ell_{0}$ norm as the
distance between ${\mathbf{S}}$ and ${\bar{\mathbf{S}}}$, and ii) the
constraint ${\mathbf{S}}{\mathbf{H}}={\mathbf{H}}{\mathbf{S}}$ was relaxed and
rewritten as a regularizer, a formulation more amenable to an alternating
optimization approach.

The next task is to solve (13) by means of an iterative algorithm that blends
techniques from alternating optimization and MM approaches. Specifically, for
a maximum of $t_{max}$ iterations, we run the following two steps at each
iteration $t=0,...,t_{max}-1$.

Step 1: GF Identification. We estimate the block of $N^{2}$ variables
collected in ${\mathbf{H}}$ while the current estimate of the GSO, denoted as
${\mathbf{S}}^{(t)}$, remains fixed. This results in the convex optimization
problem

$\displaystyle{\mathbf{H}}^{(t+1)}=\arg\min_{{\mathbf{H}}}$
$\displaystyle\|{\mathbf{Y}}\\!-\\!{\mathbf{H}}{\mathbf{X}}\|_{F}^{2}\\!+\\!\gamma\|{\mathbf{S}}^{(t)}{\mathbf{H}}\\!-\\!{\mathbf{H}}{\mathbf{S}}^{(t)}\|_{F}^{2},$
(14)

an LS minimization whose closed-form solution is

$\displaystyle{\mathrm{vec}}({\mathbf{H}}^{(t+1)})\\!=$
$\displaystyle\big{(}{\mathbf{X}}{\mathbf{X}}^{\top}\\!\\!\otimes\\!{\mathbf{I}}\\!+\\!\gamma({\mathbf{S}}{\mathbf{S}}^{\top}\\!\\!\oplus\\!{\mathbf{S}}^{\top}{\mathbf{S}}\\!-\\!{\mathbf{S}}^{\top}\\!\\!\otimes\\!{\mathbf{S}}^{\top}\\!\\!-\\!{\mathbf{S}}\\!\otimes\\!{\mathbf{S}})\big{)}^{-1}$
$\displaystyle\times({\mathbf{X}}\\!\otimes{\mathbf{I}}){\mathrm{vec}}({\mathbf{Y}}).$
(15)

Here, $\otimes$ is the Kronecker product, $\oplus$ is the Kronecker sum, and
${\mathbf{I}}$ is the identity matrix of size $N\times N$. Also note that
(IV-A) omitted the iteration superscript in ${\mathbf{S}}^{(t)}$ to alleviate
notation.

Step 2: Graph Denoising. Following an MM scheme, we optimize an upper bound of
(13) where the logarithmic penalties are linearized. Then, we estimate the
block of $N^{2}$ variables collected in ${\mathbf{S}}$ while the current
estimate of the GF ${\mathbf{H}}^{(t+1)}$ remains fixed. This yields

$\displaystyle{\mathbf{S}}^{(t+1)}=\arg\min_{{\mathbf{S}}\in{\mathcal{S}}}$
$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}\big{(}\lambda\bar{\Omega}_{ij}^{(t)}|S_{ij}-\bar{S}_{ij}|+\beta\Omega_{ij}^{(t)}|S_{ij}|\big{)}$
$\displaystyle+\gamma\|{\mathbf{S}}{\mathbf{H}}^{(t+1)}-{\mathbf{H}}^{(t+1)}{\mathbf{S}}\|_{F}^{2},$
(16)

where $\bar{\bm{\Omega}}^{(t)}$ and $\bm{\Omega}^{(t)}$ are computed in an
entry-wise fashion based on the GSO estimate from the previous iteration as

$\displaystyle\bar{\Omega}_{ij}^{(t)}=\frac{1}{|S_{ij}^{(t)}-\bar{S}_{ij}|+\delta_{1}},$
$\displaystyle\Omega_{ij}^{(t)}=\frac{1}{|S_{ij}^{(t)}|+\delta_{2}}.$ (17)

Input: ${\mathbf{X}}$, ${\mathbf{Y}}$, ${\bar{\mathbf{S}}}$

Output : ${\hat{\mathbf{H}}}$, ${\hat{\mathbf{S}}}$.

1 Initialize ${\mathbf{S}}^{(0)}$ as ${\mathbf{S}}^{(0)}={\bar{\mathbf{S}}}$.

2 for  _$t=0$ to $t_{max}-1$_ do

3 Compute ${\mathbf{H}}^{(t+1)}$ by solving (IV-A) fixing
${\mathbf{S}}^{(t)}$.

4 Update $\bm{\Omega}^{(t)}$ and $\bar{\bm{\Omega}}^{(t)}$ as in (17).

5 Compute ${\mathbf{S}}^{(t+1)}$ by solving (IV-A) using
${\mathbf{H}}^{(t+1)}$, $\bm{\Omega}^{(t)}$, and $\bar{\bm{\Omega}}^{(t)}$.

6

7 end for

${\hat{\mathbf{H}}}={\mathbf{H}}^{(t_{max})},\;{\hat{\mathbf{S}}}={\mathbf{S}}^{(t_{max})}$.

Algorithm 1 Robust GF identification with graph denoising.

The overall alternating algorithm is summarized in Alg. 1, where a fixed
number of iterations is considered. The algorithm starts by initializing the
GSO as ${\mathbf{S}}^{(0)}={\bar{\mathbf{S}}}$ (although other options could
also be appropriate), and then, it iterates between Steps 1 and 2 for a fixed
number of epochs (or until some stopping criterion is met). In this regard, a
key feature of the algorithm is that it is guaranteed to converge to a
stationary point, as is formally stated next.

######  Theorem 1.

Denote as $f({\mathbf{H}},{\mathbf{S}})$ the objective function in (13), and
let ${\mathcal{Z}}^{*}$ be the set of stationary points of $f$. Let
${\mathbf{z}}^{(t)}=[{\mathrm{vec}}({\mathbf{H}}^{\\!(t)})^{\\!\top}\\!,{\mathrm{vec}}({\mathbf{S}}^{\\!(t)})^{\\!\top}\\!]^{\\!\top}$
represent the solution provided by the iterative algorithm (IV-A)-(IV-A) after
$t$ iterations. Assuming that i) the GSO does not have repeated eigenvalues
and ii) every row of
${\tilde{\mathbf{X}}}\\!=\\!{\mathbf{V}}^{\\!-1}{\mathbf{X}}$ has at least one
nonzero entry, then ${\mathbf{z}}^{\\!(t)}$ converges to a stationary point of
$f$ as $t$ goes to infinity, i.e.,

$\lim_{t\to\infty}\mathsf{d}({\mathbf{z}}^{(t)}~{}|{\mathcal{Z}}^{*})=0,$

with
$\mathsf{d}({\mathbf{z}}~{}|{\mathcal{Z}}^{*}):=\min_{{\mathbf{z}}^{*}\in{\mathcal{Z}}^{*}}\|{\mathbf{z}}-{\mathbf{z}}^{*}\|_{2}$.

The proof relies on the convergence results shown in [39, Th. 1b] and the
details are provided in App. A. Note that the convergence of the algorithm was
not self-evident since the original optimization problem in (13) is non-convex
and Step 2 is minimizing an upper-bound of the original objective function.
The sufficient conditions in i) and ii) guarantee that every graph frequency
is excited so that the GF is identifiable and (14) has a unique solution,
which is a requirement for convergence (see Prop. 1 in App. A for details).
Clearly, condition ii) is fulfilled even for $M=1$ if all the entries of the
vector ${\tilde{\mathbf{x}}}={\mathbf{V}}^{-1}{\mathbf{x}}$ are nonzero.
Alternatively, when $M>1$ and ii) is satisfied, condition i) can be relaxed.

Another relevant element in the proposed algorithm is the weight $\gamma$. If
$\gamma$ is set to a value that is too large, the GF estimated in the first
iteration ${\mathbf{H}}^{(1)}$ will be an (almost exact) polynomial of
$\bar{{\mathbf{S}}}$ so that the algorithm will converge quickly to the same
solution as that of the non-robust design [cf. (IV) with
${\mathbf{S}}=\bar{{\mathbf{S}}}$]. On the other hand, if $\gamma$ is too
close to zero the two problems decouple and the solution converges quickly to
that of the two separated problems [cf. (14) and (IV-A) with $\gamma=0$]. As a
result, the value of the parameter must be chosen carefully. In this context,
schemes that start with a small $\gamma$ to encourage the exploration during
the warm-up phase, and then increase $\gamma$ as the iteration index grows to
guarantee that the final ${\hat{\mathbf{H}}}$ is a polynomial of
${\hat{\mathbf{S}}}$ are a suitable alternative for the setup at hand.

Finally, one drawback of the proposed robust GF identification algorithm is
that the optimization problems in (14) and (IV-A) may be slow when dealing
with large graphs. However, we will mitigate this issue by introducing an
efficient implementation that reduces the computational complexity of the
overall algorithm (see Sec. VI).

###  IV-B Leveraging stationary observations

The alternating convex approximation in Alg. 1 exploits the fact that
${\mathbf{X}}$ and ${\mathbf{Y}}$ are linearly related via ${\mathbf{H}}$,
which is a polynomial of ${\mathbf{S}}$. However, in setups where the
perturbations in ${\bar{\mathbf{S}}}$ are very large, obtaining accurate
estimates of ${\mathbf{S}}$ and ${\mathbf{h}}$ from ${\hat{\mathbf{H}}}$ may
still be challenging. One alternative to overcome this issue is to leverage
the additional structure potentially present in our data. Indeed, as detailed
in the introduction, it is common to consider setups where the signals exhibit
additional properties depending on the supporting graph, with notable examples
including graph-bandlimited signals [5, 6], diffused sparse graph signals [14,
15], or graph stationary signals [33, 40, 22]. Clearly, incorporating such
additional information into the optimization problem would enhance its
estimation performance.

This section explores this path, restricting our attention to the case where
the observed signals are stationary on ${\mathcal{G}}$. The motivation for
this decision is that, due to the tight connection between graph-stationary
signals and GFs (see Sec. II), the formulation in (13) and Alg. 1 require
relatively minor modifications to incorporate the assumption of ${\mathbf{X}}$
and ${\mathbf{Y}}$ being stationary on ${\mathbf{S}}$, leaving the
incorporation of additional signal models as future work. To formulate the
updated problem, recall that the covariance matrix of a stationary graph
signal can be expressed as a polynomial of the GSO (see Sec. II). Therefore,
incorporating stationarity calls for modifying (13) as

$\displaystyle\\!\min_{{\mathbf{S}}\in{\mathcal{S}},{\mathbf{H}}}$
$\displaystyle\;\\!\|{\mathbf{Y}}\\!-\\!{\mathbf{H}}{\mathbf{X}}\|_{F}^{2}+\lambda
r_{\delta_{1}}\\!({\mathbf{S}}\\!-\\!{\bar{\mathbf{S}}})+\beta
r_{\delta_{2}}\\!({\mathbf{S}})+\gamma\|{\mathbf{S}}{\mathbf{H}}\\!-\\!{\mathbf{H}}{\mathbf{S}}\|_{F}^{2}$
$\displaystyle\\!\mathrm{\;\;s.\;to:}$
$\displaystyle\;\|{\mathbf{C}}_{\mathbf{y}}{\mathbf{S}}\\!-\\!{\mathbf{S}}{\mathbf{C}}_{\mathbf{y}}\\!\|_{F}^{2}\\!\leq\\!\epsilon_{\mathbf{y}},\,\\!\|{\mathbf{C}}_{\mathbf{x}}{\mathbf{S}}\\!-\\!{\mathbf{S}}{\mathbf{C}}_{\mathbf{x}}\\!\|_{F}^{2}\\!\leq\\!\epsilon_{\mathbf{x}},$
(18)

where ${\mathbf{C}}_{\mathbf{y}}$ and ${\mathbf{C}}_{\mathbf{x}}$ denote the
covariance matrices of ${\mathbf{Y}}$ and ${\mathbf{X}}$, respectively. If the
covariances are perfectly known, then the corresponding parameters
$\epsilon_{\mathbf{y}}$ and $\epsilon_{\mathbf{x}}$ are set to zero.
Alternatively, if the ${\mathbf{C}}_{\mathbf{y}}$ and
${\mathbf{C}}_{\mathbf{x}}$ are the sample estimates of the true covariances,
then the values of $\epsilon_{\mathbf{y}}$ and $\epsilon_{\mathbf{x}}$ must be
selected based on the quality of the estimators (accounting, e.g., for the
number of available observations $M$).

The constraints in (IV-B) capture the graph-stationarity assumption by
promoting the commutativity with the true GSO. Therefore, such constraints are
considered in the graph denoising step [cf. (IV-A)]. In addition, since
${\mathbf{C}}_{\mathbf{y}}$, ${\mathbf{C}}_{\mathbf{x}}$ and ${\mathbf{H}}$
are all polynomials of ${\mathbf{S}}$, the equalities
${\mathbf{C}}_{\mathbf{y}}{\mathbf{H}}={\mathbf{H}}{\mathbf{C}}_{\mathbf{y}}$
and
${\mathbf{C}}_{\mathbf{x}}{\mathbf{H}}={\mathbf{H}}{\mathbf{C}}_{\mathbf{x}}$
must hold as well, so it is also possible to augment the GF identification
step [cf. (14)] with the corresponding constraints. While in the interest of
brevity, we do not spell out all the possible formulations here, the impact of
several of these alternatives is numerically analyzed in Sec. VII. Finally, it
is important to note that, since the stationarity constraints are quadratic
and convex, the convergence described in Theorem 1 also holds true for the
iterative algorithm associated with (IV-B).

##  V Joint robust identification of multiple GFs

In Sec. IV, we approached the problem of identifying a single GF
${\mathbf{H}}$ defined over a single graph ${\mathcal{G}}$. However, in a
variety of situations we encounter multiple processes (signals) over the same
graph ${\mathcal{G}}$. Consider for example a network of weather stations
measuring the temperature, humidity, and wind speed. Each of these
measurements corresponds to observations of a different process, all of them
taking place over a common graph. Intuitively, since all the GFs are related
by the underlying graph ${\mathcal{G}}$, we propose a _joint_ GF
identification approach that exploits this relationship to enhance the quality
of the estimation. We focus first on the case where the input-output signals
associated with each GF (graph process) are observed separately. Later in the
section, we address a slightly more involved case where the GFs model the (AR)
dynamics of a time-varying graph signal and, as a result, the observed signals
are intertwined.

Consider a set of $K$ unknown GFs $\\{{\mathbf{H}}_{k}\\}_{k=1}^{K}$, all
represented by $N\times N$ matrices and defined over the graph
${\mathcal{G}}$. To be consistent with 1, we assume that: i) the true
${\mathbf{S}}$ is unknown and only the perturbed version ${\bar{\mathbf{S}}}$
is available; ii) all ${\mathbf{H}}_{k}$ are polynomials of the _same_ GSO
${\mathbf{S}}$; and iii) for each $k$, matrices
${\mathbf{X}}_{k}\in{\mathbb{R}}^{N\times M_{k}}$ and
${\mathbf{Y}}_{k}\in{\mathbb{R}}^{N\times M_{k}}$ collect the observed input
and output graph signals and are related via

${\mathbf{Y}}_{k}={\mathbf{H}}_{k}{\mathbf{X}}_{k}+{\mathbf{W}}_{k},$ (19)

with ${\mathbf{H}}_{k}=\sum_{r=0}^{N-1}h_{r,k}{\mathbf{S}}^{r}$ and
${\mathbf{W}}_{k}$ being a white random matrix capturing observation noise and
model inaccuracies. Then, we aim at estimating the GFs
$\\{{\mathbf{H}}_{k}\\}_{k=1}^{K}$ in a joint fashion while taking into
account the inaccuracies in the topology of ${\mathcal{G}}$. This is
summarized in the following problem statement.

######  Problem 2.

Let ${\mathcal{G}}$ be a graph with $N$ nodes, let
${\mathbf{S}}\in{\mathbb{R}}^{N\times N}$ be the true (unknown) GSO associated
with ${\mathcal{G}}$, and let ${\bar{\mathbf{S}}}\in{\mathbb{R}}^{N\times N}$
be the perturbed (observed) GSO. Moreover, let
${\mathbf{X}}_{k}\in{\mathbb{R}}^{N\times M_{k}}$ and
${\mathbf{Y}}_{k}\in{\mathbb{R}}^{N\times M_{k}}$ be the matrices collecting
the $M_{k}$ observed input and output graph signals associated with
$k=1,...,K$ network processes, all defined over ${\mathcal{G}}$ and adhering
to the model in (19). Our goal is to use
$\\{{\mathbf{X}}_{k}\\!\\}_{k=1}^{K}$, $\\{{\mathbf{Y}}_{k}\\!\\}_{k=1}^{K}$,
and ${\bar{\mathbf{S}}}$ to learn the $K$ GFs
$\\{{\mathbf{H}}_{k}\\}_{k=1}^{K}$ that best fit the data, along with an
enhanced estimation of ${\mathbf{S}}$. To that end, we make the following
assumptions:
(AS2) ${\mathbf{S}}$ and ${\bar{\mathbf{S}}}$ are close according to some
metric $d({\mathbf{S}},{\bar{\mathbf{S}}})$, i.e., the observed perturbations
are “small” in some sense.
(AS3) Every ${\mathbf{H}}_{k}$ is a polynomial of ${\mathbf{S}}$.

Assumption (AS2), which was also considered in 1, promotes the tractability of
the problem by ensuring that ${\mathbf{S}}$ and ${\bar{\mathbf{S}}}$ are
related. As discussed in Sec. III-A, the distance function $d(\cdot,\cdot)$
must be selected depending on the perturbation model at hand. (AS3) captures
the key fact that all the matrices ${\mathbf{H}}_{k}$ are GFs of the _same_
GSO, establishing a link that can be leveraged via a joint estimation
(optimization) of the $K$ GFs. Implementing an approach similar to that in
Sec. IV (i.e., working on the vertex domain, considering the true GSO as an
explicit optimization variable, accounting for the GF structure via a
commutativity constraint, and assuming that the graph perturbations create and
destroy links), the multi-filter counterpart to (13) that codifies 2 is

$\displaystyle\\!\min_{{\mathbf{S}}\in{\mathcal{S}},\\{{\mathbf{H}}_{k}\\}_{k=1}^{K}}$
$\displaystyle\sum_{k=1}^{K}\alpha_{k}\|{\mathbf{Y}}_{k}-\\!{\mathbf{H}}_{k}{\mathbf{X}}_{k}\|_{F}^{2}\\!+\\!\lambda
r_{\delta_{1}}({\mathbf{S}}-{\bar{\mathbf{S}}})$ $\displaystyle+\beta
r_{\delta_{2}}({\mathbf{S}})+\sum_{k=1}^{K}\gamma\|{\mathbf{S}}{\mathbf{H}}_{k}\\!\\!-\\!{\mathbf{H}}_{k}{\mathbf{S}}\|_{F}^{2}.$
(20)

Ideally, the value of the positive weight $\alpha_{k}$ must be selected based
on the norm of ${\mathbf{W}}_{k}$ (e.g., prior information on the noise level
and the number of signal pairs $M_{k}$). If none is available, then
$\alpha_{k}=1$ for all $k$. Equally important, the fact of pursuing a joint
optimization implies that each ${\mathbf{H}}_{k}$ contributes with a
regularization term
$\|{\mathbf{S}}{\mathbf{H}}_{k}-{\mathbf{H}}_{k}{\mathbf{S}}\|_{F}^{2}$
promoting the commutativity of the $k$-th GF with the _single_ ${\mathbf{S}}$.
Intuitively, having the same ${\mathbf{S}}$ in all these terms couples the
optimization across $k$ and contributes to reduce the uncertainty over
${\mathbf{S}}$, leading to enhanced estimates of both ${\mathbf{S}}$ and
$\\{{\mathbf{H}}_{k}\\}_{k=1}^{K}$. As a result, the joint GF identification
approach is expected to provide better results than estimating each
${\mathbf{H}}_{k}$ separately by solving $K$ instances of (13). We validate
this hypothesis numerically via the experiments in Sec. VII.

Following a motivation similar to that in the previous section, we deal with
the non-convex minimization in (V) designing an alternating optimization
algorithm that breaks the bilinear terms ${\mathbf{S}}{\mathbf{H}}_{k}$ and
${\mathbf{H}}_{k}{\mathbf{S}}$, and approximates the logarithmic terms with a
linear upper-bound. The resulting algorithm solves iteratively the following
two subproblems for $t=1,...,t_{\max}$ iterations.

Step 1: Multiple GF Identification. Given the current estimate
${\mathbf{S}}^{(t)}$, we solve the optimization problem in (V) with respect to
each ${\mathbf{H}}^{(k)}$. This yields

$\displaystyle\\!\\!{\mathbf{H}}_{k}^{(t+1)\\!}\\!\\!=\\!\operatornamewithlimits{argmin}_{{\mathbf{H}}_{k}}\alpha_{k}\|\\!{\mathbf{Y}}_{k}\\!\\!-\\!\\!{\mathbf{H}}_{k}{\mathbf{X}}_{k}\\!\|_{\\!F}^{\\!2}\\!\\!+\\!\gamma\|\\!{\mathbf{S}}^{\\!(t)\\!}{\mathbf{H}}_{k}\\!\\!-\\!\\!{\mathbf{H}}_{k}{\mathbf{S}}^{\\!(t)}\\!\|_{\\!F}^{\\!2},\\!$
(21)

whose closed-form solution can be found using (IV-A) replacing $\gamma$ with
$\gamma/\alpha_{k}$, ${\mathbf{X}}$ with ${\mathbf{X}}_{k}$, and
${\mathbf{Y}}$ with ${\mathbf{Y}}_{k}$. Note that since the only coupling
across GFs is via the GSO, (21) estimates each ${\mathbf{H}}_{k}^{(t+1)}$
separately from the other GFs, solving $K$ LS problems (each with $N^{2}$
unknowns). Furthermore, if multiple processors are available, (21) can be run
in parallel across $k$.

Step 2: Graph Denoising. Given the current estimates of the GFs
$\\{{\mathbf{H}}_{k}^{(t+1)}\\}_{k=1}^{K}$, we follow an MM scheme that,
minimizing a linear upper-bound of the logarithmic penalties, yields the
estimate of the GSO via

$\displaystyle{\mathbf{S}}^{(t+1)}=\operatornamewithlimits{argmin}_{{\mathbf{S}}\in{\mathcal{S}}}$
$\displaystyle\sum_{ij=1}^{N}\big{(}\lambda\bar{\Omega}_{ij}^{(t)}|S_{ij}-\bar{S}_{ij}|+\beta\Omega_{ij}^{(t)}|S_{ij}|\big{)}$
$\displaystyle+\sum_{k=1}^{K}\gamma\|{\mathbf{S}}{\mathbf{H}}_{k}^{(t+1)}-{\mathbf{H}}_{k}^{(t+1)}{\mathbf{S}}\|_{F}^{2},$
(22)

where $\bm{\Omega}$ and $\bar{\bm{\Omega}}$ are obtained as in (17).

The solution to 2 is simply given by
${\hat{\mathbf{S}}}={\mathbf{S}}^{(t_{max})}$ and
${\hat{\mathbf{H}}}_{k}={\mathbf{H}}_{k}^{(t_{max})}$ for every $k$. Similar
to (13), convergence to a stationary point of (V) is guaranteed, as formally
stated next.

######  Corollary 1.

Denote as $f(\\{{\mathbf{H}}_{k}\\}_{k=1}^{K},{\mathbf{S}})$ the objective
function in (V). If
${\mathbf{z}}^{(t)}=[{\mathrm{vec}}({\mathbf{H}}_{1}^{(t)})^{\top},...,{\mathrm{vec}}({\mathbf{H}}_{K}^{(t)})^{\top},{\mathrm{vec}}({\mathbf{S}})^{\top}]^{\top}$
represents the solution provided by the iterative algorithm (21)-(V) after $t$
iterations and every ${\mathbf{X}}_{k}$ excites all graph frequencies, then
${\mathbf{z}}^{(t)}$ converges to a stationary point of $f$ as the number of
iterations $t$ goes to infinity.

The key to prove Theorem 1, which established the convergence to a stationary
point for the robust estimation of a single GF, was to show that the
optimization problem in (13) and the proposed algorithm satisfied the
conditions in [39, Th. 1b]. The formulation we put forth for the multi-filter
case resembles closely that of the single-filter case, and, as a result, it is
not difficult to show that those conditions also hold true for the problem in
(V) (see App. A for details).

The discussion and formulations in Sec. IV-B dealing with incorporating
additional information about the input-output signals into the optimization
are also pertinent for the setup in this section. The details of such a
formulation are omitted for brevity, but it will be explored in the
experimental section.

###  V-A Joint GF identification for time series

A slightly different, practically relevant, setup where multiple GFs need to
be estimated is that of graph-based multivariate time series. In that setup,
each variable is associated with a node of the graph and the multiple graph-
signal observations correspond to different instants of a time-varying graph
signal. AR and moving-average (MA) modeling of time series has a long
tradition, with common approaches to decrease the degrees of freedom including
limiting the memory of the series and assuming that matrices of coefficients
relating different time instants are low rank [41]. In the context of graph
signals and network processes, a natural approach is to constrain the matrices
of coefficients to be GFs, all defined over the same graph [42, 43]. This
section introduces a variation of the problem in (V) tailored to this setup.

To introduce the multiple-graph identification problem formally, let
${\mathbf{X}}_{\kappa}$ and ${\mathbf{Y}}_{\kappa}$ denote a collection of
$M_{\kappa}$ graph signals corresponding to measurements of a network process
for $\kappa=1,...,\kappa_{max}$ time instants. Suppose now that
${\mathbf{Y}}_{\kappa}$ can be accurately modeled by an AR dynamics with
memory $K$ so, at every instant $\kappa$, the observations
${\mathbf{Y}}_{\kappa}$ satisfy the equation

${\mathbf{Y}}_{\kappa}=\sum_{k=1}^{K}{\mathbf{H}}_{k}{\mathbf{Y}}_{\kappa-k}+{\mathbf{X}}_{\kappa},\;\mathrm{with}\;{\mathbf{H}}_{k}=\sum_{r=0}^{N-1}h_{r,k}{\mathbf{S}}^{r},$
(23)

where ${\mathbf{X}}_{\kappa}$ is the exogenous input, and the GF
${\mathbf{H}}_{k}$ models the influence that the signal observations from the
time instant $\kappa-k$ exert on the (current) signal at time $\kappa$.

Suppose now that: i) we have access to an estimated (imperfect) graph
${\bar{\mathbf{S}}}$, ii) the value of the graph signals at different time
instants is available, and iii) our goal is to estimate the set of matrices
(GFs) $\\{{\mathbf{H}}_{k}\\}_{k=1}^{K}$ in (23) that describe the dynamics of
the multivariate time series. This can be accomplished as

$\displaystyle\\!\min_{{\mathbf{S}}\in{\mathcal{S}},\\{{\mathbf{H}}_{k}\\}_{k=1}^{K}}\sum_{\kappa=K+1}^{\kappa_{max}}\Big{\|}{\mathbf{Y}}_{\kappa}-{\mathbf{X}}_{\kappa}-\\!\sum_{k=1}^{K}{\mathbf{H}}_{k}{\mathbf{Y}}_{\kappa-k}\Big{\|}_{F}^{2}$
$\displaystyle\\!\hskip 17.07182pt+\\!\lambda
r_{\delta_{1}}({\mathbf{S}}-{\bar{\mathbf{S}}})\\!+\\!\beta
r_{\delta_{2}}({\mathbf{S}})\\!+\\!\sum_{k=1}^{K}\gamma\|{\mathbf{S}}{\mathbf{H}}_{k}\\!\\!-\\!{\mathbf{H}}_{k}{\mathbf{S}}\|_{F}^{2}.\\!$
(24)

The main difference relative to (V) is in the first term, which accounts for
the new observation model [cf. (19) vs. (23)]. Note that we assume that the
exogenous input ${\mathbf{X}}_{\kappa}$ is observed. If that were not the
case, it would suffice to remove ${\mathbf{X}}_{\kappa}$ from the objective
(possibly updating the Frobenius norm in case statistical knowledge about
${\mathbf{X}}_{\kappa}$ were available). Albeit the differences, the problem
in (V-A) is closely related to (V), with the sources of non-convexities being
the same. As a result, we approach its solution with a modified version of
Alg. 1 which, at each iteration $t$, runs two steps. In the first one, we
estimate each of the $K$ GFs by solving

$\displaystyle{\mathbf{H}}_{k}^{(t+1)}\\!$
$\displaystyle\\!=\\!\operatornamewithlimits{argmin}_{{\mathbf{H}}_{k}}\\!\\!\\!\sum_{\kappa=K+1}^{\kappa_{max}}\\!\\!\Big{\|}{\mathbf{Y}}_{\kappa}\\!-\\!{\mathbf{X}}_{\kappa}\\!-\\!{\mathbf{H}}_{k}{\mathbf{Y}}_{\kappa-k}\\!-\\!\\!\\!\sum_{k^{\prime}<k}\\!\\!{\mathbf{H}}_{k^{\prime}}^{(t+1)}{\mathbf{Y}}_{\kappa-k^{\prime}}$
$\displaystyle\\!\\!-\\!\\!\sum_{k^{\prime}<k}\\!\\!{\mathbf{H}}_{k^{\prime}}^{(t)}{\mathbf{Y}}_{\kappa-k^{\prime}}\Big{\|}_{\\!F}^{\\!2}+\\!\\!\sum_{k=1}^{K}\\!\gamma\Big{\|}{\mathbf{S}}^{(t)}{\mathbf{H}}_{k}-{\mathbf{H}}_{k}{\mathbf{S}}^{(t)}\Big{\|}_{\\!F}^{\\!2},\\!\\!$
(25)

which is different from the previous GF identification step [cf. (21)]. In
contrast, the graph-denoising step in (V) remains the same. Note that (V-A)
updates each GF separately in a cyclic way by solving an LS problem with
$N^{2}$ unknowns. Alternative implementations include using
${\mathbf{H}}_{k^{\prime}}^{(t)}$ in lieu of
${\mathbf{H}}_{k^{\prime}}^{(t+1)}$ for all $k^{\prime}<k$ (so that a parallel
implementation is enabled) as well as considering a single LS problem with
$KN^{2}$ unknowns.

Finally, it is worth emphasizing that the formulation introduced in this
section can be used as a starting point to design more general robust schemes
for multivariate time series defined over a graph. Dealing with both AR and MA
matrices, assuming that the memory of the system is not known, having only
partial/statistical information on the exogenous input, and observing the
signals at only a subset of nodes are all examples of setups of interest.
Since our goal in this section was to demonstrate the relevance of a robust
multiple GF formulation in the context of multivariate time series, to
facilitate exposition we restricted our discussion to the relatively simple
case in (23), but many other setups (including those previously listed) will
be subject of our future work.

##  VI Efficient implementation of the robust GF identification algorithm

The algorithms proposed up to this point are able to find a solution to the
robust GF identification problem in polynomial time. However, their
computational complexity scales with the number of nodes as $N^{7}$. To
facilitate the deployment in setups where $N$ is large, this section puts
forth an efficient implementation that reduces the number of operations.

The new algorithm (summarized in Alg. 2 ) preserves the core structure of Alg.
1, with an outer loop that, at each iteration, runs two steps: one involving
the estimation of the GF(s) and another one dealing with the estimation of the
GSO. The main difference is that now, instead of finding the exact solution to
those two problems, we obtain an approximate solution. While the details,
which are step-dependent, will be specified in the next paragraphs, the
overall idea is that for each of the steps we run a few simple
(gradient/proximal) iterations. Although Alg. 2 involves two nested loops, the
complexity of the problems in the inner loop is cut down significantly, so
that the overall computational overhead is reduced.

To be specific, we describe next the two steps that, at each iteration of the
outer loop $t\\!=\\!0,...,t_{max}\\!-\\!1$, Alg. 2 runs.

Step 1: Efficient GF Identification. Solving the GF-identification step with
the closed-form solution presented in (IV-A) involves inverting a matrix of
size $N^{2}\times N^{2}$, which requires ${\mathcal{O}}(N^{6})$ operations. To
explain our alternative implementation, let
$f_{1}({\mathbf{H}}|{\mathbf{S}}^{(t)})$ denote the objective function in
(14). Since $f_{1}$ is strictly convex and smooth, it can be efficiently
optimized using a gradient descent approach [44].

To that end, for each iteration $t$ of the outer loop, we define the inner
iteration index $\tau$ as well as the sequence of variables
${\mathbf{\check{H}}}^{(\tau)}$ with $\tau=0,...,\tau_{max_{1}}$, which is
initialized as ${\mathbf{\check{H}}}^{(0)}={\mathbf{H}}^{(t)}$. With this
notation at hand, at each iteration $\tau=0,...,\tau_{max_{1}}-1$ of the inner
loop, we update ${\mathbf{\check{H}}}^{(\tau+1)}$ via

${\mathbf{\check{H}}}^{(\tau+1)}={\mathbf{\check{H}}}^{(\tau)}-\mu\nabla
f_{1}({\mathbf{\check{H}}}^{(\tau)}|{\mathbf{S}}^{(t)}).$ (26)

Here, $\mu>0$ is the step size and $\nabla f_{1}$ denotes the gradient of
$f_{1}$ with respect to ${\mathbf{H}}$, which is given by

$\displaystyle\\!\nabla
f_{1}\\!({\mathbf{H}}|{\mathbf{S}}^{(t)})\\!=\\!2\Big{(}{\mathbf{H}}{\mathbf{X}}{\mathbf{X}}^{\top}\\!\\!-\\!\\!{\mathbf{Y}}{\mathbf{X}}^{\top}\Big{)}\\!\\!$
$\displaystyle\hskip
5.69046pt+\\!2\gamma\Big{(}\\!{\mathbf{S}}^{(t)^{\\!\top}}\\!\\!({\mathbf{S}}^{(t)}{\mathbf{H}}\\!-\\!{\mathbf{H}}{\mathbf{S}}^{(t)})\\!-\\!({\mathbf{S}}^{(t)}{\mathbf{H}}\\!-\\!{\mathbf{H}}{\mathbf{S}}^{(t)}){\mathbf{S}}^{(t)^{\top}}\\!\Big{)}.$
(27)

When the $\tau_{max_{1}}$ gradient updates are computed, we conclude the GF-
identification step by setting
${\mathbf{H}}^{(t+1)}={\mathbf{\check{H}}}^{(\tau_{max_{1}})}$.

Since each gradient calculation involves the multiplication of $N\times N$
matrices, the resultant computational complexity is
${\mathcal{O}}(\tau_{max_{1}}N^{3})$, which may go down to
${\mathcal{O}}(\tau_{max_{1}}N^{2.4})$ if an efficient multiplication
algorithm is employed [45]. For large values of $N$, this complexity is
substantially smaller than that required to find the inverse of an
$N^{2}\times N^{2}$ matrix.

Step 2: Efficient graph denoising. Since the optimization in (IV-A) involves
$N^{2}$ variables (the entries in ${\mathbf{S}}$), using an off-the-shelf
convex solver incurs a computational complexity of ${\mathcal{O}}(N^{7})$
[44]. Inspired by the Lasso regression algorithm [46], we optimize
individually over each entry $S_{ij}$ in an iterative manner. The main idea is
running multiple rounds of $N^{2}$ efficient scalar optimizations rather than
dealing with a single but demanding $N^{2}$-dimensional problem. To provide
the details of the scheme developed to estimate ${\mathbf{S}}$, we need to
specify the set of constraints ${\mathcal{S}}$ and introduce some definitions.
Let us focus on the set of adjacency matrices
${\mathcal{S}}_{\mathcal{A}}:=\\{{\mathbf{S}}|S_{ij}\geq 0,\;S_{ii}=0\\}$ and
define the vectors ${\mathbf{s}}:={\mathrm{vec}}({\mathbf{S}})$, vector
${\bar{\mathbf{s}}}:={\mathrm{vec}}({\bar{\mathbf{S}}})$, and the matrix
$\bm{\Sigma}^{(t)}:={\mathbf{H}}^{(t+1)^{\top}}\oplus-{\mathbf{H}}^{(t+1)}$.
With these definitions in place, the minimization in (IV-A) is equivalent to
solving

$\displaystyle\\!\min_{{\mathbf{s}}}$
$\displaystyle\sum_{i=1}^{N^{2}}\left(\lambda\bar{\omega}^{(t)}_{i}|s_{i}-\bar{s}_{i}|+\beta\omega^{(t)}_{i}s_{i}\right)+\gamma\|\bm{\Sigma}^{(t)}{\mathbf{s}}\|_{2}^{2},$
$\displaystyle\\!\mathrm{\;\;s.\;to:}$ $\displaystyle\;\;{\mathbf{s}}\geq
0,\;\;{\mathbf{s}}_{\mathcal{D}}=0,$ (28)

where ${\mathbf{s}}_{\mathcal{D}}$ collects the elements in the diagonal of
${\mathbf{S}}$, and the vectors $\bar{\bm{\omega}}^{(t)}$ and
$\bm{\omega}^{(t)}$ are computed according to (17) but with
${\bar{\mathbf{s}}}^{(t)}$ and ${\mathbf{s}}^{(t)}$ in lieu of
${\bar{\mathbf{S}}}^{(t)}$ and ${\mathbf{S}}^{(t)}$. The constraint
${\mathbf{s}}_{\mathcal{D}}=0$, implies that only the $N^{2}-N$ elements of
${\mathbf{s}}$ representing the off-diagonal entries of ${\mathbf{S}}$ need to
be optimized. The key point to find those $N^{2}-N$ values is to leverage that
the non-differentiable part of the cost in (VI) is separable across $s_{i}$,
postulate $N^{2}-N$ scalar optimization problems (coupled via the $\ell_{2}$
term in the cost), and address the optimization following a projected cyclic
coordinate descent scheme.

To define clearly the operation of Step 2 at each iteration $t$ of the outer
loop, we need to introduce some notation. First, let us denote as $\tau$ the
iteration index for the inner loop, define the sequence of variables
${\mathbf{\check{s}}}^{(\tau)}$ where $\tau=0,...,\tau_{max_{2}}$, and
initialize the sequence as ${\mathbf{\check{s}}}^{(0)}={\mathbf{s}}^{(t)}$.
Moreover, with $\ell\not\in{\mathcal{D}}$ denoting an index of the off-
diagonal elements of the GSO, let $\bm{\sigma}_{\ell}\in{\mathbb{R}}^{N^{2}}$
denote the associated $\ell$-th column of $\bm{\Sigma}^{(t)}$,
$\omega_{\ell}\geq 0$ and $\bar{\omega}_{\ell}\geq 0$ the associated entries
of $\bm{\omega}^{(t)}$ and $\bar{\bm{\omega}}^{(t)}$, and
$\check{s}_{\ell}^{(\tau)}\in{\mathbb{R}}$ the associated entry of
${\mathbf{\check{s}}}^{(\tau)}$ (note that dependence on $t$ was omitted to
facilitate readability). Then, at every iteration
$\tau=0,...,\tau_{max_{2}}-1$ of the inner loop, Alg. 2 optimizes over each
$\check{s}_{\ell}$ separately in a cyclic (successive) way. The advantage of
this approach is that the solution to the _scalar_ optimization over
$\check{s}_{\ell}$ is given in closed form by the following projected soft-
thresholding operation

$\check{s}_{\ell}^{(\tau+1)}=\left\\{\begin{array}[]{cl}\left(-\bar{\lambda}_{\ell}+u^{(\tau)}_{\ell}\right)^{+}&\mathrm{if}\;\bar{s}_{\ell}<-\bar{\lambda}_{\ell}+u^{(\tau)}_{\ell},\\\
\left(\bar{\lambda}_{\ell}+u^{(\tau)}_{\ell}\right)^{+}&\mathrm{if}\;\bar{s}_{\ell}>\bar{\lambda}_{\ell}+u^{(\tau)}_{\ell},\\\
\bar{s}_{\ell}&\mathrm{otherwise},\end{array}\right.$ (29)
$\mathrm{with}\;\;\;\bar{\lambda}_{\ell}=\frac{\lambda\bar{\omega}_{\ell}}{\gamma\bm{\sigma}_{\ell}^{\top}\bm{\sigma}_{\ell}}\;\;\mathrm{and}\;\;\;u^{(\tau)}_{\ell}=\frac{-\beta\omega_{\ell}-\gamma\bm{\sigma}_{\ell}^{\top}{\mathbf{r}}_{\ell}^{(\tau)}}{\gamma\bm{\sigma}_{\ell}^{\top}\bm{\sigma}_{\ell}}.$

Here, $(\cdot)^{+}$ denotes the operation $(x)^{+}=\max(0,x)$, and

${\mathbf{r}}_{\ell}^{(\tau)}:=\sum_{j<\ell}\bm{\sigma}_{j}\check{s}_{j}^{(\tau+1)}+\sum_{j>\ell}\bm{\sigma}_{j}\check{s}_{j}^{(\tau)}.$
(30)

Note that (29) is a soft-thresholding operation with respect to the term
$|s_{i}-\bar{s}_{i}|$. Also, the constraints in ${\mathcal{S}}_{\mathcal{A}}$
are satisfied due to the projection operator
$(\cdot)^{+}\\!:=\\!\max\\{\cdot,\\!0\\}$, and because we do not optimize over
the elements of the diagonal of ${\mathbf{S}}$.

At first sight, computing each $\check{s}_{\ell}$ requires roughly $N^{2}$
operations, so estimating the whole vector ${\mathbf{s}}$ would entail a
computational complexity of ${\mathcal{O}}(N^{4})$. However, a closer
inspection of the vectors $\bm{\sigma}_{\ell}$ reveals that no more than $2N$
of their entries are non-zero because $\bm{\sigma}_{\ell}$ are the columns of
the Kronecker sum of two $N\times N$ matrices. We exploit this sparsity and
reduce the number of operations required to compute each $s_{\ell}$ to
approximately $2N$, rendering the final computational complexity of the graph
denoising step to ${\mathcal{O}}(2\tau_{max_{2}}N^{3})$.

Input: ${\mathbf{X}}$, ${\mathbf{Y}}$, ${\bar{\mathbf{S}}}$

Output : ${\hat{\mathbf{H}}}$, ${\hat{\mathbf{S}}}$.

1 Initialize ${\mathbf{H}}^{(0)}$ and ${\mathbf{S}}^{(0)}$

2 ${\bar{\mathbf{s}}}={\mathrm{vec}}({\bar{\mathbf{S}}})$

3 for  _$t=0$ to $t_{max}-1$_ do

// GF-identification step

4 ${\mathbf{\check{H}}}^{(0)}={\mathbf{H}}^{(t)}$

5 for  _$\tau=0$ to $\tau_{max_{1}}-1$_ do

6 ${\mathbf{\check{H}}}^{(\tau+1)}={\mathbf{\check{H}}}^{(\tau)}+\mu\nabla
f_{1}({\mathbf{\check{H}}}^{(\tau)}|{\mathbf{S}}^{(t)})$

7

8 end for

9 ${\mathbf{H}}^{(t+1)}={\mathbf{\check{H}}}^{(\tau_{max_{1}})}$

10

// Graph denoising step

11
$[\bm{\sigma}_{1},...,\bm{\sigma}_{N^{2}}]={\mathbf{H}}^{(t+1)^{\top}}\oplus{\mathbf{H}}^{(t+1)}$

12 ${\mathbf{\check{s}}}^{(0)}={\mathrm{vec}}({\mathbf{S}}^{(t)})$

13 Update $\bar{\bm{\omega}}^{(t)}$, $\bm{\omega}^{(t)}$ via (17) using
${\bar{\mathbf{s}}}$ and ${\mathbf{\check{s}}}^{(0)}$

14

15 for  _$i=0$ to $\tau_{max_{2}}-1$_ do

16 for  _$\ell\not\in{\mathcal{D}}$_ do

17 Obtain ${\mathbf{r}}_{\ell}^{(\tau)}$ via (30)

18 Obtain $\check{s}_{\ell}^{(\tau+1)}$ via (29) using $\bm{\sigma}_{\ell}$,
${\mathbf{r}}_{\ell}^{(\tau)}$, $\omega_{\ell}$, $\bar{\omega}_{\ell}$

19

20 end for

21

22 end for

23
${\mathbf{S}}^{(t+1)}=\mathrm{unvec}({\mathbf{\check{s}}}^{(\tau_{max_{2}})})$

24 end for

${\hat{\mathbf{H}}}={\mathbf{H}}^{(t_{max})},\;{\hat{\mathbf{S}}}={\mathbf{S}}^{(t_{max})}$.

Algorithm 2 Reduced-complexity robust GF identification.

The pseudocode describing the efficient implementation of Steps 1 and 2 is
provided in Alg. 2. The summary is as follows. We postulate a nested algorithm
with two loops. The outer loop runs $t_{\max}$ iterations. The inner loop runs
two steps: Step 1, with $\tau_{max_{1}}$ iterations, and Step 2, with
$\tau_{\max_{2}}$ iterations. While the complexity for Alg. 1 scaled as
${\mathcal{O}}(t_{\max}N^{7})$, with $t_{\max}$ being typically small, the
overall computational complexity of Alg. 2 is roughly
${\mathcal{O}}(t_{\max}(\tau_{max_{1}}+\tau_{max_{2}})N^{3})$, which is
encouraging, since $2N^{2}$ variables are optimized and it scales with $N$
significantly better than Alg. 1. Solving Steps 1 and 2 optimally requires
setting large values for $\tau_{\max_{1}}$ and $\tau_{\max_{2}}$. Nonetheless,
we observe that in most tested setups the approach of setting small values for
$\tau_{\max_{1}}$ and $\tau_{\max_{2}}$ (at the cost of setting a slightly
higher value for $t_{\max}$) typically yields a faster convergence. Finally,
implementations where the number of iterations is not fixed but selected based
on some convergence criterion are also sensible alternatives.

We close the section noting that we developed Alg. 2 for the setting described
in 1 because the notation was simpler and facilitated the discussion.
Nonetheless, an analogous approach may be followed for the joint estimation of
$K$ GFs (cf. Sec. V), resulting in an algorithm with complexity per GF similar
to that for Alg. 2.

Figure 2: Assessing the performance of the robust GF identification algorithm
and the impact of perturbations in the topology. (a) shows the error of
estimating ${\hat{\mathbf{h}}}$ as the order of the GF increases; (b) and (c)
respectively show the error of estimating ${\hat{\mathbf{H}}}$ and
${\hat{\mathbf{S}}}$ using a robust or a non-robust approach for several types
of perturbations.

##  VII Numerical results

This section discusses several numerical experiments to gain insights and
assess the performance of the robust GF identification algorithms. Unless
specified otherwise, for a variable of interest $\bm{\Theta}$, we report its
normalized estimation error defined as

$nerr(\hat{\bm{\Theta}},\bm{\Theta}):={\|\hat{\bm{\Theta}}-\bm{\Theta}\|_{F}^{2}}/{\|\bm{\Theta}\|_{F}^{2}},$
(31)

where $\hat{\bm{\Theta}}$ and $\bm{\Theta}$ denote the estimated and the true
value, respectively. The code implementing our algorithms and the experiments
presented next is available on
GitHub111https://github.com/reysam93/graph_denoising. The interested reader is
referred there for additional details and tests.

###  VII-A Synthetic experiments

We start by evaluating our algorithms with synthetic data, which is key to
gain intuition. Unless otherwise stated, graphs are sampled from an Erdős
Rényi (ER) random graph model with a link probability of $p=0.2$ and $N=20$
nodes; ${\bar{\mathbf{S}}}$ is obtained by randomly creating and destroying
10% of the links in ${\mathbf{S}}$; $M=50$ signals ${\mathbf{X}}$ and
${\mathbf{Y}}$ are generated according to (2), with the columns of
${\mathbf{X}}$ being drawn from a multivariate Gaussian distribution
${\mathcal{N}}(\mathbf{0},{\mathbf{I}})$, so the signals ${\mathbf{Y}}$ are
stationary on ${\mathbf{S}}$; signals in ${\mathbf{Y}}$ are corrupted with
white Gaussian noise with a normalized power of $\eta_{\mathbf{W}}=0.05$; and
the reported error corresponds to the median of $nerr$ across 64 realizations
of graphs and graph signals.

Test case 1. The first experiment evaluates the influence of perturbations as
the order of the GF $R$ increases. The number of observed pairs of signals
considered is $M=100$ and 10% of the edges in ${\mathbf{S}}$ are perturbed.
Results are reported in Fig. 2(a), where the x-axis shows $R$ and the y-axis
$nerr({\hat{\mathbf{h}}},{\mathbf{h}})$. The algorithms considered are: (i)
the GF identification algorithm that ignores perturbations [see (4)], denoted
as “FI”; (ii) the robust GF identification algorithm from Alg. 1 (“RFI”);
(iii) a variation of “RFI” where the reweighted $\ell_{1}$ norm is replaced by
the standard $\ell_{1}$ norm (“RFI-$\ell_{1}$”); and (iv) the robust GF
identification algorithm accounting for the stationarity of ${\mathbf{Y}}$
(“RFI-ST”). First, we observe that the error of the “FI” algorithm, while
small for low values of $R$, increases rapidly as $R$ grows. This is aligned
with the discussion of high-order polynomials in Sec. III and illustrates the
merits of the robust algorithms. Moreover, “RFI-ST” presents the best
performance illustrating the importance of exploiting additional structure
when it is available. Finally, comparing the error of “RFI” and
“RFI-$\ell_{1}$” showcases the benefits of replacing the $\ell_{1}$ norm with
its reweighted version.

Test case 2. The next experiment tests the influence of different types of
perturbations in the robust and non-robust GF identification algorithms. Figs.
2(b) and 2(c) illustrate the error of the estimated GF ${\hat{\mathbf{H}}}$
and the denoised GSO ${\hat{\mathbf{S}}}$ as the ratio of perturbed links in
${\bar{\mathbf{S}}}$ increases. Graphs are sampled from the small world [47]
random graph model and ${\bar{\mathbf{S}}}$ is obtained by creating new links,
destroying existing links, or simultaneously creating and destroying links,
which are respectively denoted as “C”, “D”, and “C/D” in the legend. Since the
non-robust “FI” algorithm does not perform graph denoising we show the error
$nerr({\bar{\mathbf{S}}},{\mathbf{S}})$, denoted as “${\bar{\mathbf{S}}}$” in
Fig. 2(c). Furthermore, because the number of perturbed links is fixed, the
error of ${\bar{\mathbf{S}}}$ is the same for the considered perturbations and
it is only plotted once. From the figures, we observe that destroying links is
the most harmful perturbation, especially when the focus is on
${\hat{\mathbf{S}}}$. This may be explained because destroying links is prone
to produce non-connected graphs. Nevertheless, the results show the resilience
of the “RFI” algorithm, which provides low-error estimates
${\hat{\mathbf{H}}}$ and ${\hat{\mathbf{S}}}$ even when more than 20% of the
links are perturbed.

Figure 3: Comparing the performance of several robust GF identification
algorithms. (a) shows the error of ${\hat{\mathbf{H}}}$ when estimated with
the proposed algorithm and with other baselines as the ratio of perturbed
links increases. Different graph-signal models are considered. (b) and (c)
respectively show the running time and error of ${\hat{\mathbf{H}}}$ using
Alg. 1 and Alg. 2 as the number of nodes increases. Different values for the
maximum number of iterations of the inner loops are considered.

Test case 3. Next, we compare the performance of our algorithms with other
robust alternatives. Fig. 3(a) reports, for each algorithm,
$nerr({\hat{\mathbf{H}}},{\mathbf{H}})$ as the ratio of perturbed links
increases. The baselines considered are the TLS-SEM algorithm from [26], and
LLS-SCP from [28]. We note that the TLS-SEM algorithm is tailored to graph
signals following a SEM of the form

${\mathbf{Y}}={\mathbf{A}}{\mathbf{Y}}+{\mathbf{X}}=({\mathbf{I}}-{\mathbf{A}})^{-1}{\mathbf{X}},$
(32)

where the observations at the $i$-th node are represented by the values of the
neighbors of $i$ and an exogenous input. As a result, the TLS-SEM algorithm
may not be well suited to deal with signals generated according to the more
general model in (2). Taking this into account, to offer a more favorable
comparison we consider two types of graph signals: (i) signals generated
according to (32), denoted as “SEM”; and (ii) signals generated according to
(2), denoted as “H”. It is worth noting that the “SEM” can be considered as a
particular case of the model “H” when the GF
${\mathbf{H}}_{SEM}=({\mathbf{I}}-{\mathbf{A}})^{-1}$ is employed.

Looking at the results in Fig. 3(a) we observe the following. When the “SEM”
model is considered, TLS-SEM (denoted as “TLS-SEM”) obtains the best
performance when the perturbation probability is small, and then, the
performance of “TLS-SEM” and that of the “RFI” algorithm become comparable.
This illustrates that our algorithm is especially suitable to deal with a
large number of perturbed links. On the other hand, when the “H” model is
considered, we observe that the “RFI” algorithm consistently outperforms the
baselines in the presence of perturbations. The good performance of the “RFI”
algorithm on both signal models highlights the flexibility of the proposed
formulation since it considers more lenient assumptions than the other
alternatives.

Test case 4. Now, we compare the performance of the standard and the efficient
implementation of the robust identification algorithm, as described in Algs. 1
and 2. The results are shown in Figs. 3(b) and 3(c), where the figures depict
the running time measured in seconds and
$nerr({\hat{\mathbf{H}}},{\mathbf{H}})$ as $N$ increases. The legend
identifies first the algorithm employed, then the number of iterations of the
outer loop ($t_{max}$), and finally the iterations of the inner loops (with
$\tau_{{max_{1}}}=\tau_{{max_{2}}}$). As expected, Fig. 3(b) shows that Alg. 2
is remarkably faster than Alg. 1 even with medium-sized graphs, achieving a
running time $10^{3}$ times smaller when $N=100$. On the other hand, in Fig.
3(c) we observe that “Eff-5-50” has an error that is close to the standard
implementation (“Stand-5”) even though it is considerably faster. Furthermore,
the trade-off between speed and estimation accuracy is also evident.
“Eff-5-10” is the fastest implementation but the quality of its estimated GF
may not be enough for graphs with more than 40 nodes.

Test case 5. The last experiment with synthetic data studies the benefits of
the joint GF estimation. All the GFs are polynomials of the same
${\mathbf{S}}$, and for each ${\mathbf{H}}_{k}$ we consider $M_{k}=15$ noisy
observations with $\eta_{\mathbf{w}}=0.01$. Fig. 4 shows the results, with the
y-axis being the normalized error averaged across the $K$ graphs, i.e.,
$\frac{1}{K}\sum_{k=1}^{K}nerr({\hat{\mathbf{H}}}_{k},{\mathbf{H}}_{k})$, and
the x-axis representing $K$. We compare the performance of estimating the GFs
jointly (marked as “J” in the legend) or separately for the three algorithms
(“RFI-$\ell_{1}$”, “RFI”, and “RFI-st”) described in Test case 1. Note that
“RFI-J” corresponds to the formulation in (V). The first thing we observe from
the results in Fig. 4 is that the error decreases as $K$ increases when a
joint algorithm is employed. This is aligned with the discussion in Sec. V and
illustrates the benefit of exploiting the common structure. In addition,
algorithms accounting for the stationary of ${\mathbf{Y}}$ outperform the non-
stationary alternatives even though we only have $M=15$ signals to estimate
the covariance ${\hat{\mathbf{C}}}_{\mathbf{y}}$.

Figure 4: Error performance when estimating $K$ GFs using the separate and
joint approach for different values of $K$.

###  VII-B Real-world datasets

To close the numerical evaluation, we test our robust GF identification
algorithms over two real-world datasets.

Figure 5: Performance of the algorithms predicting ozone levels in the AirData station network, as the time horizon of the prediction increases. Models | 1-Step | 3-Step
---|---|---
TTS=0.25 | TTS = 0.5 | TTS=0.25 | TTS = 0.5
LS | $6.9\cdot 10^{-3}$ | $3.1\cdot 10^{-3}$ | $2.1\cdot 10^{-2}$ | $9.1\cdot 10^{-3}$
LS-GF | $3.3\cdot 10^{-3}$ | $3.3\cdot 10^{-3}$ | $8.4\cdot 10^{-3}$ | $8.5\cdot 10^{-3}$
TLS-SEM | $4.0\cdot 10^{1}$ | $3.7\cdot 10^{-2}$ | $6.8\cdot 10^{-1}$ | $5.5\cdot 10^{-2}$
RFI | $3.4\cdot 10^{-3}$ | $3.1\cdot 10^{-3}$ | $8.5\cdot 10^{-3}$ | $7.5\cdot 10^{-3}$
AR(3)-RFI | $3.2\cdot 10^{-3}$ | $2.8\cdot 10^{-3}$ | $7.8\cdot 10^{-3}$ | $6.9\cdot 10^{-3}$
TABLE I: Performance of the algorithms in predicting the temperature for 2
prediction horizons (1 and 3) and 2 values (25% and 50%) of train-test split
(TTS). The metrics shown are the average of the normalized error at each
timestep
$\frac{1}{M}\sum_{\kappa=1}^{M}nerr({\hat{\mathbf{y}}}_{\kappa},{\mathbf{y}}_{\kappa})$
for all samples.

Weather station network. This test case evaluates the ability of our
algorithms to predict the temperature measured by a network of stations using
the data from previous days. The data comes from the “Global Summary of the
Day” dataset of the National Centers for Environmental
Information222https://www.ncei.noaa.gov/data/global-summary-of-the-
day/archive/ and we used daily temperature measurements from $N=17$ stations
in California during 2017 & 2018\. Specifically, with
${\mathbf{y}}_{\kappa}\in{\mathbb{R}}^{N}$ collecting the measurements of the
17 stations at day $\kappa$, we consider an AR model without exogenous inputs,
so that
${\mathbf{y}}_{\kappa}\approx\sum_{k=1}^{K}{\mathbf{H}}_{k}{\mathbf{y}}_{\kappa-k}$.
The data samples were divided into two subsets, the first one (training) was
used to obtain the GFs ${\mathbf{H}}_{k}$ and the second one (evaluation) was
used to assess the performance and the generalization power of the GFs
obtained. Also, the data is normalized so that the signal at each station for
all time samples has unitary norm.

The underlying ${\mathcal{G}}$ was constructed as the unweighted 5-nearest
neighbors graph, using the geographical distance between stations. Since
temperature relations across stations are likely to be due to a range of
factors (including, e.g., altitude), the considered adjacency (based only on
geographical positions) may be imperfect, rendering our robust algorithms
better suited for this task.

The estimation performance of the different algorithms is shown in Table I.
Since in this case the ground-truth GF is not known, we use the signal
denoising error $nerr({\mathbf{y}}_{\kappa},{\hat{\mathbf{y}}}_{\kappa})$ to
assess the quality of the schemes. In this specific experiment, the error is
measured over all samples (both training and test subsets), to see a clear
downward trend when increasing the number of training samples, or
equivalently, the train-test split (TTS) value. The algorithms evaluated are
“LS”, “LS-GF” (which postulates a GF with coefficients
${\hat{\mathbf{h}}}=\operatornamewithlimits{argmin}_{\mathbf{h}}\|{\mathbf{Y}}-\sum_{\ell}h_{\ell}{\mathbf{S}}^{\ell}{\mathbf{X}}\|_{F}^{2}$),
“TLS-SEM”, “RFI” (which assumes an AR(1) process) and “AR(3)-RFI”. Two values
of TTS (0.25 and 0.50) and two prediction horizons (1 and 3) are considered.
The main observation is that “AR(3)-RFI” yields the best performance in all
settings. Additionally, the results for TTS=0.25 demonstrate the benefits of
considering the underlying graph in the low-sample regime, since even “LS-GF”,
which relies on the imperfect $\bar{{\mathbf{S}}}$, outperforms “LS”. On the
other hand, “LS-GF” does not seem to improve its prediction as TTS increases,
while our two algorithms yield a lower prediction error.

Air quality station network. We consider an experimental setup (AR model,
graph creation method…) similar to that for the weather station data but, in
this case, we use 2018 & 2019 data from the United States Environmental
Protection Agency333https://www.epa.gov/outdoor-air-quality-data to predict
the ozone levels in a network of 17 outdoor stations in California. The
stations chosen were those with at least 330 measurements each year for a
selection of pollutants, and missing data was filled via first-order
interpolation.

The goal here is to analyze how the prediction horizon affects the prediction
error. The value of TTS chosen was 0.5, i.e. evaluation data represented 50%
of the samples. Fig. 5 shows the performance of the algorithms when predicting
ozone levels. As a baseline, “LS-Eval-(LB)” shows the error measured on the
evaluation data when obtaining the GF also using evaluation data, therefore
representing a lower bound for the LS error using AR models of order 1. Also,
“Copy-Prev-Day” represents the error obtained by the “identity GF”, which
copies the previous day’s measurement. As in the previous example, the best
performing algorithm is “AR(3)-RFI”, whose performance is close to the
baseline, followed by “RFI”.

##  VIII Concluding remarks

This paper put forth a framework dealing with estimation problems in GSP where
the information about (the links of) the supporting graph is uncertain.
Specifically, we addressed the problem of estimating a GF (i.e., a polynomial
of the GSO) from input and output graph signals under the key assumption that
only a perturbed version of the true GSO was available. In contrast to the
majority of existing approaches that operate on the spectral domain, we recast
the true graph as an additional estimation variable and formulated an
optimization problem that _jointly_ estimated the GF and the true (unknown)
GSO. We focused first on the case where only one GF needed to be estimated
and, then, shifted to (multi-feature and AR graph signal) setups where
multiple GFs have to be jointly identified. The formulated optimizations
operated completely in the vertex domain and bypassed the problem of computing
high-order polynomials, avoiding the challenges of dealing with the influence
of perturbations in the graph spectrum as well as the numerical instability
and error propagation associated with high-order matrix polynomials. While
non-convex, upon blending techniques from alternating optimization and MM, the
proposed algorithm was shown to be capable to find a stationary point in
polynomial time. This algorithm was later modified so that the scaling of the
computational complexity with respect to the number of nodes in the graph is
reduced. Future work includes delving into the robust estimation of ARMA time-
varying graph signals, consideration of additional graph perturbation models,
and application of our robust estimation framework to other GSP problems, to
name a few.

## Appendix A: Proof of Th. 1

The proof relies on the results presented in [39, Th. 1b], so it suffices to
show that our formulation and algorithm fulfill the required conditions in
[39]. To that end, recall that $f({\mathbf{z}})$ is the objective function in
(13), and let ${\mathbf{z}}_{1}:={\mathrm{vec}}({\mathbf{H}})$ and
${\mathbf{z}}_{2}:={\mathrm{vec}}({\mathbf{S}})$ denote the $B=2$ blocks of
variables considered in our algorithm. Moreover, at each step, the function
$f({\mathbf{z}})$ is approximated by $u_{1}({\mathbf{z}}_{1})$ and
$u_{2}({\mathbf{z}}_{2})$, corresponding to the objective functions in (14)
and (IV-A). Then, to ensure the convergence of our iterative algorithm the
following conditions are required.

(C1) Each function $u_{b}({\mathbf{z}}_{b})$ must be a global upper bound of
$f({\mathbf{z}})$ and the first-order behavior of $u_{b}({\mathbf{z}}_{b})$
and $f({\mathbf{z}})$ must be the same.

(C2) $f({\mathbf{z}})$ must be regular (cf. [39]) at every point in
${\mathcal{Z}}^{*}$.

(C3) The level set
${\mathcal{Z}}^{(0)}=\\{{\mathbf{z}}\;|\;f({\mathbf{z}})\leq
f({\mathbf{z}}^{(0)})\\}$ is compact.

(C4) At least one of the problems in (14) and (IV-A) must have a unique
solution.

Next, we address each of the four conditions separately, proving that our
approach satisfies all of them.

Condition (C1) requires the surrogate functions $u_{b}({\mathbf{z}}_{b})$ to
be global upper bounds of $f({\mathbf{z}})$. For the first block ($b=1$), it
is easy to see that $u_{1}({\mathbf{z}}_{1})=f({\mathbf{z}})$ when the block
${\mathbf{z}}_{2}$ remains constant, so it satisfies the requirements.
Regarding $u_{2}({\mathbf{z}}_{2})$, we approximate $f({\mathbf{z}})$ with the
first-order Taylor series of the logarithmic penalty. Because the $\log$ is a
concave differentiable function, it follows that its Taylor series of order
one constitutes a global upper bound. Moreover, because
$u_{2}({\mathbf{z}}_{2})$ is a first-order Taylor series approximation of
$f({\mathbf{z}})$, it also follows that the first-order behavior of
$f({\mathbf{z}})$ and $u_{2}({\mathbf{z}}_{2})$ is the same. Therefore,
$u_{2}$ also satisfies the requirement, and hence, (C1) is fulfilled.

To prove (C2), according to [39], a function $f({\mathbf{z}})$ is regular if
its non-smooth components are separable across the different blocks of
variables. To show this, we decompose $f$ as $f=g_{A}+g_{B}$, with functions
$g_{A}$ and $g_{B}$ being defined as

$\displaystyle
g_{A}({\mathbf{H}},{\mathbf{S}})=\|{\mathbf{Y}}-{\mathbf{H}}{\mathbf{X}}\|_{F}^{2}+\gamma\|{\mathbf{H}}{\mathbf{S}}-{\mathbf{S}}{\mathbf{H}}\|_{F}^{2},$
$\displaystyle
g_{B}({\mathbf{S}})=\lambda\sum_{i,j=1}^{N}\log(|S_{ij}\\!-\\!\bar{S}_{ij}|+\delta_{2})+\beta\sum_{i,j=1}^{N}\log(|S_{ij}|+\delta_{1}).$

Since $g_{A}$ is a smooth function and the non-smooth function $g_{B}$ only
depends on the variables on the second block,
${\mathbf{z}}_{2}={\mathrm{vec}}({\mathbf{S}})$, it follows that
$f({\mathbf{z}})$ is a regular function for all feasible points.

Next, we show that the level set
${\mathcal{Z}}^{(0)}=\\{{\mathbf{z}}\;|\;f({\mathbf{z}})\leq
f({\mathbf{z}}^{(0)})\\}$ is compact as required by (C3). We start by noting
that the entries of ${\mathbf{S}}$ are continuous subsets of ${\mathbb{R}}$,
(e.g., $S_{ij}\in{\mathbb{R}}_{+}$ when ${\mathbf{S}}={\mathbf{A}}$), and that
${\mathbf{H}}\in{\mathbb{R}}^{N\times N}$, so $f({\mathbf{z}})$ is continuous.
Moreover, $f({\mathbf{z}})\leq f({\mathbf{z}}^{(0)})$ implies that the
functions $\|{\mathbf{Y}}-{\mathbf{H}}{\mathbf{X}}\|_{F}^{2}$ and
$\log(|S_{ij}|+\delta_{1})$ are all bounded, rendering the domain of
$f({\mathbf{z}})$ bounded. It follows then that the level set
${\mathcal{Z}}^{(0)}$ is compact.

Finally, we need to prove that either (14) or (IV-A) has a unique solution, so
that (C4) is fulfilled. Prop. 1 (see below) states that, under the two
conditions required by Th. 1 (i.e., ${\mathbf{S}}$ does not have repeated
eigenvalues, and the graph signals ${\mathbf{X}}$ excite every graph
frequency), the solution to (14) is unique. This confirms that (C4) is
satisfied, concluding the proof.

######  Proposition 1.

Let ${\mathbf{H}}\in{\mathbb{R}}^{N\times N}$,
${\mathbf{S}}={\mathbf{V}}\text{diag}(\bm{\lambda}){\mathbf{V}}^{-1}\in{\mathbb{R}}^{N\times
N}$, and ${\mathbf{X}}\in{\mathbb{R}}^{N\times M}$ be the GF, the GSO, and the
input signals in (14). Then, (14) has a unique solution w.r.t. ${\mathbf{H}}$
if the following conditions are satisfied:

  1. 1.

$\lambda_{i}\neq\lambda_{i^{\prime}}$, for all $i\neq i^{\prime}$ and
$(i,i^{\prime})\in\\{1,...,N\\}^{2}$.

  2. 2.

Every row of ${\tilde{\mathbf{X}}}={\mathbf{V}}^{-1}{\mathbf{X}}$ has at least
one non-zero entry.

###### Proof.

To simplify exposition, we focus first on the (most restrictive) setup of
having only $M=1$ input-output pairs. Defining
${\hat{\mathbf{h}}}:={\mathrm{vec}}({\mathbf{H}})$, we can reformulate (14) as

$\text{min}_{{\hat{\mathbf{h}}}\in{\mathbb{R}}^{N^{2}}}\gamma\|({\mathbf{I}}\otimes{\mathbf{S}}-{\mathbf{S}}^{\top}\otimes{\mathbf{I}}){\hat{\mathbf{h}}}\|_{2}^{2}+\|{\mathbf{y}}-({\mathbf{x}}^{\top}\otimes{\mathbf{I}}){\hat{\mathbf{h}}}\|_{2}^{2},$
(33)

where lowercase symbols ${\mathbf{y}}$ and ${\mathbf{x}}$ are used to
emphasize that the output and input signals are a single $N$-dimensional
vector. Upon defining
${\mathbf{D}}:={\mathbf{I}}\otimes{\mathbf{S}}-{\mathbf{S}}^{\top}\otimes{\mathbf{I}}$,
and ${\mathbf{E}}:={\mathbf{x}}^{\top}\otimes{\mathbf{I}}$, solving (33) is
equivalent to solving

$\text{min}_{{\hat{\mathbf{h}}}\in{\mathbb{R}}^{N^{2}}}\Big{\|}\begin{bmatrix}{\mathbf{0}}_{N^{2}}\\\
{\mathbf{y}}~{}\end{bmatrix}-{\mathbf{F}}{\hat{\mathbf{h}}}\Big{\|}_{2}^{2}~{}\text{with}~{}{\mathbf{F}}:=\begin{bmatrix}\gamma{\mathbf{D}}\\\
~{}{\mathbf{E}}\end{bmatrix}$ (34)

To prove that (34) has a unique solution, it suffices to show that
${\mathbf{F}}$ is full column rank, i.e.
$\nexists\;{\mathbf{n}}\in{\mathbb{R}}^{N^{2}}$ such that
${\mathbf{F}}{\mathbf{n}}={\mathbf{0}}_{N+N^{2}}$. To show this, we first
identify ${\mathcal{N}}({\mathbf{D}})$, the null space of ${\mathbf{D}}$, and
then show that ${\mathbf{E}}{\mathbf{n}}\neq{\mathbf{0}}_{N}\;\text{for
all}\;{\mathbf{n}}\in{\mathcal{N}}({\mathbf{D}})\setminus\\{{\mathbf{0}}_{N^{2}}\\}$.

We start with the characterization of ${\mathcal{N}}({\mathbf{D}})$. Given the
Kronecker structure of ${\mathbf{D}}$, each of its $N^{2}$ eigenvalues has the
form $\lambda_{k}-\lambda_{k^{\prime}}$, with
$({\mathbf{V}}^{-1})^{\top}\otimes{\mathbf{V}}$ being the associated
eigenvectors. Leveraging that $\lambda_{i}\neq\lambda_{i^{\prime}}$ for $i\neq
i^{\prime}$, it follows that only when $i=i^{\prime}$ the eigenvalue of
${\mathbf{D}}$ is zero. As a result, $\text{rank}({\mathbf{D}})=N^{2}-N$ and
$\text{dim}({\mathcal{N}}({\mathbf{D}}))=N$. Equally important, the $N$
eigenvectors associated with the $N$ zero eigenvalues are given by
$({\mathbf{V}}^{-1})^{\top}\odot{\mathbf{V}}$, which, as a result, constitutes
a basis spanning ${\mathcal{N}}({\mathbf{D}})$. More formally, we concluded
that
${\mathcal{N}}({\mathbf{D}})=\\{(({\mathbf{V}}^{-1})^{\top}\odot{\mathbf{V}})\bm{\theta}\;|\forall\;\bm{\theta}\in{\mathbb{R}}^{N}\\}$.

Thus, to show that ${\mathbf{F}}$ in (34) is full column rank we just need to
prove that the only element ${\mathbf{n}}\in{\mathcal{N}}({\mathbf{D}})$ that
renders ${\mathbf{E}}{\mathbf{n}}={\mathbf{0}}_{N}$ is the all-zero vector
${\mathbf{0}}_{N^{2}}$. To do so, we leverage the characterization of
${\mathcal{N}}({\mathbf{D}})$ and write ${\mathbf{E}}{\mathbf{n}}$ as

$\displaystyle{\mathbf{E}}{\mathbf{n}}$
$\displaystyle=({\mathbf{x}}^{\top}\otimes{\mathbf{I}})(({\mathbf{V}}^{-1})^{\top}\odot{\mathbf{V}})\bm{\theta}=({\mathbf{x}}^{\top}({\mathbf{V}}^{-1})^{\top}\odot{\mathbf{V}})\bm{\theta}$
$\displaystyle={\mathbf{V}}\text{diag}(\bm{\theta})({\mathbf{x}}^{\top}({\mathbf{V}}^{-1})^{\top})^{\top}={\mathbf{V}}\text{diag}(\bm{\theta}){\mathbf{V}}^{-1}{\mathbf{x}}$
$\displaystyle={\mathbf{V}}\text{diag}(\bm{\theta}){\tilde{\mathbf{x}}}={\mathbf{V}}(\bm{\theta}\circ{\tilde{\mathbf{x}}}),$
(35)

where we used the property $(a\otimes b)(c\odot d)=ac\odot bd$. Since
${\mathbf{V}}$ is invertible, the first and last terms in (35) demonstrate
that ${\mathbf{E}}{\mathbf{n}}={\mathbf{0}}_{N}$ requires
$\bm{\theta}\circ{\tilde{\mathbf{x}}}={\mathbf{0}}_{N}$. However, condition 2)
in Prop. 1 states that $\tilde{x}_{i}\neq 0~{}\text{for all}~{}i$; hence,
$\bm{\theta}\circ{\tilde{\mathbf{x}}}={\mathbf{0}}_{N}$ requires
$\bm{\theta}={\mathbf{0}}_{N}$. This implies that the only element in
${\mathcal{N}}({\mathbf{D}})$ that renders
${\mathbf{E}}{\mathbf{n}}={\mathbf{0}}_{N}$ is
${\mathbf{n}}=({\mathbf{V}}^{-1})^{\top}\odot{\mathbf{V}}){\mathbf{0}}_{N}={\mathbf{0}}_{N^{2}}$,
concluding the proof.

The proof can be generalized for $M>1$. In that case, the matrix
${\mathbf{E}}$ has size $MN\times N^{2}$ and the counterpart to (35)
establishes that having ${\mathbf{E}}{\mathbf{n}}={\mathbf{0}}$ requires
${\mathrm{vec}}(\text{diag}(\bm{\theta}){\tilde{\mathbf{X}}})={\mathbf{0}}_{MN}$.
Since Prop. 1 assumes that each row of ${\tilde{\mathbf{X}}}$ has at least one
nonzero entry, it follows that $\bm{\theta}={\mathbf{0}}_{MN}$, concluding the
proof. ∎

## References

  * [1] S. Rey and A. G. Marques, “Robust graph-filter identification with graph denoising regularization,” in _IEEE Intl. Conf. Acoustics, Speech Signal Process. (ICASSP)_. IEEE, 2021, pp. 5300–5304.
  * [2] K. Nodop, R. Connolly, and F. Girardi, “The field campaigns of the European Tracer Experiment (etex): Overview and results,” _Atmospheric Environment_ , vol. 32, no. 24, pp. 4095–4108, 1998.
  * [3] E. D. Kolaczyk, _Statistical Analysis of Network Data: Methods and Models_. New York, NY: Springer, 2009.
  * [4] O. Sporns, _Discovering the Human Connectome_. MIT Press, 2012.
  * [5] D. Shuman et al., “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, 2013\.
  * [6] A. Sandryhaila and J. M. F. Moura, “Discrete signal processing on graphs: Frequency analysis,” _IEEE Trans. Signal Process._ , vol. 62, no. 12, pp. 3042–3054, 2014.
  * [7] P. Djuric and C. Richard, _Cooperative and Graph Signal Processing: Principles and Applications_. Academic Press, 2018.
  * [8] F. J. Iglesias et al., “Demixing and blind deconvolution of graph-diffused sparse signals,” in _IEEE Intl. Conf. Acoustics, Speech Signal Process. (ICASSP)_. IEEE, 2018, pp. 4189–4193.
  * [9] A. Ortega et al., “Graph signal processing: Overview, challenges, and applications,” _Proc. IEEE_ , vol. 106, no. 5, pp. 808–828, 2018.
  * [10] A. G. Marques et al., “Graph signal processing: Foundations and emerging directions (editorial),” _IEEE Signal Process. Mag._ , vol. 37, Nov. 2020\.
  * [11] A. Sandryhaila and J. M. F. Moura, “Discrete signal processing on graphs,” _IEEE Trans. Signal Process._ , vol. 61, no. 7, pp. 1644–1656, 2013.
  * [12] 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.
  * [13] J. Liu, E. Isufi, and G. Leus, “Filter design for autoregressive moving average graph filters,” _IEEE Trans. Signal Process. Inf. Netw._ , vol. 5, no. 1, pp. 47–60, 2018.
  * [14] 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, 2017.
  * [15] S. Rey, F. J. 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.
  * [16] Y. Zhu, F. J. Iglesias-Garcia, 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.
  * [17] Y. He and H. Wai, “Detecting central nodes from low-rank excited graph signals via structured factor analysis,” _IEEE Trans. Signal Process._ , 2022.
  * [18] E. Isufi, A. Loukas, A. Simonetto, and G. Leus, “Filtering random graph processes over random time-varying graphs,” _IEEE Trans. Signal Process._ , vol. 65, no. 16, pp. 4406–4421, 2017.
  * [19] J. Friedman, T. Hastie, and R. Tibshirani, “Sparse inverse covariance estimation with the graphical lasso,” _Biostatistics_ , vol. 9, no. 3, pp. 432–441, 2008.
  * [20] X. Dong, D. Thanou, P. Frossard, and P. Vandergheynst, “Learning Laplacian matrix in smooth graph signal representations,” _IEEE Trans. Signal Process._ , vol. 64, no. 23, pp. 6160–6173, 2016.
  * [21] 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.
  * [22] A. Buciulea, S. Rey, and A. G. Marques, “Learning graphs from smooth and graph-stationary signals with hidden variables,” _IEEE Trans. Signal Inf. Process. Netw._ , vol. 8, pp. 273–287, 2022.
  * [23] E. Ceci and S. Barbarossa, “Graph signal processing in the presence of topology uncertainties,” _IEEE Trans. Signal Process._ , vol. 68, pp. 1558–1573, 2020.
  * [24] S. Segarra and A. Ribeiro, “Stability and continuity of centrality measures in weighted graphs,” _IEEE Trans. Signal Process._ , vol. 64, no. 3, pp. 543–555, 2015.
  * [25] J. Miettinen, S. A. Vorobyov, and E. Ollila, “Modelling graph errors: Towards robust graph signal processing,” _arXiv preprint arXiv:1903.08398_ , 2019.
  * [26] E. Ceci, Y. Shen, G. B. Giannakis, and S. Barbarossa, “Graph-based learning under perturbations via total least-squares,” _IEEE Trans. Signal Process._ , vol. 68, pp. 2870–2882, 2020.
  * [27] J. Miettinen, S. A. Vorobyov, and E. Ollila, “Graph error effect in graph signal processing,” in _IEEE Intl. Conf. Acoustics, Speech Signal Process. (ICASSP)_. IEEE, 2018, pp. 4164–4168.
  * [28] A. Natali, M. Coutino, and G. Leus, “Topology-aware joint graph filter and edge weight identification for network processes,” in _IEEE Intl. Wrkshp. Mach. Learn. Signal Process. (MLSP)_. IEEE, 2020, pp. 1–6.
  * [29] V. M. Tenorio et al., “A robust alternative for graph convolutional neural networks via graph neighborhood filters,” in _Conf. Signals, Syst., Computers (Asilomar)_. IEEE, 2021, pp. 1573–1578.
  * [30] R. Levie, E. Isufi, and G. Kutyniok, “On the transferability of spectral graph filters,” in _Int. Conf. Sampling Theory Appl. (SampTA)_. IEEE, 2019, pp. 1–5.
  * [31] R. Levie et al., “Transferability of spectral graph convolutional neural networks.” _J. Mach. Learn. Res._ , vol. 22, pp. 272–1, 2021.
  * [32] L. Ruiz, F. Gama, and A. Ribeiro, “Graph neural networks: Architectures, stability, and transferability,” _Proc. IEEE_ , vol. 109, no. 5, pp. 660–682, 2021.
  * [33] 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.
  * [34] H. Zou and T. Hastie, “Regularization and variable selection via the elastic net,” _J. Roy. Statistical Soc.: Ser. B (Statistical Methodology)_ , vol. 67, no. 2, pp. 301–320, 2005.
  * [35] B. Dai, S. Ding, and G. Wahba, “Multivariate Bernoulli distribution,” _Bernoulli_ , vol. 19, no. 4, pp. 1465–1483, 2013.
  * [36] E. J. Candes, M. B. Wakin, and S. P. Boyd, “Enhancing sparsity by reweighted $\ell_{1}$ minimization,” _J. Fourier Anal. Appl._ , vol. 14, no. 5-6, pp. 877–905, 2008.
  * [37] Y. Sun, P. Babu, and D. P. Palomar, “Majorization-minimization algorithms in signal processing, communications, and machine learning,” _IEEE Trans. Signal Process._ , vol. 65, no. 3, pp. 794–816, 2016.
  * [38] J. Gorski, F. Pfeuffer, and K. Klamroth, “Biconvex sets and optimization with biconvex functions: A survey and extensions,” _Math. Methods Oper. Res._ , vol. 66, no. 3, pp. 373–407, 2007.
  * [39] M. Hong, M. Razaviyayn, Z. Luo, and J. Pang, “A unified algorithmic framework for block-structured optimization involving big data: With applications in machine learning and signal processing,” _IEEE Signal Process. Mag._ , vol. 33, no. 1, pp. 57–77, 2015.
  * [40] R. Shafipour and G. Mateos, “Online topology inference from streaming stationary graph signals with partial connectivity information,” _Algorithms_ , vol. 13, no. 9, p. 228, 2020.
  * [41] G. C. Reinsel, _Elements of multivariate time series analysis_. Springer Science & Business Media, 2003.
  * [42] J. Mei and J. M. F. Moura, “Signal processing on graphs: Causal modeling of unstructured data,” _IEEE Trans. Signal Process._ , vol. 65, no. 8, pp. 2077–2092, 2017.
  * [43] E. Isufi, A. Loukas, N. Perraudin, and G. Leus, “Forecasting time series with varma recursions on graphs,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 18, pp. 4870–4885, 2019.
  * [44] S. Boyd and L. Vandenberghe, _Convex Optimization_. Cambridge University Press, 2004.
  * [45] D. Coppersmith and S. Winograd, “Matrix multiplication via arithmetic progressions,” in _Proc. ACM Symp. Theory Comput._ , 1987, pp. 1–6.
  * [46] T. Hastie, R. Tibshirani, and M. Wainwright, “Statistical learning with sparsity,” _Monographs Statistics Appl. Probability_ , vol. 143, p. 143, 2015\.
  * [47] M. E. J. Newman and D. J. Watts, “Renormalization group analysis of the small-world network model,” _Phys. Lett. A_ , vol. 263, no. 4-6, pp. 341–346, 1999.