# Polynomial Graphical Lasso: Learning Edges from Gaussian Graph-Stationary
Signals

Andrei Buciulea, , Jiaxi Ying, , Antonio G. Marques, , and Daniel P. Palomar,
Work partially supported by the Spanish Spanish AEI Grants
PID2019-105032GB-I00 and PID2022-136887NB-I00, and by the Autonomous Community
of Madrid within the ELLIS Unit Madrid framework, URJC/CAM grant F861 and
URJC/CAM grant PREDOC20-003. Work partially supported by the Hong Kong GRF
16206123 research grant. An early preliminary version of this work was
presented as a conference paper in [1]. (Corresponding author: Antonio G.
Marques) Andrei Buciulea and Antonio G. Marques are with the Dept. of Signal
Theory and Comms., King Juan Carlos University, 28933 Madrid, Spain
(e-mail:andrei.buciulea, antonio.garcia.marques@urjc.es). Jiaxi Ying is with
the Department of Mathematics, Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong (e-mail:jx.ying@connect.ust.hk). Daniel P.
Palomar is with the Department of Electronic and Computer Engineering and the
Department of Industrial Engineering and Decision Analytics, Hong Kong
University of Science and Technology, Kowloon, Hong Kong (e-mail:
palomar@ust.hk).

###### Abstract

This paper introduces Polynomial Graphical Lasso (PGL), a new approach to
learning graph structures from nodal signals. Our key contribution lies in
modeling the signals as Gaussian and stationary on the graph, enabling the
development of a graph-learning formulation that combines the strengths of
graphical lasso with a more encompassing model. Specifically, we assume that
the precision matrix can take any polynomial form of the sought graph,
allowing for increased flexibility in modeling nodal relationships. Given the
resulting complexity and nonconvexity of the resulting optimization problem,
we (i) propose a low-complexity algorithm that alternates between estimating
the graph and precision matrices, and (ii) characterize its convergence. We
evaluate the performance of PGL through comprehensive numerical simulations
using both synthetic and real data, demonstrating its superiority over several
alternatives. Overall, this approach presents a significant advancement in
graph learning and holds promise for various applications in graph-aware
signal analysis and beyond.

###### Index Terms:

Network-topology inference, graph learning, Gaussian signals, graph stationary
signals, covariance estimation.

##  I Introduction

Modern datasets often exhibit an irregular non-Euclidean support. In such
scenarios, graphs have emerged as a pivotal tool, facilitating the
generalization of classical information processing and structured learning
techniques to irregular domains. Today, there is a wide range of applications
that leverage graphs when processing, learning, and extracting knowledge from
their associated datasets (see, e.g., problems in the context of electrical,
communication, social, geographic, financial, genetic, and brain networks [2,
3, 4, 5, 6], to name a few). When using graphs to process structured non-
Euclidean data, it is usually assumed that the underlying network topology is
known. Unfortunately, this is not always the case. In many cases, the
structure of the graph is not well defined, either because there is no
underlying physical network or because the (best) metric to assess the level
of association between the nodes is not known.

Since in most cases the existing relationships are not known beforehand, the
standard approach is to infer the structure of the network from a set of
available nodal observations/signals/features. To estimate the interactions
between the existing nodes, the first step is to formally define the
relationship between the topology of the graph and the properties of the
signals defined on top of it. Early graph topology inference methods [7, 8]
adopted a statistical approach, such as the correlation network [2, Ch.
7.3.1], partial correlations, or Gaussian Markov random fields (GMRF), with
the latter leading to the celebrated graphical Lasso (GL) scheme [9, 2].
Partial correlation methods have been generalized to nonlinear settings [10].
Also in the nonlinear realm, less rigorous approaches simply postulate a
similarity scores, with links being drawn if their score exceeds a given
threshold. In recent years, graph signal processing (GSP) based models [11,
12, 13] have brought new ideas to the field, considering more complex
relationships between the signals and the sought graph. These approaches have
been generalized to deal with more complex scenarios that often arise in
practice, such as the presence of hidden variables [14, 15, 16] or the
simultaneous inference of multiple networks [17, 18].

The existing graph-learning methods exhibit different pros and cons, with
relevant tradeoffs including computational complexity, expressiveness, model
accuracy, or sample complexity, to name a few. For instance, correlation
networks need very few samples and can be run in parallel for each pair of
nodes, but fail to capture intermediation nodal effects. On the other hand, GL
(a maximum likelihood estimator for GMRF) can handle the intermediation effect
while still requiring a relatively small number of samples compared to the
size of the network. Some disadvantages of GL include assuming a relatively
simple signal model (failing to deal with, e.g., linear autoregressive network
models) and forcing the learned graph to be a positive definite matrix. To
overcome some of these issues, [13] proposed a more general model that assumed
that the signals were stationary in the network (GSR) or, equivalently, that
the covariance matrix can be defined as a polynomial of the adjacency of the
graph [19]. Since GSR is a more general model, it is less restrictive in terms
of the signals it can handle. However, it has the disadvantage of requiring a
significantly larger number of observations than GL [13].

Our proposal is to combine the advantages of assuming Gaussianity, which
implies solving a maximum likelihood problem that requires fewer node
observations, with the larger generality of graph-stationary approaches. Our
ultimate goal is to generalize the range of scenarios where GL can be used,
while keeping the number of observations and computational complexity under
control. To be more precise, we introduce Polynomial Graphical Lasso (PGL), a
new scheme to learn graphs from signals that works under the assumption that
the samples are Gaussian and graph stationary, so that the covariance
(precision) matrix of the observations can be written as a polynomial of a
sparse graph. These assumptions give rise to a constrained log-likelihood
minimization that is jointly optimized over the precision and adjacency
matrices, with GL being a particular instance of our problem. The price to pay
is that the postulated optimization, even after relaxing the sparsity
constraints, is more challenging, leading to a biconvex problem. To mitigate
this issue we provide an efficient alternating algorithm with convergence
guarantees.

Contributions. To summarize, our main contributions are:

  * •

Introducing PGL, a novel graph-learning scheme that, by assuming that the
observations are Gaussian and graph-stationary, generalizes GL and is able to
learn a meaningful graph structure in scenarios where the precision/covariance
matrices are polynomials of the sparse matrix that represents the graph.

  * •

Formulating the inference problem as a biconvex constrained optimization, with
the variables to optimize being the precision matrix and the graph. While our
focus is on learning the graphs, note that this implies that PGL can also be
used in the context of covariance estimation.

  * •

Developing an efficient algorithm to solve the proposed optimization, together
with convergence guarantees to a stationary point (block coordinatewise
minimizer).

  * •

Evaluating the performance of the proposed approach through comparisons with
alternatives from the literature on synthetic and real-world datasets.

Outline. The remainder of the paper is organized as follows. Section II
surveys the basic graph and GSP background, placing special emphasis on the
concepts of graph stationarity and Gaussian graph signals. In Section III, the
problem of learning (inferring) graphs from signals under different
assumptions on the observations is formally stated. Section IV presents a
computationally tractable relaxation of the graph-learning problem, along with
an efficient algorithm and its associated convergence guarantees. Section V
quantifies and compares the recovery performance of the proposed approach with
other methods from the literature using both synthetic and real data. Section
VI closes the paper with concluding remarks.

##  II Graphs and Random Graph Signals

This section introduces the notation used in the paper and fundamental
concepts of graphs and GSP that will help to explain the relationship between
the available signals and the topology of the underlying graph.

Notation. We represent scalar, vectors and matrices using lowercase
($x\in{\mathbb{R}}$), bold lowercase (${\mathbf{x}}\in{\mathbb{R}}^{N}$), and
bold uppercase (${\mathbf{X}}\in{\mathbb{R}}^{N\times R}$) letters,
respectively. $X_{ij}$ represents the values of the entry $[i,j]$ of matrix
${\mathbf{X}}$. ${\mathbf{I}}_{N}$ denotes the identity matrix of size
$N\times N$. The expression diag(${\mathbf{x}}$) represents a diagonal matrix
with the values of ${\mathbf{x}}$ in the diagonal. The expressions
$\|\cdot\|_{0}$ and $\|\cdot\|_{1}$ refer to the entry-wise (vectorized)
$\ell_{0}$ and $\ell_{1}$ matrix norms, and , and $\|\cdot\|_{F}$ is the
Frobenious norm.

Graphs. Let ${\mathcal{G}}=\\{{\mathcal{V}},{\mathcal{E}}\\}$ be a weighted
and undirected graph with $N$ nodes, where ${\mathcal{V}}$ and ${\mathcal{E}}$
denote the vertex and edge set, respectively. The topology of the graph
${\mathcal{G}}$ is encoded in the sparse, weighted adjacency matrix
${\mathbf{A}}\in{\mathbb{R}}^{N\times N}$, where $A_{ij}$ represents the
weight of the edge between nodes $j$ and $i$, and $A_{ij}=0$ if $i$ and $j$
are not connected. A more general representation of the graph is the Graph
Shift Operator (GSO), denoted by ${\mathbf{S}}\in\mathbb{R}^{N\times N}$,
where $S_{ij}\neq 0$ if and only if $i=j$ or $(i,j)\in{\mathcal{E}}$. Common
choices for the GSO are the adjacency matrix ${\mathbf{A}}$ [20], the
combinatorial graph Laplacian
${\mathbf{L}}:=\text{diag}({\mathbf{A}}\textbf{1})-{\mathbf{A}}$ [4], and
their degree-normalized variants. Since the graph is undirected, the GSO is
symmetric and can be diagonalized as
${\mathbf{S}}={\mathbf{V}}\boldsymbol{\Lambda}{\mathbf{V}}^{\top}$, where the
orthogonal matrix ${\mathbf{V}}\in{\mathbb{R}}^{N\times N}$ contains the
eigenvectors of the GSO, and the diagonal matrix $\boldsymbol{\Lambda}$
contains its eigenvalues.

Graph signals and filters. A graph signal is represented by a vector
${\mathbf{x}}\in\mathbb{R}^{N}$, where $x_{i}$ denotes the signal value
observed at node $i$. Alternatively, a graph signal can be understood as a
${\mathcal{V}}\rightarrow{\mathbb{R}}$ mapping, so that the signal values can
represent features associated with the node at hand. Graph filters provide a
flexible tool for either processing or modeling the relationship between a
signal ${\mathbf{x}}$ and its underlying graph ${\mathcal{G}}$. Succinctly, a
graph filter is a linear operator ${\mathbf{H}}\in{\mathbb{R}}^{N\times N}$
that takes into account the topology of the graph and is defined as a
polynomial of the GSO ${\mathbf{S}}$ of the form

${\mathbf{H}}=\sum_{l=0}^{L-1}h_{l}{\mathbf{S}}^{l}=\sum_{l=0}^{L-1}h_{l}{\mathbf{V}}\boldsymbol{\Lambda}^{l}{\mathbf{V}}^{\top}={\mathbf{V}}\Big{(}\sum_{l=0}^{L-1}h_{l}\boldsymbol{\Lambda}^{l}\Big{)}{\mathbf{V}}^{\top},$
(1)

where $\\{h_{l}\\}_{l=0}^{L-1}$ represent the filter coefficients, and
${\mathbf{V}}$ and $\boldsymbol{\Lambda}$ are the eigenvectors and eigenvalues
of the GSO respectively. Since ${\mathbf{H}}$ is a polynomial of
${\mathbf{S}}$, it readily follows that both matrices have the same
eigenvectors ${\mathbf{V}}$.

Stationary graph signals. A graph signal ${\mathbf{x}}$ is defined to be
stationary on the graph ${\mathcal{G}}$ if it can be expressed as the output
of a graph filter ${\mathbf{H}}$ excited by a zero-mean white signal
${\mathbf{w}}\in\mathbb{R}^{N}$ [19, 21, 22]. In other words, if
${\mathbf{w}}$ has a covariance of
$\mathbb{E}[{\mathbf{w}}{\mathbf{w}}^{\top}]={\mathbf{I}}_{N}$, and
${\mathbf{x}}={\mathbf{H}}{\mathbf{w}}$, then the covariance of ${\mathbf{x}}$
is
$\boldsymbol{\Sigma}=\mathbb{E}[{\mathbf{x}}{\mathbf{x}}^{\top}]={\mathbf{H}}\mathbb{E}[{\mathbf{w}}{\mathbf{w}}^{\top}]{\mathbf{H}}^{\top}={\mathbf{H}}{\mathbf{H}}^{\top}={\mathbf{H}}^{2}$.
This reveals that assuming that ${\mathbf{x}}$ is stationary on
${\mathcal{G}}$ is equivalent to saying that the covariance matrix
$\boldsymbol{\Sigma}$ can be written as a polynomial of the GSO
${\mathbf{S}}$, and that, as a result, $\boldsymbol{\Sigma}$ and
${\mathbf{S}}$ have the same eigenvectors ${\mathbf{V}}$ [cf. (1)]. Hence,
graph stationarity implies that the matrices ${\mathbf{S}}$ and
$\boldsymbol{\Sigma}$ commute, which is a significant property to be exploited
later on.

Gaussian graph signals. A Gaussian graph signal is an $N$-dimensional _random_
vector ${\mathbf{x}}$ whose (vertex-indexed) entries are jointly Gaussian.
Hence, if ${\mathbf{x}}\in{\mathbb{R}}^{N}$ is distributed as
$\mathcal{N}({\mathbf{0}},\boldsymbol{\Sigma})$, we have that both the
covariance matrix $\boldsymbol{\Sigma}$ and the precision matrix
$\boldsymbol{\Theta}=\boldsymbol{\Sigma}^{-1}$ are $N\times N$ matrices whose
entries capture statistical relationships between the nodal values. In
particular, using the distribution of the Gaussian, it holds that:

$f_{\boldsymbol{\Theta}}({\mathbf{x}})=(2\pi)^{-N/2}\cdot\det(\boldsymbol{\Theta})^{\frac{1}{2}}\cdot\exp^{-\frac{1}{2}{\mathbf{x}}^{T}\boldsymbol{\Theta}{\mathbf{x}}},$
(2)

demonstrating that, in the context of graph signals, the entries of
$\boldsymbol{\Theta}\in{\mathbb{R}}^{N\times N}$ account for conditional
dependence relationships between the (features of the) nodes of the graph [7].

Suppose now that we have a collection of $R$ Gaussian signals
${\mathbf{X}}=[{\mathbf{x}}_{1},...,{\mathbf{x}}_{R}]\in\mathbb{R}^{N\times
R}$ each of them independently drawn from the distribution in (2). The log-
likelihood associated with ${\mathbf{X}}$ is

$L({\mathbf{X}}|\boldsymbol{\Theta})=\prod_{r=1}^{R}f_{\boldsymbol{\Theta}}({\mathbf{x}}_{r}),\
\ \
{\mathcal{L}}({\mathbf{X}}|\boldsymbol{\Theta})=\sum_{r=1}^{R}\log(f_{\boldsymbol{\Theta}}({\mathbf{x}}_{r})).$
(3)

This expression will be exploited when formulating our proposed inference
approach and establishing links with classical methods.

##  III Graph learning problem formulation

This section begins with a formal definition of the graph learning problem,
followed by an explanation of some common approaches used in the literature to
tackle this problem. Afterwards, we proceed to formalize the learning problem
we aim to address and cast it as an optimization problem. We then provide an
overview of the key features of our formulation and conduct a comparative
analysis with the two closest approaches available in the literature.

To formally state the graph learning problem, let us recall that we assume: i)
${\mathcal{G}}$ is an undirected graph with $N$ nodes, ii) there is a random
process associated with ${\mathcal{G}}$, and iii) we denote by
${\mathbf{X}}=[{\mathbf{x}}_{1},...,{\mathbf{x}}_{R}]\in{\mathbb{R}}^{N\times
R}$ a collection of $R$ independent realizations of such a process. The goal
in graph learning is to use a given set of nodal observations to find the
(estimates of the) links/associations between the nodes in the graph (i.e.,
use ${\mathbf{X}}$ to estimate the ${\mathbf{S}}$ associated with
${\mathcal{G}}$).

This problem has been addressed under different approaches [23, 24, 19].
Differences among these models typically arise from the underlying assumptions
that are made regarding 1) the graph, which almost universally entails just
considering that the graph is sparse (see, e.g., [25] for a recent exception),
2) the signals, which assume certain properties related to the nature of the
signals such as smoothness [24] or Gaussianity, [23] among others, and 3) the
relationship between the graph and the signal, such as the stationarity
property [19].

The model we propose aims to incorporate several assumptions about both the
graph and the graph signals. To that end, we propose an approach for which we
assume that: 1) the graph is sparse, 2) the signals are Gaussian and, 3) the
signals are stationary in the underlying graph.

Having established these assumptions, we now proceed to formally state our
graph learning problem as follows

Problem 1 Given a set of signals ${\mathbf{X}}\in{\mathbb{R}}^{N\times R}$,
find the underlying sparse graph structure encoded in ${\mathbf{S}}$ under the
assumptions:

(AS1): The graph signals ${\mathbf{X}}$ are i.i.d. realizations of
${\mathcal{N}}({\mathbf{0}},\boldsymbol{\Sigma})$.

(AS2): The graph signals ${\mathbf{X}}$ are stationary in ${\mathbf{S}}$.

Our approach is to recast Problem 1 as the following optimization

$\displaystyle\\!\underset{\boldsymbol{\Theta}\succeq
0,{\mathbf{S}}\in{\mathcal{S}}}{\mathsf{minimize}}\ $
$\displaystyle-\log(\det(\boldsymbol{\Theta}))+\text{tr}(\hat{\boldsymbol{\Sigma}}\boldsymbol{\Theta})+\rho\|{\mathbf{S}}\|_{0}.$
$\displaystyle\\!\mathsf{subject~{}to}$
$\displaystyle\;\;\boldsymbol{\Theta}{\mathbf{S}}={\mathbf{S}}\boldsymbol{\Theta},$
(4)

where ${\mathcal{S}}$ is a generic set representing additional constraints
that ${\mathbf{S}}$ is known to satisfy (e.g., the GSO being symmetric and its
entries being between zero and one). The minimization takes as input the
sample covariance matrix
$\hat{\boldsymbol{\Sigma}}=\frac{1}{R}{\mathbf{X}}{\mathbf{X}}^{\mathsf{T}}$
and generates as output the estimate for ${\mathbf{S}}$ and, as a byproduct,
the estimate for $\boldsymbol{\Theta}$. For the problem in (III), we require
$\boldsymbol{\Theta}$ to be positive semidefinite. This constraint arises
because $\boldsymbol{\Theta}$ is the inverse of $\hat{\boldsymbol{\Sigma}}$
which is symmetric and positive definite by construction. Consequently,
$\boldsymbol{\Theta}$ inherits the properties of being symmetric and positive
semidefinite.

Next, we explain the motivation for each term in (III) with special emphasis
on the constraint
$\boldsymbol{\Theta}{\mathbf{S}}={\mathbf{S}}\boldsymbol{\Theta}$, which is a
fundamental component of our approach.

  * •

The first two terms in the objective function are due to (AS1) and arise from
minimizing the negative log-likelihood expression in (3). Indeed, it is clear
that substituting (2) into (3) yields

$\displaystyle\sum_{r=1}^{R}\left(-\frac{N}{2}\log(2\pi)-\log(\det(\boldsymbol{\Theta}))+\text{tr}({\mathbf{x}}_{r}^{T}\boldsymbol{\Theta}{\mathbf{x}}_{r})\right).$

Since constants are irrelevant for the optimization, we drop the first term
and divide the other two by $R$, yielding

$\displaystyle-\log(\det(\boldsymbol{\Theta}))+\frac{1}{R}\sum_{r=1}^{R}\text{tr}({\mathbf{x}}_{r}{\mathbf{x}}_{r}^{T}\boldsymbol{\Theta})$
$\displaystyle=-\log(\det(\boldsymbol{\Theta}))+\text{tr}(\hat{\boldsymbol{\Sigma}}\boldsymbol{\Theta}).$
(5)

  * •

The term $\rho\|{\mathbf{S}}\|_{0}$ accounts for the fact of ${\mathbf{S}}$
being a GSO (hence, sparse), with $\rho>0$ being a regularization parameter
that determines the desired level of sparsity in the graph.

  * •

Finally, the equality constraint serves to embody (AS2). It is important to
highlight that the polynomial relationship between $\boldsymbol{\Sigma}$ and
${\mathbf{S}}$, as implied by (AS2), is typically addressed in estimation and
optimization problems through either: i) extracting the eigenvectors of
$\boldsymbol{\Sigma}$ and enforcing them to be the eigenvectors of
${\mathbf{S}}$ [13], or ii) imposing the constraint
$\boldsymbol{\Sigma}{\mathbf{S}}={\mathbf{S}}\boldsymbol{\Sigma}$ [16]. In
contrast, our approach encodes the polynomial relation implied by (AS2) by
enforcing commutativity between $\boldsymbol{\Theta}$ and ${\mathbf{S}}$. Note
that if $\boldsymbol{\Sigma}$ and ${\mathbf{S}}$ are full-rank matrices and
commute, it follows that $\boldsymbol{\Theta}$ and ${\mathbf{S}}$ also
commute. In other words, $\boldsymbol{\Theta}=\boldsymbol{\Sigma}^{-1}$ can be
represented as a polynomial in ${\mathbf{S}}$, which can be verified by the
Cayley-Hamilton theorem.

It is important to note that the assumption of stationarity may seem stringent
since it implies commutativity between ${\mathbf{S}}$ and
$\boldsymbol{\Theta}$. However, it provides more degrees of freedom than many
existing methods. For example, in partial correlation methods, ${\mathbf{S}}$
is restricted to be ${\mathbf{S}}=\boldsymbol{\Theta}$, while in our case,
assuming stationarity allows $\boldsymbol{\Theta}$ to be _any_ polynomial in
${\mathbf{S}}$. This leads to a more general approach, including partial
correlation as a particular case. To further illustrate this point, consider
the sparse structural equation model
${\mathbf{X}}={\mathbf{A}}{\mathbf{X}}+{\mathbf{W}}$, with ${\mathbf{W}}$
being white noise. GL will identify
${\mathbf{S}}=\boldsymbol{\Theta}=({\mathbf{I}}_{N}-{\mathbf{A}})^{2}$, while
PGL will identify ${\mathbf{S}}={\mathbf{A}}$.

To better understand the features of PGL, let us briefly discuss the main
differences relative to its two closest competitors: GSR and GL.

GSR handles Problem 1 without considering the Gaussian assumption in (AS1). As
a result, the first two terms in (III), which are associated with the log-
likelihood function, are not present. This reduces the problem to inferring a
sparse graph under the stationarity constraint. The stationarity assumption is
incorporated into (III) through the expression
$\boldsymbol{\Theta}{\mathbf{S}}={\mathbf{S}}\boldsymbol{\Theta}$, which is
equivalent to
$\boldsymbol{\Sigma}{\mathbf{S}}={\mathbf{S}}\boldsymbol{\Sigma}$ if
$\boldsymbol{\Sigma}$ is a full-rank matrix. This property enables learning
the graph by solving the following optimization problem with the commutativity
constraint between $\hat{\boldsymbol{\Sigma}}$ and ${\mathbf{S}}$:

$\displaystyle\\!\underset{{\mathbf{S}}\in{\mathcal{S}}}{\mathsf{minimize}}\quad\rho\|{\mathbf{S}}\|_{0}\;\;\;\mathsf{subject~{}to}\;\;\;\hat{\boldsymbol{\Sigma}}{\mathbf{S}}={\mathbf{S}}\hat{\boldsymbol{\Sigma}},$
(6)

where the constraint is typically relaxed as
$\|\hat{\boldsymbol{\Sigma}}{\mathbf{S}}-{\mathbf{S}}\hat{\boldsymbol{\Sigma}}\|_{F}^{2}\leq\epsilon$
to account for the fact that we have
$\hat{\boldsymbol{\Sigma}}\approx\boldsymbol{\Sigma}$. By assuming
stationarity, the formulation in (6) allows the sample covariance to be
modeled as any polynomial in ${\mathbf{S}}$, making it more general than the
formulation in (III). However, the absence of Gaussianity in (6) means that it
is no longer a maximum likelihood estimation. As a result, correctly
identifying the ground truth ${\mathbf{S}}$ requires very reliable estimates
of $\hat{\boldsymbol{\Sigma}}$, which usually entails having access to a large
number of signals to set $\epsilon$ close to zero. This is indeed a challenge,
especially in setups with a large number of nodes.

For the second scenario, suppose we simplify (AS2) and instead of considering
$\boldsymbol{\Theta}$ as any polynomial in ${\mathbf{S}}$, we restrict it to a
particular case with the following structure:
$\boldsymbol{\Theta}=\boldsymbol{\Sigma}^{-1}=\sigma{\mathbf{I}}+\delta{\mathbf{S}}$.
Then, up to the diagonal values and scaling issues, the sparse matrix
${\mathbf{S}}$ to be estimated and
$\boldsymbol{\Theta}=\boldsymbol{\Sigma}^{-1}$ are the same. Consequently, it
suffices to optimize over one of them, leading to the well-known GL
formulation:

$\displaystyle\\!\underset{\boldsymbol{\Theta}\succeq
0,\boldsymbol{\Theta}\in{\mathcal{S}}_{\boldsymbol{\Theta}}}{\mathsf{minimize}}\
$
$\displaystyle-\log(\det(\boldsymbol{\Theta}))+\text{tr}(\hat{\boldsymbol{\Sigma}}\boldsymbol{\Theta})+\rho\|\boldsymbol{\Theta}\|_{0}.$
(7)

The main advantages of (7) relative to (III) are that the number of variables
is smaller and the resulting problem (after relaxing the $\ell_{0}$ norm) is
convex. The main drawback is that by forcing the support of ${\mathbf{S}}$ and
$\boldsymbol{\Theta}$ to be the same, the set of feasible graphs (and their
topological properties) is more limited. Indeed, GL can only estimate graphs
that are positive definite, while the problem in (III) can yield any symmetric
matrix. Remarkably, when the model assumed in (7) holds true (i.e., data is
Gaussian and $\boldsymbol{\Theta}$ is sparse), GL is able to find reliable
estimates of ${\mathbf{S}}$ even when the number of samples $R$ is fairly low.
On the other hand, simulations will show that GL does a poor job estimating
${\mathbf{S}}$ when the relation between the precision matrix and
${\mathcal{G}}$ is more involved.

In conclusion, from a conceptual point of view, our formulation reaches a
favorable balance between GL and graph-stationarity approaches. This leads to
the following two main advantages i) a more general model than GL since our
approach models $\boldsymbol{\Theta}$ as any polynomial in ${\mathbf{S}}$ and
ii) a model with more structure than the graph-stationarity approaches due to
the incorporation of (AS1). However, it is important to note that the
optimization in (III), even if the $\ell_{0}$ norm is relaxed, lacks convexity
due to the presence of a bilinear constraint that couples the optimization
variables $\boldsymbol{\Theta}$ and ${\mathbf{S}}$. These challenges will be
addressed in the subsequent section.

##  IV Biconvex relation and algorithm design

As explained in the previous section, the problem in (III) is not convex and
this challenges designing an algorithm to find a good solution. This section
reformulates (III), develops an iterative algorithm, referred to as PGL, to
estimate ${\mathbf{S}}$ and $\boldsymbol{\Theta}$, and characterizes its
convergence to a coordinate-wise minimum point. The proposed approach involves
several modifications: 1) we replace the $\ell_{0}$-norm with an elastic net
regularizer, which is convex [26]; and 2) we relax the commutativity
constraint using an inequality instead of an equality. Next, we explain step
by step the resulting formulation.

###  IV-A Biconvex relaxation

The first modification to reformulate (III) is to relax the constraint that
imposes commutativity between ${\mathbf{S}}$ and $\boldsymbol{\Theta}$. Such a
constraint is stringent and significantly narrows the feasible solution set of
(III), which may not be practical in real-world scenarios. Furthermore,
considering that our access to the covariance (or precision) matrix is limited
to its sampled estimates, enforcing exact commutativity is excessively
restrictive. To mitigate this, we relax the original constraint by replacing
the matrix equality
$\boldsymbol{\Theta}{\mathbf{S}}={\mathbf{S}}\boldsymbol{\Theta}$ with the
scalar Frobenius norm-based inequality
$\|\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{S}}\boldsymbol{\Theta}\|_{F}^{2}\leq\delta$.
This modification not only expands the feasible region but also endows the
model with a greater degree of robustness.

The second modification addresses the non-convexity of the objective in (III),
which originates from the use of the $\ell_{0}$-norm. To alleviate this issue,
we relax the problem using an elastic net regularizer. Specifically, we
replace the $\rho\|{\mathbf{S}}\|_{0}$ penalty with
$\rho\big{(}\|{\mathbf{S}}\|_{1}+\frac{\eta}{2\rho}\|{\mathbf{S}}\|_{F}^{2}\big{)}$,
where the parameter $\eta$ controls the trade-off between the $\ell_{1}$-norm
and the Frobenius norm components, and is typically set to a very small value.
Although elastic net regularizers have demonstrated practical effectiveness,
alternative methods for relaxing the $\ell_{0}$-norm exist (see, for example,
[27, 28]), each offering distinct trade-offs in computational complexity,
convergence speed, and theoretical underpinnings.

With the incorporation of these two modifications, we reformulate the original
graph learning problem presented in (III) as follows

$\begin{array}[]{rll}&\underset{\boldsymbol{\Theta}\succeq
0,{\mathbf{S}}\in{\mathcal{S}}}{\mathsf{minimize}}&-\log(\det(\boldsymbol{\Theta}))+\text{tr}(\hat{\boldsymbol{\Sigma}}\boldsymbol{\Theta})+\rho\|{\mathbf{S}}\|_{1}+\frac{\eta}{2}\|{\mathbf{S}}\|_{F}^{2},\\\
&\mathsf{subject~{}to}&\|\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{S}}\boldsymbol{\Theta}\|_{F}\leq\delta,\end{array}$
(8)

In this setup, $\delta$ serves as a hyperparameter chosen according to the
quality of the estimation of $\hat{\boldsymbol{\Sigma}}$ which affects the
estimation of $\boldsymbol{\Theta}$. A smaller value of $\delta$ is
appropriate when the quality of $\hat{\boldsymbol{\Sigma}}$ is high, which
typically corresponds to having a sample size $R$ that is substantially larger
than the number of nodes. While the relaxation of the commutativity constraint
enhances the robustness of our formulation to data quality and simplifies the
optimization by reducing the number of Lagrange multipliers, the product of
$\boldsymbol{\Theta}$ and ${\mathbf{S}}$ still introduces nonconvexity into
the problem. The way we propose for dealing with the (updated) biconvex
constraint is to solve (8) using an alternating optimization algorithm. This
family of algorithms is widely used to approximate nonconvex problems by
dividing the original problem into several convex subproblems and solving them
with respect to each of the variables by fixing all the others. In our
particular case, this methodology involves alternately solving for
$\boldsymbol{\Theta}$ with ${\mathbf{S}}$ held fixed, and then updating
${\mathbf{S}}$ using the newly updated $\boldsymbol{\Theta}$, at each
iteration.

In the subsequent two subsections, we delve into the detailed methodologies
employed to solve each of the two subproblems. Following this, we outline the
overall algorithm and discuss its convergence properties. To simplify
exposition, in the remainder of the section, we will assume that
${\mathbf{S}}$ represents the _adjacency_ matrix of an undirected graph.
Consequently, the feasible solution set for ${\mathbf{S}}$ is defined as:

${\mathcal{S}}\\!:=\\!\big{\\{}{\mathbf{S}}\\!\in\\!{\mathbb{R}}^{N\times
N}|{\mathbf{S}}\\!=\\!{\mathbf{S}}^{T};\ {\mathbf{S}}\\!\geq\\!\mathbf{0};\
\mathrm{diag}({\mathbf{S}})\\!=\\!{\mathbf{0}};\
{\mathbf{S}}{\mathbf{1}}\\!\geq\\!{\mathbf{1}}\big{\\}},$ (9)

where ${\mathbf{S}}$ is constrained to be symmetric with zero diagonal entries
and non-negative off-diagonal elements. The additional condition
${\mathbf{S}}\mathbf{1}\geq\mathbf{1}$ is imposed to preclude the trivial
solution, i.e., ${\mathbf{S}}={\mathbf{0}}$. Nonetheless, the techniques
presented next can be readily extended to accommodate different forms of
${\mathbf{S}}$.

###  IV-B Solving subproblem for ${\mathbf{S}}$

We begin by addressing the subproblem with respect to ${\mathbf{S}}$, while
holding $\boldsymbol{\Theta}$ fixed. The subproblem is formulated as:

$\begin{array}[]{rll}&\underset{{\mathbf{S}}\in{\mathcal{S}}}{\mathsf{minimize}}&\rho\|{\mathbf{S}}\|_{1}+\frac{\eta}{2}\|{\mathbf{S}}\|_{F}^{2},\\\
&\mathsf{subject~{}to}&\|\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{S}}\boldsymbol{\Theta}\|_{F}\leq\delta.\end{array}$
(10)

To solve (10), we adopt a linearized ADMM approach, which introduces an
auxiliary variable ${\mathbf{P}}$ and leads to the following equivalent
formulation:

$\begin{array}[]{rll}&\underset{{\mathbf{S}}\in{\mathcal{S}},{\mathbf{P}}}{\mathsf{minimize}}&\rho\|{\mathbf{S}}\|_{1}+\frac{\eta}{2}\|{\mathbf{S}}\|_{F}^{2},\\\
&\mathsf{subject~{}to}&\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{S}}\boldsymbol{\Theta}={\mathbf{P}},\|{\mathbf{P}}\|_{F}\leq\delta.\end{array}$
(11)

The augmented Lagrangian associated with (11) is then given by

$\displaystyle\\!L({\mathbf{S}},{\mathbf{P}},{\mathbf{Z}})=\rho\|{\mathbf{S}}\|_{1}+$
$\displaystyle\frac{\eta}{2}\|{\mathbf{S}}\|_{F}^{2}+\langle{\mathbf{Z}},\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{S}}\boldsymbol{\Theta}-{\mathbf{P}}\rangle$
$\displaystyle+\frac{\beta}{2}\|\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{S}}\boldsymbol{\Theta}-{\mathbf{P}}\|_{F}^{2},$
(12)

where ${\mathbf{Z}}$ is the Lagrange multiplier.

To update ${\mathbf{S}}$ at the $t$-th iteration, we address the following
minimization problem:

$\underset{{\mathbf{S}}\in{\mathcal{S}}}{\mathsf{minimize}}\quad\rho\|{\mathbf{S}}\|_{1}+\frac{\eta}{2}\|{\mathbf{S}}\|_{F}^{2}+\frac{\beta}{2}\|\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{S}}\boldsymbol{\Theta}-{\mathbf{P}}+\frac{1}{\beta}{\mathbf{Z}}\|_{F}^{2},$
(13)

where, for the sake of simplicity, we omit the iteration subscript from
$\boldsymbol{\Theta}^{(t)}$ and ${\mathbf{Z}}^{(t)}$. Problem (13) does not
admit a closed-form solution due to the term
$\frac{1}{2}\|\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{S}}\boldsymbol{\Theta}-{\mathbf{P}}+\frac{1}{\beta}{\mathbf{Z}}\|_{F}^{2}$.
To deal conveniently with this problem we resort to the majorization-
minimization (MM) algorithm [29]. We denote this term as $g({\mathbf{S}})$ and
proceed to majorize both $g({\mathbf{S}})$ and
$\frac{\eta}{2}\|{\mathbf{S}}\|_{F}^{2}$ at the point ${\mathbf{S}}^{(t)}$,
resulting in the following problem:

$\begin{split}\underset{{\mathbf{S}}\in{\mathcal{S}}}{\mathsf{minimize}}\quad&\langle\rho{\mathbf{I}}_{N\times
N}+\eta{\mathbf{S}}^{(t)}+\beta\nabla
g({\mathbf{S}}^{(t)}),{\mathbf{S}}-{\mathbf{S}}^{(t)}\rangle\\\
&\qquad+\frac{L_{1}}{2}\|{\mathbf{S}}-{\mathbf{S}}^{(t)}\|_{F}^{2},\end{split}$
(14)

where $\nabla g({\mathbf{S}})$ represents the gradient of $g({\mathbf{S}})$,
detailed in the equation:

$\displaystyle\\!\nabla g({\mathbf{S}})=\;\;$
$\displaystyle\boldsymbol{\Theta}\boldsymbol{\Theta}{\mathbf{S}}+{\mathbf{S}}\boldsymbol{\Theta}\boldsymbol{\Theta}-2\boldsymbol{\Theta}{\mathbf{S}}\boldsymbol{\Theta}+{\mathbf{P}}\boldsymbol{\Theta}$
$\displaystyle-\boldsymbol{\Theta}{\mathbf{P}}+\frac{1}{\beta}(\boldsymbol{\Theta}{\mathbf{Z}}-{\mathbf{Z}}\boldsymbol{\Theta}).$
(15)

Now, Problem (14) has a closed-form solution, allowing for the update of
${\mathbf{S}}^{(t+1)}$ as follows:

${\mathbf{S}}^{(t+1)}={\mathcal{P}}_{{\mathcal{S}}}\Big{(}{\mathbf{S}}^{(t)}-\frac{1}{L_{1}}\big{(}\rho{\mathbf{I}}_{N}+\eta{\mathbf{S}}^{(t)}+\beta\nabla
g\bigl{(}{\mathbf{S}}^{(t)}\bigl{)}\big{)}\Big{)},$ (16)

where ${\mathcal{P}}_{{\mathcal{S}}}$ is the projection onto the set
$\mathcal{S}$ with respect to the Frobenius norm, which can be computed
efficiently by the Dykstra’s projection algorithm [30]. More specifically, the
set $\mathcal{S}$ can be written as the intersection of two closed convex sets
as follows:

$\mathcal{S}=\mathcal{S}_{A}\cap\mathcal{S}_{B},$ (17)

where $\mathcal{S}_{A}:=\left\\{\mathbf{S}\in\mathbb{R}^{N\times
N}\,|\,\mathbf{S}=\mathbf{S}^{T}\right\\}$ and
$\mathcal{S}_{B}:=\left\\{\mathbf{S}\in\mathbb{R}^{N\times
N}\,|\,\mathbf{S}\geq
0;\,\mathrm{diag}(\mathbf{S})=\mathbf{0};\,\mathbf{S}\mathbf{1}\geq\mathbf{1}\right\\}$.
We employ Dykstra’s projection algorithm [30] to compute the nearest point
projection of a given point onto the intersection of sets $\mathcal{S}_{A}$
and $\mathcal{S}_{B}$. Dykstra’s algorithm achieves this by alternately
projecting the point onto $\mathcal{S}_{A}$ and $\mathcal{S}_{B}$ until the
solution is reached. For a more comprehensive understanding of Dykstra’s
projection algorithm, the reader is directed to [30]. Detailed descriptions of
the projection computations onto sets $\mathcal{S}_{A}$ and $\mathcal{S}_{B}$
are provided in Appendix A.

Input: $\hat{\boldsymbol{\Theta}}^{(k)}$, ${\hat{\mathbf{S}}}^{(k)}$,
${\hat{\mathbf{P}}}^{(k)}$, ${\hat{\mathbf{Z}}}^{(k)}$, $\rho$, $\eta$,
$\beta$, $\delta$

Outputs : ${\hat{\mathbf{S}}}^{(k+1)}$, ${\hat{\mathbf{P}}}^{(k+1)}$,
${\hat{\mathbf{Z}}}^{(k+1)}$

1

2Initialize ${\mathbf{S}}^{(0)}\\!=\\!{\hat{\mathbf{S}}}^{(k)}$,
${\mathbf{P}}^{(0)}\\!=\\!{\mathbf{P}}^{(k)}$,
${\mathbf{Z}}^{(0)}\\!=\\!{\mathbf{Z}}^{(k)}$

3

4for  _$t=0$ to $T-1$_ do

5 Update ${\mathbf{S}}^{(t+1)}$ by (16);

6 Update ${\mathbf{P}}^{(t+1)}$ by (19);

7 Update ${\mathbf{Z}}^{(t+1)}$ by (23);

8

9 end for

10 ${\hat{\mathbf{S}}}^{(k+1)}={\mathbf{S}}^{(T)}$,
${\hat{\mathbf{P}}}^{(k+1)}={\mathbf{P}}^{(T)}$,
${\hat{\mathbf{Z}}}^{(k+1)}={\mathbf{Z}}^{(T)}$.

11

12

Algorithm 1 Inner loop for ${\mathbf{S}}$ update.

Returning to the augmented Lagrangian in (IV-B), we update ${\mathbf{P}}$ at
the $t$-th iteration by solving the following problem:

$\begin{array}[]{rll}&\underset{{\mathbf{P}}}{\mathsf{minimize}}&\frac{\beta}{2}\|{\mathbf{P}}-\boldsymbol{\Theta}{\mathbf{S}}+{\mathbf{S}}\boldsymbol{\Theta}-\frac{1}{\beta}{\mathbf{Z}}\|_{F}^{2},\\\
&\mathsf{subject~{}to}&\|{\mathbf{P}}\|_{F}\leq\delta,\end{array}$ (18)

where we have simplified the notation by omitting the iteration subscripts
from ${\mathbf{S}}^{(t+1)}$ and ${\mathbf{Z}}^{(t)}$. Problem (18) has a
closed-form solution. As a result, ${\mathbf{P}}^{(t+1)}$ can be updated by

${\mathbf{P}}^{(t+1)}={\mathcal{P}}_{\delta}\left(\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{S}}\boldsymbol{\Theta}+\frac{1}{\beta}{\mathbf{Z}}\right),$
(19)

where ${\mathcal{P}}_{\delta}$ denotes the projection defined by:

$\displaystyle{\mathcal{P}}_{\delta}({\mathbf{A}})=\left\\{\begin{array}[]{ll}\frac{\delta}{\|{\mathbf{A}}\|_{F}}{\mathbf{A}}&\mbox{if
$\|{\mathbf{A}}\|_{F}>\delta$}\\\
{\mathbf{A}}&\mbox{otherwise.}\end{array}\right.$ (22)

Finally, the dual variable ${\mathbf{Z}}$ is updated according to:

${\mathbf{Z}}^{(t+1)}={\mathbf{Z}}^{(t)}+\beta(\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{S}}\boldsymbol{\Theta}-{\mathbf{P}}),$
(23)

where the iteration subscripts from ${\mathbf{S}}^{(t+1)}$ and
${\mathbf{P}}^{(t+1)}$ have been omitted for simplicity. A pseudocode of the
steps to be performed for the update of ${\mathbf{S}}$ is summarized in
Algorithm 1.

If the parameter $L_{1}$ in (14) is larger than the Lipschitz constant of the
gradient of $\beta
g\bigl{(}{\mathbf{S}}\bigl{)}+\frac{\eta}{2}\|{\mathbf{S}}\|_{F}^{2}$, then
the sequence
$\bigl{\\{}\big{(}{\mathbf{S}}^{(t)},{\mathbf{P}}^{(t)}\big{)}\bigl{\\}}$
converges to the optimal solution of Problem (11), and
$\bigl{\\{}{\mathbf{Z}}^{(t)}\bigl{\\}}$ converges to the optimal solution of
the dual of problem (11), which follows from the existing convergence result
of majorized ADMM [31]. To enhance empirical convergence rates, adopting a
more proactive strategy for selecting the parameter $L_{1}$ is beneficial. For
example, utilizing a backtracking line search to determine the stepsize in
(16) can help to accelerate convergence.

We note that the choice of the penalty parameter $\beta$ can affect the
convergence speed of the ADMM algorithm. A poorly chosen $\beta$ may lead to
very slow convergence in practice. Adaptive schemes that adjust $\beta$ have
been shown to often result in better practical performance. For example, we
can adopt the adaptive update rule presented in [32]:

$\displaystyle\beta^{(t+1)}=\left\\{\begin{array}[]{ll}\tau^{\mathrm{inc}}\beta^{(t)},&\text{if
}\big{\|}\mathbf{r}^{(t)}\big{\|}_{F}>\mu\big{\|}\mathbf{s}^{(t)}\big{\|}_{F},\\\
\beta^{(t)}/\tau^{\mathrm{dec}},&\text{if
}\big{\|}\mathbf{s}^{(t)}\big{\|}_{F}>\mu\big{\|}\mathbf{r}^{(t)}\big{\|}_{F},\\\
\beta^{(t)},&\text{otherwise},\end{array}\right.$ (27)

where $\mu>1$, $\tau^{\mathrm{inc}}>1$, and $\tau^{\mathrm{dec}}>1$ are
predefined parameters. Here, $\mathbf{r}^{(t)}$ and $\mathbf{s}^{(t)}$
represent the primal and dual residuals at iteration $t$, respectively. They
are defined as

$\mathbf{r}^{(t)}=\boldsymbol{\Theta}{\mathbf{S}}^{(t)}-{\mathbf{S}}^{(t)}\boldsymbol{\Theta}-\mathbf{P}^{(t)},$

and

$\mathbf{s}^{(t)}=\beta^{(t)}\boldsymbol{\Theta}\big{(}\mathbf{P}^{(t)}-\mathbf{P}^{(t-1)}\big{)}-\beta^{(t)}\big{(}\mathbf{P}^{(t)}-\mathbf{P}^{(t-1)}\big{)}\boldsymbol{\Theta}.$

Although it can be challenging to prove the convergence of ADMM when $\beta$
varies by iteration, the convergence theory established for a fixed $\beta$
remains applicable if one assumes that $\beta$ becomes constant after a finite
number of iterations.

###  IV-C Solving subproblem for $\boldsymbol{\Theta}$

Using the formulation from (8), we now turn our attention to the subproblem
for $\boldsymbol{\Theta}$

$\begin{array}[]{rll}&\underset{\boldsymbol{\Theta}\succeq
0}{\mathsf{minimize}}&-\log(\det(\boldsymbol{\Theta}))+\text{tr}(\hat{\boldsymbol{\Sigma}}\boldsymbol{\Theta}),\\\
&\mathsf{subject~{}to}&\|{\mathbf{S}}\boldsymbol{\Theta}-\boldsymbol{\Theta}{\mathbf{S}}\|_{F}\leq\delta.\end{array}$
(28)

Similarly to the approach taken for the ${\mathbf{S}}$ subproblem, we
reformulate the subproblem (28) for $\boldsymbol{\Theta}$ as follows

$\begin{array}[]{rll}&\underset{\boldsymbol{\Theta}\succeq 0,\
{\mathbf{Q}}}{\mathsf{minimize}}&-\log(\det(\boldsymbol{\Theta}))+\text{tr}(\hat{\boldsymbol{\Sigma}}\boldsymbol{\Theta}),\\\
&\mathsf{subject~{}to}&{\mathbf{S}}\boldsymbol{\Theta}-\boldsymbol{\Theta}{\mathbf{S}}={\mathbf{Q}},\
\|{\mathbf{Q}}\|_{F}\leq\delta.\end{array}$ (29)

The augmented Lagrangian associated with (29) is given by

$\displaystyle L(\boldsymbol{\Theta},{\mathbf{Q}},{\mathbf{Y}})=$
$\displaystyle-\log(\det(\boldsymbol{\Theta}))+\text{tr}(\hat{\boldsymbol{\Sigma}}\boldsymbol{\Theta})$
(30)
$\displaystyle+\frac{\beta}{2}\|{\mathbf{S}}\boldsymbol{\Theta}-\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{Q}}+\frac{1}{\beta}{\mathbf{Y}}\|_{F}^{2}.$

To update $\boldsymbol{\Theta}$ at the $t$-th iteration, we solve the
following optimization problem

$\displaystyle\underset{\boldsymbol{\Theta}\succeq 0}{\mathsf{minimize}}\ $
$\displaystyle-\log(\det(\boldsymbol{\Theta}))+\text{tr}(\hat{\boldsymbol{\Sigma}}\boldsymbol{\Theta})$
$\displaystyle\qquad+\frac{\beta}{2}\|{\mathbf{S}}\boldsymbol{\Theta}-\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{Q}}+\frac{1}{\beta}{\mathbf{Y}}\|_{F}^{2},$
(31)

where we have omitted the iteration subscripts from ${\mathbf{Q}}^{(t)}$ and
${\mathbf{Y}}^{(t)}$ for simplicity.

Let
$f(\boldsymbol{\Theta})=\frac{1}{2}\|{\mathbf{S}}\boldsymbol{\Theta}-\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{Q}}+\frac{1}{\beta}{\mathbf{Y}}\|_{F}^{2}$.
We then construct the majorizer of the objective function in (IV-C) at the
point $\boldsymbol{\Theta}^{(t)}$ and obtain

$\displaystyle\underset{\boldsymbol{\Theta}\succeq 0}{\mathsf{minimize}}\ $
$\displaystyle-\log(\det(\boldsymbol{\Theta}))+\langle\beta\nabla
f(\boldsymbol{\Theta}^{(t)})+\hat{\boldsymbol{\Sigma}},\boldsymbol{\Theta}-\boldsymbol{\Theta}^{(t)}\rangle$
$\displaystyle\qquad+\frac{L_{2}}{2}\|\boldsymbol{\Theta}-\boldsymbol{\Theta}^{(t)}\|_{F}^{2},$
(32)

where $\nabla f(\boldsymbol{\Theta})$ denotes the gradient of
$f(\boldsymbol{\Theta})$

$\displaystyle\nabla f(\boldsymbol{\Theta})$
$\displaystyle={\mathbf{S}}{\mathbf{S}}\boldsymbol{\Theta}+\boldsymbol{\Theta}{\mathbf{S}}{\mathbf{S}}-2{\mathbf{S}}\boldsymbol{\Theta}{\mathbf{S}}$
$\displaystyle\quad+{\mathbf{Q}}{\mathbf{S}}-{\mathbf{S}}{\mathbf{Q}}+\frac{1}{\beta}({\mathbf{S}}{\mathbf{Y}}-{\mathbf{Y}}{\mathbf{S}}).$
(33)

Input: $\boldsymbol{\Sigma}$, $\hat{\boldsymbol{\Theta}}^{(k)}$,
${\hat{\mathbf{S}}}^{(k+1)}$, ${\hat{\mathbf{Q}}}^{(k)}$,
${\hat{\mathbf{Y}}}^{(k)}$, $\beta$, $\delta$

Outputs : $\hat{\boldsymbol{\Theta}}^{(k+1)}$, ${\hat{\mathbf{Q}}}^{(k+1)}$,
${\hat{\mathbf{Y}}}^{(k+1)}$

1

2Initialize $\boldsymbol{\Theta}^{(0)}=\hat{\boldsymbol{\Theta}}^{(k)}$,
${\mathbf{Q}}^{(0)}={\hat{\mathbf{Q}}}^{(k)}$,
${\mathbf{Y}}^{(0)}={\hat{\mathbf{Y}}}^{(k)}$

3

4for  _$t=0$ to $T-1$_ do

5 Update $\boldsymbol{\Theta}^{(t+1)}$ by (34);

6 Update ${\mathbf{Q}}^{(t+1)}$ by (36);

7 Update ${\mathbf{Y}}^{(t+1)}$ by (37);

8

9 end for

10 $\hat{\boldsymbol{\Theta}}^{(k+1)}=\boldsymbol{\Theta}^{(T)}$,
${\hat{\mathbf{Q}}}^{(k+1)}={\mathbf{Q}}^{(T)}$,
${\hat{\mathbf{Y}}}^{(k+1)}={\mathbf{T}}^{(T)}$.

11

Algorithm 2 Inner loop for $\boldsymbol{\Theta}$ update.

Lemma 1. The optimal solution of problem (IV-C) is [33]

$\displaystyle\boldsymbol{\Theta}^{(t+1)}={\mathbf{U}}\left(\frac{\boldsymbol{\Lambda}+\sqrt{\boldsymbol{\Lambda}^{2}+\frac{4}{L_{2}}}}{2}\right){\mathbf{U}}^{\top},$
(34)

where $\boldsymbol{\Lambda}$ and ${\mathbf{U}}$ contain the eigenvalues and
eigenvectors of $\boldsymbol{\Theta}^{(t)}-\frac{1}{L_{2}}\left(\beta\nabla
f(\boldsymbol{\Theta}^{(t)})+\hat{\boldsymbol{\Sigma}}\right)$, respectively.

Then, ${\mathbf{Q}}$ is updated by addressing the following problem

$\begin{array}[]{rll}&\underset{{\mathbf{Q}}}{\mathsf{minimize}}&\frac{\beta}{2}\|{\mathbf{Q}}-{\mathbf{S}}\boldsymbol{\Theta}+\boldsymbol{\Theta}{\mathbf{S}}-\frac{1}{\beta}{\mathbf{Y}}\|_{F}^{2},\\\
&\mathsf{subject~{}to}&\|{\mathbf{Q}}\|_{F}\leq\delta,\end{array}$ (35)

where the iteration subscripts from $\boldsymbol{\Theta}^{(t+1)}$ and
${\mathbf{Y}}^{(t)}$ have been omitted for simplicity. Similar to the case of
updating ${\mathbf{P}}$, we update ${\mathbf{Q}}$ as

$\displaystyle{\mathbf{Q}}^{(t+1)}={\mathcal{P}}_{\delta}\left({\mathbf{S}}\boldsymbol{\Theta}-\boldsymbol{\Theta}{\mathbf{S}}+\frac{1}{\beta}{\mathbf{Y}}\right).$
(36)

Finally, the Lagrange multiplier ${\mathbf{Y}}$ is updated as follows

$\displaystyle{\mathbf{Y}}^{(t+1)}={\mathbf{Y}}^{(t)}+\beta({\mathbf{S}}\boldsymbol{\Theta}-\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{Q}}),$
(37)

where the iteration subscripts from $\boldsymbol{\Theta}^{(t+1)}$ and
${\mathbf{Q}}^{(t+1)}$ have been similarly omitted for clarity.

We can also use the adaptive strategy in (27) to adjust $\beta$ during
iterations, with $\mathbf{r}^{(t)}$ and $\mathbf{s}^{(t)}$ defined as

$\mathbf{r}^{(t)}={\mathbf{S}}\boldsymbol{\Theta}^{(t)}-\boldsymbol{\Theta}^{(t)}{\mathbf{S}}-\mathbf{Q}^{(t)},$

and

$\mathbf{s}^{(t)}=\beta^{(t)}{\mathbf{S}}\big{(}\mathbf{Q}^{(t)}-\mathbf{Q}^{(t-1)}\big{)}-\beta^{(t)}\big{(}\mathbf{Q}^{(t)}-\mathbf{Q}^{(t-1)}\big{)}{\mathbf{S}}.$

A pseudocode of the steps to be performed for the update of
$\boldsymbol{\Theta}$ is summarized in Algorithm 2.

###  IV-D Graph-learning algorithm and convergence analysis

Leveraging the results presented in the previous two subsections, the steps to
run our iterative scheme to find a solution
$(\hat{\boldsymbol{\Theta}},{\hat{\mathbf{S}}})$ to (8) are summarized in
Algorithm 3.

Before presenting the associated theoretical analysis, several comments
regarding the implementation of Algorithm 3 are in order:

  * •

For simplicity, the algorithm considers a fixed number of iterations, but a
prudent approach is to monitor the cost reduction at each iteration and
implement an early exit approach if no meaningful improvement is achieved.

  * •

The value of the hyperparameters $\rho$, $\eta$ and $\delta$ is an input to
the algorithm. We note that $\eta$ is typically set to a small value to
guarantee that the (sparsity promoting) $\ell_{1}$ norm plays a more prominent
role. Moreover, the value of $\delta$ should be chosen based on the quality of
the estimate of the sample covariance matrix $\hat{\boldsymbol{\Sigma}}$. The
higher the number of observations $R$ (hence, the better the quality of
$\hat{\boldsymbol{\Sigma}}$), the smaller the value of $\delta$. Similarly, if
the number of nodes $N$ is high, the value of $\delta$ should be re-scaled
accordingly, so that the constraint does not become too restrictive.

  * •

The update of ${\mathbf{S}}$ poses the primary computational challenge, mainly
due to the complex nature of its estimation. In contrast to
$\boldsymbol{\Theta}$, which is primarily estimated from the data matrix
$\hat{\boldsymbol{\Sigma}}$, the estimation of ${\mathbf{S}}$ relies on its
interplay with $\boldsymbol{\Theta}$ through the relaxed commutativity
constraint. This indirect relationship adds to the complexity of the
estimation, as it does not directly benefit from data-driven insights, often
necessitating greater precision in solving the corresponding subproblem.
Furthermore, the constraints imposed on ${\mathbf{S}}$ are substantially more
complex than those on $\boldsymbol{\Theta}$, thereby increasing the
computational load to obtain a solution that lies within the feasible set. To
alleviate the computational demand, we may employ a relatively loose stopping
criterion for the $\boldsymbol{\Theta}$ subproblem, which can expedite
convergence without significantly affecting the quality of the solution.

To establish the theoretical convergence of Algorithm 3, we begin by
introducing several definitions and (mild) assumptions pertinent to Problem
(8).

Let $f(\boldsymbol{\Theta},{\mathbf{S}})$ denote the objective function and
$\mathcal{X}$ the feasible set of Problem (8), respectively. We define
$(\bar{\boldsymbol{\Theta}},\bar{{\mathbf{S}}})$ as a block coordinatewise
minimizer of Problem (8) if:

$f(\bar{\boldsymbol{\Theta}},\bar{{\mathbf{S}}})\leq
f(\boldsymbol{\Theta},\bar{{\mathbf{S}}}),\quad\forall\,\boldsymbol{\Theta}\in\bar{\mathcal{X}}_{\boldsymbol{\Theta}},$
(38)

and

$f(\bar{\boldsymbol{\Theta}},\bar{{\mathbf{S}}})\leq
f(\bar{\boldsymbol{\Theta}},{\mathbf{S}}),\quad\forall\,{\mathbf{S}}\in\bar{\mathcal{X}}_{{\mathbf{S}}},$
(39)

where
$\bar{\mathcal{X}}_{\boldsymbol{\Theta}}=\\{\boldsymbol{\Theta}\in\mathbb{R}^{N\times
N}\,|\,(\boldsymbol{\Theta},\bar{{\mathbf{S}}})\in\mathcal{X}\\}$, and
$\bar{\mathcal{X}}_{{\mathbf{S}}}=\\{{\mathbf{S}}\in\mathbb{R}^{N\times
N}\,|\,(\bar{\boldsymbol{\Theta}},{\mathbf{S}})\in\mathcal{X}\\}$.

Input: $\hat{\boldsymbol{\Sigma}}$, $\rho$, $\eta$, $\beta$, and $\delta$

Outputs : $\hat{\boldsymbol{\Theta}}$ and ${\hat{\mathbf{S}}}$

1

2Initialize $\hat{\boldsymbol{\Theta}}^{(0)}=\hat{\boldsymbol{\Sigma}}^{-1}$

3 Initialize ${\hat{\mathbf{S}}}^{(0)}$ by solving (6)

4 Initialize ${\mathbf{P}}^{(0)}$, ${\mathbf{Q}}^{(0)}$, ${\mathbf{Y}}^{(0)}$,
and ${\mathbf{Z}}^{(0)}$ to zero

5for  _$k=0$ to $K-1$_ do

6

7 Update  ${\hat{\mathbf{S}}}^{(k+1)}$  by running Algorithm 1;

8

9 Update  $\hat{\boldsymbol{\Theta}}^{(k+1)}$  by running Algorithm 2;

10

11 end for

12 $\hat{\boldsymbol{\Theta}}=\hat{\boldsymbol{\Theta}}^{(K)}$,
${\hat{\mathbf{S}}}={\hat{\mathbf{S}}}^{(K)}$.

13

Algorithm 3 Polynomial Graphical Lasso (PGL)

Furthermore, we introduce the following assumptions for our analysis:

######  Assumption 1.

The parameter $\delta$ in Problem (8) is sufficiently large to ensure that the
feasible set of the subproblem (10) is nonempty at every iteration.

We require Assumption 1 because the intersection
$\mathcal{S}\cap\\{{\mathbf{S}}\in\mathbb{R}^{N\times
N}|,\|\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{S}}\boldsymbol{\Theta}\|_{F}\leq\delta\\}$
may otherwise be empty, implying that the feasible set of subproblem (10)
could be nonexistent. However, Assumption 1 is relatively mild, as we can
always choose a sufficiently large $\delta$ to ensure that the feasible set of
subproblem (10) remains nonempty at every iteration. Given Assumption 1, our
algorithm is guaranteed to find a minimizer of subproblem (10) throughout its
iterations. Additionally, the feasibility of subproblem (28) is inherently
assured.

######  Assumption 2.

The matrix $\hat{\boldsymbol{\Sigma}}$ in Problem (8) is positive definite.

Assumption 2, which requires all the eigenvalues to be nonzero, guarantees
that subproblem (28) is well defined. Without this assumption, the objective
function in (28) may fail to achieve a finite minimum value. In cases where
Assumption 2 may not hold, incorporating an norm regularizer for
$\boldsymbol{\Theta}$ would bound the solution, thereby ensuring the existence
of a minimizer.

(a)

(b)

(c)

Figure 1: Graph estimation error
$\mathrm{nme}({\mathbf{S}}^{*},{\hat{\mathbf{S}}})$ vs. number of samples $R$
for different graph learning approaches (PGL, GL, and GSR) and simulations
setups. Specifically, the six lines reported in each subplot correspond to the
combination of (a) 3 different graph learning methods and 2 different
covariance setups (SSEM and Poly) for a noise-free scenario; (b) 3 different
graph learning schemes and 2 noise levels $\sigma\in\\{0.05,0.2\\}$ for a Poly
setup; and (c) 3 graph generation models and 2 covariance models (SSEM and
Poly) for the PGL algorithm in a noise-free scenario.

######  Theorem 1.

Let
$\big{\\{}\big{(}\hat{\boldsymbol{\Theta}}^{(k)},{\hat{\mathbf{S}}}^{(k)}\big{)}\big{\\}}_{k\in\mathbb{N}}$
be a sequence generated by Algorithm 3. Under Assumptions 1 and 2, every limit
point of
$\big{\\{}\big{(}\hat{\boldsymbol{\Theta}}^{(k)},{\hat{\mathbf{S}}}^{(k)}\big{)}\big{\\}}_{k\in\mathbb{N}}$
is a block coordinatewise minimizer of Problem (8).

The proof of Theorem 1 is deferred to Appendix C. This theorem asserts the
subsequence convergence of our algorithm to a block coordinatewise minimizer
of Problem (8). When the block coordinatewise minimizer lies within the
interior of the feasible set, it becomes a stationary point. The theorem’s
assertions are significant both theoretically and practically. As discussed
earlier in the paper, GMRFs are a particular case of Gaussian graph-stationary
processes. Taking this into account, we can always initialize our algorithm
using the solution estimated by GL (which is optimal for GMRF) and then, run
iteratively the updates over $\boldsymbol{\Theta}$ and ${\mathbf{S}}$ in
Algorithm 3, to get an (enhanced) coordinatewise minimum estimate.

##  V Numerical experiments

This section evaluates quantitatively the performance of PGL. Since PGL can be
understood as a generalization of the widely adopted GL (with
$\boldsymbol{\Theta}$ being any polynomial in ${\mathbf{S}}$), in most test
cases, we will test both algorithms. Similarly, we also test the learning
performance of GSR [13], which assumes stationarity but not Gaussianity. The
performance results obtained from both synthetic and real-data experiments are
summarized in Figs. 1-4. Unless otherwise stated, to assess the quality of the
estimated GSO, we use the normalized mean error between the estimated and true
${\mathbf{S}}$. Mathematically, this entails computing111Results for other
recovery metrics (including accuracy and F1 score) as well as additional
simulations can be found both in our conference precursor [1] and in our
online repository https://github.com/andreibuciulea/GaussSt˙TopoID.

$\mathrm{nme}({\mathbf{S}}^{*},{\hat{\mathbf{S}}})=\frac{\|{\mathbf{S}}^{*}-{\hat{\mathbf{S}}}\|_{F}^{2}}{\|{\mathbf{S}}^{*}\|_{F}^{2}},$
(40)

where ${\hat{\mathbf{S}}}$ and ${\mathbf{S}}^{*}$ represent the estimated and
true ${\mathbf{S}}$ respectively. Moreover, for the synthetic experiments we
test the graph learning algorithms over $100$ realizations of random graphs
$\\{{\mathcal{G}}_{i}\\}_{i=1}^{100}$ and report the average normalized mean
error
$\frac{1}{100}\sum_{i=1}^{100}\mathrm{nme}({\mathbf{S}}_{i}^{*},{\hat{\mathbf{S}}}_{i})$.

If not specified otherwise, in our synthetic experiments, we consider graphs
with $N=20$ nodes generated using the Erdös-Rényi (ER) model with a link
probability of $p=0.1$. Regarding the generation of the graph signals, three
different setups for the covariance matrices have been studied. For the setup
referred to as “Poly”, the covariance matrix $\boldsymbol{\Sigma}$ is
generated as a random polynomial of the GSO of the form
$\boldsymbol{\Sigma}=(\sum_{l=0}^{L-1}h_{l}{\mathbf{S}}^{l})^{2}$, where
$h_{l}$ are random coefficients drawn from a normalized zero-mean Gaussian
distribution, and the square operator guarantees matrix $\boldsymbol{\Sigma}$
to be positive definite. The setup referred to as “SSEM” constructs the
covariance matrices following the sparse structural equation model [34] for
graph signal generation as
$\boldsymbol{\Sigma}=({\mathbf{I}}-{\mathbf{S}})^{2}$, where ${\mathbf{S}}$ is
selected to ensure that $\boldsymbol{\Sigma}$ is positive definite. The setup
referred to as “MRF” constructs the covariance matrices following the
assumptions made by GL as
$\boldsymbol{\Sigma}=(\mu{\mathbf{I}}+\nu{\mathbf{S}})^{-1}$, where $\mu$ is
some positive number large enough to assure that $\boldsymbol{\Sigma}^{-1}$ is
positive definite and $\nu$ is some positive random number.

###  V-A Test case 1: Estimation error vs. number of samples for multiple
synthetic scenarios.

In this first test case, we employ synthetic scenarios for testing the
performance of our approach in terms of
nme$({\mathbf{S}}^{*},{\hat{\mathbf{S}}})$ vs. $R$ and also compare the
results with other methods from the literature. The different scenarios
considered are detailed below.

Error vs. number of samples for different graph learning methods and signal
models. The results of the experiment depicted in Fig. 1a, compare the
nme$({\mathbf{S}}^{*},{\hat{\mathbf{S}}})$ (y-axis) of various algorithms with
respect to the number of available samples $R$ (x-axis). Moreover, we consider
SSEM and Poly setups for data generation. The results shown in Fig. 1a reveal
that: i) PGL outperforms its competitors, ii) the error decreases as $R$
increases, and iii) estimation is more accurate for SSEM than for Poly. Next,
we discuss these three main findings in greater detail. All algorithms do a
better job estimating the graph for the SSEM model. Since
$\boldsymbol{\Theta}$ for SSEM is a specific second-order polynomial in
${\mathbf{S}}$ it can be seen as a particular case of Poly, and consequently,
a scenario from which the graph structure is easier to estimate. Indeed, while
GL fails to estimate the graph for the Poly model, it is able to estimate some
of the links for the SSEM. However, the estimation error of GL is quite large
and does not change with the number of samples $R$, demonstrating that the
poor performance is due to a model mismatch. This will be further confirmed in
Section V-B, where we simulate a GMRF data generation setup that aligns
perfectly with the assumptions made by GL. Regarding PGL and GSR, we observe
that the error decreases almost linearly with the number of samples $R$.
Perhaps more importantly, we also observe that, as $R$ increases, the gap
between PGL and GSR diminishes. For example, while for the Poly case GSR needs
10 times more samples than PGL to achieve an error of $10^{-1}$, GSR only
needs 3 times more samples than PGL to achieve an error of $10^{-3}$. This
illustrates that, as expected, the gains associated with assuming Gaussianity
are stronger when the number of observations $R$ is small, vanishing as $R$
grows very large.

Error vs. number of samples for noisy observations. Next, we assess the
performance of the graph-learning algorithms in the presence of additive white
Gaussian observation noise. The results are shown in Fig. 1b. As in the
previous test case, we report nme$({\mathbf{S}}^{*},{\hat{\mathbf{S}}})$ vs
$R$ for PGL, GL and GSR. The difference here is that we consider only the more
intricate signal generation model (Poly) and two normalized noise levels
($\sigma$ = 0.05, $\sigma$ = 0.2). The main observations in this case are: i)
PGL outperforms GSR and GL, ii) the error for PGL and GSR decreases as $R$
increases, while the one for GL is flat and close to $1$, iii) the estimation
performance for PGL and GSR worsens as the noise level $\sigma$ increases,
deteriorating noticeably with respect to the setup in Fig. 1a, and iv) the gap
between PGL and GSR grows. The findings in i) and ii) are consistent with
those found in the previous test case. Finding iii) is expected and common in
all graph-learning approaches. Finally, the larger gap in finding iv) is due
to the fact that this is a more challenging scenario (high-order polynomial
covariances and observation noise), and, as a result, the higher level of
sophistication of PGL relative to GSR translates into more noticeable gains.

Error vs. number of samples for different graph models. This test case
considers network models other than ER. In particular, three different types
of graphs are considered: 1) Small World (SW) graphs with mean node degree $4$
and rewiring probability $0.15$; 2) Stochastic block model (SBM) graphs with
$4$ clusters, and intra and inter-cluster edge probability of 0.8 and 0.05,
respectively; and 3) Barabási-Albert (BA) graphs with $2$ edges to attach at
every step. As in the first test case, we consider two types of signal
generation models: SSEM and Poly. Fig. 1c reports the error vs. the number of
samples for the six combinations considered (3 types of graphs and 2 types of
signal generation models). The results show that there is a significant
difference in performance between SW, SBM, and BA, which is part due to the
sparsity level present in each graph. One of the assumptions codified in our
model is that the graph is sparse, and, as a result, our algorithm does a
better job estimating BA (the one with the lowest average degree) than SBM and
SW (the one with the highest average degree). Finally, we also note that the
estimation error achieved with SSEM consistently outperforms that of the Poly
setup. These findings align with the results presented in Fig. 1a for ER
graphs and are in accordance with the theoretical discussion that postulated
SSEM as a specific instance of Poly.

###  V-B Test case 2: Noisy GMRF graph signals.

In this test case, the goal is to assess the behavior of PGL for GMRF
observations, which is the setup that motivated the development of the GL
algorithm.

Figure 2: Graph estimation error $\mathrm{nme}({\mathbf{S}}^{*},{\hat{\mathbf{S}}})$ vs. number of samples $R$ for data generated according to a GMRF. We consider three different graph learning schemes (PGL, GL, and GSR) and two levels of additive white noise ($sigma\in\\{0.05,0.2\\}$), giving rise to the six lines in the figure. Alg. $\backslash$ $N$ | 20 | 30 | 40 | 50 | 60 | 70 | 80 | Metric
---|---|---|---|---|---|---|---|---
PGL-CVX | $2.37\cdot 10^{1}$ | $3.88\cdot 10^{1}$ | $1.29\cdot 10^{2}$ | $1.11\cdot 10^{3}$ | $3.13\cdot 10^{3}$ | —- | —- | Time (s)
PGL-Alg.3 | $2.49\cdot 10^{0}$ | $2.72\cdot 10^{0}$ | $2.88\cdot 10^{0}$ | $3.81\cdot 10^{0}$ | $4.45\cdot 10^{0}$ | $5.86\cdot 10^{0}$ | $6.75\cdot 10^{0}$
PGL-CVX | $7.98\cdot 10^{-4}$ | $3.12\cdot 10^{-3}$ | $1.35\cdot 10^{-2}$ | $2.59\cdot 10^{-2}$ | $9.16\cdot 10^{-2}$ | —- | —- | nme$({\mathbf{S}}^{*},{\hat{\mathbf{S}}})$
PGL-Alg.3 | $4.66\cdot 10^{-4}$ | $1.30\cdot 10^{-3}$ | $2.18\cdot 10^{-3}$ | $5.07\cdot 10^{-3}$ | $1.59\cdot 10^{-2}$ | $1.43\cdot 10^{-1}$ | $7.94\cdot 10^{-1}$
TABLE I: Values of the running time and
$\mathrm{nme}({\mathbf{S}}^{*},{\hat{\mathbf{S}}})$ for two different
implementations of PGL:one using an off-the-shelf convex solver (PGL-CVX) and
the other one using the efficient method in Algorithm 3 (PGL-Alg.3). The
results are shown for different graph sizes $N$ and $-$ indicates that the
running time was more than 2 hours.

Estimation error considering noisy GMFR signals. In this experiment, we
replicate the scenario from Fig.1b, utilizing a GMFR model to generate the
signals. The performance of PGL, GL and GSR for two different noise levels,
$\sigma\in\\{0.05,0.2\\}$, is depicted in Fig. 2. The main observations are:
i) across all considered approaches, increasing the noise level $\sigma$ leads
to a deterioration in terms of nme$({\mathbf{S}}^{*},{\hat{\mathbf{S}}})$; ii)
GSR always performs worse than PGL, providing very poor results when
$R<10^{4}$; and iii) GL outperforms PGL when the number of observations $R$ is
small. Findings i) and ii) are unsurprising and consistent with the behavior
observed in the previous experiments, showcasing the benefits of considering
the log-likelihood regularization in the optimization run by PGL. Regarding
iii), GL outperforming PGL is expected, since the latter needs to ”learn” the
particular polynomial between $\boldsymbol{\Theta}$ and ${\mathbf{S}}$.

However, as the value of $R$ increases, the error in PGL gradually decreases,
while the error in GL remains relatively constant, leading to PGL
outperforming GL for large values of $R$. This behavior is more surprising and
can be attributed to the fact that GL focuses on learning the precision matrix
$\boldsymbol{\Theta}$, while PGL balances the accuracy in terms of both the
precision $\boldsymbol{\Theta}$ and the graph ${\mathbf{S}}$. Since the
reported error focuses on the estimation ${\mathbf{S}}$, the values of the
diagonal of $\boldsymbol{\Theta}$ are not relevant for
nme$({\mathbf{S}}^{*},{\hat{\mathbf{S}}})$ and this can explain that GL, which
is a maximum likelihood estimate for $\boldsymbol{\Theta}$, does not yield the
minimum nme$({\mathbf{S}}^{*},{\hat{\mathbf{S}}})$.

###  V-C Test case 3: Computational complexity.

Here, we compare the runtime obtained by the efficient implementation of our
PGL scheme provided in Algorithm 3 and a generic block-coordinate alternating
minimization algorithm that uses CVX [35], the most common off-the-shelf
solver for convex problems.

Computational complexity and estimation error. The objective of this
experiment is to evaluate the running time and estimation error for different
versions of our algorithm as the number of nodes increases. In particular, the
experiment focuses on the Poly setup, utilizing $R=10^{6}$ graph signals, and
averages the results over 50 graph realizations. Table I lists the elapsed
time required to obtain the graph and precision estimates for problems with
different numbers of nodes $N$ using two algorithms: 1) solving the
optimization in (8) with a block coordinate approach where the minimization
over $\boldsymbol{\Theta}$ given ${\mathbf{S}}$ and the minimization over
${\mathbf{S}}$ given $\boldsymbol{\Theta}$ are run using CVX (this algorithm
is labeled as PGL-CVX) and 2) employing the efficient scheme outlined in
Algorithm 3 (this choice is labeled as PGL-Alg.3). To guarantee that the
results are comparable, the nme$({\mathbf{S}}^{*},{\hat{\mathbf{S}}})$ is also
reported. Examining the listed running times, we observe that as the number of
nodes increases, both solvers require more time to estimate the graph (note
that the number of variables scales with $N^{2}$). More importantly, there
exists a noticeable difference between PGL-CVX and PGL-Alg.3. Not only the
latter is faster for small graphs, but the gains grow significantly as $N$
increases. Note that the results for graphs with more than $N=60$ nodes are
not reported for PGL-CVX, since the computation time exceeded two hours. In
terms of nme$({\mathbf{S}}^{*},{\hat{\mathbf{S}}})$, our approach achieves
similar (slightly better) results than PGL-CVX. In conclusion, the experiment
demonstrates that the efficient implementation described in Section IV is more
efficient than readily available convex solvers, rendering it particularly
well-suited for large graphs while yielding comparable
nme$({\mathbf{S}}^{*},{\hat{\mathbf{S}}})$.

###  V-D Test case 4: Real data scenarios.

Finally, we compare different graph-learning algorithms (including PGL) in the
context of two different graph-aware applications. The details and results are
provided next.

Stock graph-based clustering from returns. For this real-data experiment, we
tested our graph learning algorithm on financial data and further performed a
clustering task. Specifically, 40 companies from 4 different sectors of the
S&P 500 were selected (10 companies from each sector) and the market data (log
returns) of each company in the period 2010-2015 was retrieved. This gives
rise to a data matrix of size ${\mathbf{X}}\in{\mathbb{R}}^{40\times 1510}$.
In this application, as in many others dealing with graph learning, we do not
have access to the ground truth graph. Hence, we cannot quantify the quality
of the estimated network directly by using a metric like $\mathrm{nme}$ or F1
score. As a result, we need to assess the quality of the estimated graph
indirectly, using the graph as input for an ulterior task. In this experiment,
we use the graph to estimate the community each company belongs to in an
unsupervised manner. More specifically, we implement the following pipeline:
1) estimate several graphs from a subset of the available graph signals, 2)
use spectral clustering to group the nodes into 4 communities (as many as
sectors were selected), and 3) compute the ratio of incorrectly clustered
nodes. To obtain more reliable results, we averaged the clustering errors over
50 realizations in which the subset of graph signals was chosen uniformly at
random. The vertical axis of Fig.3 represents the normalized clustering error
$\frac{1}{50}\sum_{i=1}^{50}\frac{N_{w}^{(i)}}{N}$ which is computed as the
average of the fraction between the number of wrongly clustered nodes $N_{w}$
and the total number of nodes $N$ over $50$ graph realizations. The horizontal
axis of Fig. 3 represents the percentages of graph signals used to estimate
${\hat{\mathbf{S}}}$. Results are provided for 4 graph-learning approaches:
PGL, GL, GSR, and “Corr”, which estimates the graph as a thresholded version
of the correlation matrix. The idea is that companies from the same sector
have stronger ties among them and, as a result, when running a graph-based
clustering method, the 4 sectors should arise. Based on the results presented
in Fig.3, it can be observed that PGL outperforms the other alternatives and
additionally, as the number of available signals increases, the clustering
error for PGL drops significantly. The superiority of PGL may indicate that
considering more complex relationships among stocks (beyond the simple
correlations considered in “Corr” or the partial correlations considered in
GL) is a better model to understand the dependencies between log-returns in
the stock market. On the other hand, GSR shows high clustering error with a
limited number of samples, but it improves as the number of available samples
increases. This observation aligns with our earlier discussion in Section III,
where we highlighted that this particular model offers greater generality but
necessitates a substantial number of samples to accurately estimate the graph.

Figure 3: The experiment considers 40 companies of the S&P 500 from 4
different sectors, using as nodal signals the daily returns in the period
2010-2015. Four different graph-learning schemes are considered: correlation
networks (Corr), PGL, GL and GSR. After learning the graph, a spectral
clustering method is implemented. The figure shows the normalized node
clustering error (fraction of wrongly clustering nodes) as the percentage of
available signals to learn the graph increases.

Learning sequential graphs for investing. This experiment deals with a
different real-world problem and dataset. We still look at stocks, but
consider now the close price of the 7 FAAMUNG222Facebook, Amazon, Apple,
Microsoft, Uber, Netflix, and Google companies from Jul 2019 to May 2020. The
goal is to design an investment strategy to maximize the benefits using as
input a graph describing the relationships among the companies. Inspired by
the approach in [36], we first build a graph, analyze its connectivity and
then, invest (or not) in a stock according to the graph connectivity. To be
more specific, we use the close price to estimate multiple $7\times 7$
adjacency matrices. We estimate a total of 200 matrices, where each adjacency
matrix (graph) is estimated in a rolling window fashion. The window consists
of 30 consecutive days and for each graph estimation, we shift the window one
day at a time. These graph estimations help to visualize how the graph
topology changes during this time period. The finding in [36] is that big
changes in the graph connectivity indicate opportunities to invest. To that
end, we keep track of the algebraic connectivity value, which is the second
smallest eigenvalue ($\lambda_{2}$) of the estimated Laplacian matrix. The
lower the value, the less connected the graph is and, as a result, the easier
breaking the graph into multiple components is. Fig.4 shows the value of
$\lambda_{2}$ for each of the 200 considered graphs (each associated with a
30-day period). Then, using the approach in [36] we invest only if
$\lambda_{2}$ is below a fixed threshold. We learn the graph using 3 methods
(GL, GSR, and PGL) and, for each of them, we select the best possible
threshold (the one that maximizes the benefits). Correlation was not used here
due to its poor performance. The results of applying the graph-connectivity-
based investment strategy to the graphs estimated with each of the algorithms
are shown in Fig. 4. The blue line labeled as “Strategy I” represents the
benefits of investing every day the available amount and is used as a
baseline. By analyzing the results obtained, we can observe that: i) the
graph-based strategies outperform (gain more money than) the baseline; ii) the
strategies based on GL and GSR provide similar gains; and iii) the strategy
based on GSSR yields the highest gains. This provides additional validation
for the graph-learning methodology proposed in this paper.

(a)

(b)

Figure 4: This experiment learns the graph connecting the 7 FAAMUNG stocks
using as signals the closing price from July 2019 to May 2020. Three different
graph learning methods are used (PGL, GL, and GSR) and a different graph is
learned for every day (using the signals of the previous 30 days). Subplot (a)
shows the value of the algebraic connectivity ($\lambda_{2}$) associated with
each one of the $200\times 3$ estimated graphs. Subplot (b) shows the value of
the portfolio for 4 different investment strategies, 3 of which are based on
the algebraic connectivity estimated in the subplot (a).

##  VI Conclusions

This paper has introduced PGL, a novel scheme for learning a graph from nodal
signals, with our key contribution being the modeling of the signals as
Gaussian and stationary on the graph. This approach opens the door to a graph-
learning formulation that leverages the advantages of GL (needing a relatively
small number of signals to get a good estimation of the graph structure) while
encompassing a more comprehensive model (because it handles cases where the
precision matrix can be any polynomial form of the sought graph). Given the
increased complexity and nonconvex nature of the resulting optimization
problem, we have developed a low-complexity algorithm that alternates between
estimating the graph and precision matrices and have characterized its
convergence to a block coordinatewise minimum. To assess its efficacy, we have
conducted numerical simulations comparing PGL with various alternatives, using
both synthetic and real data. The results have showcased the benefits of our
approach, motivating the adoption and further investigation of the proposed
graph-learning methodology.

##  Appendix A Computations of projections

We present the details about how to compute the projections
$\mathcal{P}_{\mathcal{S}_{A}}$ and $\mathcal{P}_{\mathcal{S}_{B}}$.

The computation of $\mathcal{P}_{\mathcal{S}_{A}}$ is straightforward as
follows:

$\mathcal{P}_{\mathcal{S}_{A}}(\mathbf{A})=\big{(}\mathbf{A}+\mathbf{A}^{\top}\big{)}/2.$
(41)

The projection $\mathcal{P}_{\mathcal{S}_{B}}$ is defined as the minimizer of
the optimization problem as follows:

$\mathcal{P}_{\mathcal{S}_{B}}(\mathbf{A}):=\underset{\mathbf{X}\in\mathcal{S}_{B}}{\mathsf{arg\,min}}\,\frac{1}{2}\|\mathbf{X}-\mathbf{A}\|_{F}^{2}.$
(42)

To compute the projection $\mathcal{P}_{\mathcal{S}_{B}}$, we solve Problem
(42) row by row. For the $j$-th row, we solve the following problem,

$\begin{array}[]{rll}&\underset{\mathbf{x}\in\mathbb{R}^{N}}{\mathsf{minimize}}&\frac{1}{2}\|\mathbf{x}-\mathbf{a}\|^{2},\\\
&\mathsf{subject~{}to}&\mathbf{x}^{\top}\mathbf{1}\geq
1,\,x_{j}=0,\,\mathbf{x}_{\setminus j}\geq\mathbf{0},\end{array}$ (43)

where $\mathbf{a}\in\mathbb{R}^{N}$ contains all entries of the $j$-th row of
$\mathbf{A}$, $x_{j}$ denotes the $j$-th entry of $\mathbf{x}$, and
$\mathbf{x}_{\setminus j}\in\mathbb{R}^{N-1}$ contains all entries of
$\mathbf{x}$ except the $j$-th one.

Let $\hat{\mathbf{x}}$ denote the optimal solution of Problem (43).
Proposition 1 below, proved in Appendix B, presents the optimal solution of
Problem (43).

######  Proposition 1.

The optimal solution $\hat{\mathbf{x}}$ to Problem (43) can be obtained as
follows:

  * •

If $\sum_{i\neq j}\max(a_{i},\,0)\geq 1$, then $\hat{x}_{j}=0$, and
$\hat{x}_{i}=\max(a_{i},\,0)$, for $i\neq j$.

  * •

If $\sum_{i\neq j}\max(a_{i},\,0)<1$, then $\hat{x}_{j}=0$, and
$\hat{x}_{i}=\max(a_{i}+\phi,\,0)$, for $i\neq j$, where $\phi$ satisfies
$\sum_{i\neq j}\max(a_{i}+\phi,\,0)=1$.

Several efficient approaches have been developed to tackle the piecewise
linear equation $\sum_{i\neq j}\max(a_{i}+\phi,0)=1$. Among these, the
sorting-based method described in [37] is noteworthy. Central to this method
is the sorting of the vector $\mathbf{a}$, which constitutes the most
computationally intensive step, generally necessitating $\mathcal{O}(N\log N)$
operations.

##  Appendix B Proof of Proposition 1

###### Proof.

The Lagrangian of the optimization in (43) is

$L(\mathbf{x},u,\mathbf{v})=\frac{1}{2}\|\mathbf{x}-\mathbf{a}\|^{2}-u\big{(}\mathbf{x}^{\top}\mathbf{1}-1\big{)}-\langle\mathbf{v}_{\setminus
j},\,\mathbf{x}_{\setminus j}\rangle+v_{j}x_{j},$

where $u\in\mathbb{R}$ and $\mathbf{v}\in\mathbb{R}^{N}$ are KKT multipliers.
Let $(\hat{\mathbf{x}},\hat{u},\hat{\mathbf{v}})$ be the primal and dual
optimal point. Then $(\hat{\mathbf{x}},\hat{u},\hat{\mathbf{v}})$ must satisfy
the KKT system:

$\displaystyle\hat{x}_{i}-a_{i}-\hat{u}-\hat{v}_{i}$
$\displaystyle=0,\quad\mathrm{for}\ i\neq j;$ (44)
$\displaystyle\hat{x}_{j}-a_{j}-\hat{u}+\hat{v}_{j}$ $\displaystyle=0;$ (45)
$\displaystyle\ \hat{x}_{i}\geq 0,\ \hat{v}_{i}\geq 0,\
\hat{v}_{i}\hat{x}_{i}$ $\displaystyle=0,\quad\mathrm{for}\ i\neq j;$ (46)
$\displaystyle\hat{x}_{j}=0,\ \hat{u}\geq 0,\ \hat{\mathbf{x}}^{T}\mathbf{1}$
$\displaystyle\geq 1;$ (47)
$\displaystyle\hat{u}\big{(}\hat{\mathbf{x}}^{T}\mathbf{1}-1\big{)}$
$\displaystyle=0;$ (48)

Therefore, for any $i\neq j$, it holds that
$\hat{x}_{i}=a_{i}+\hat{u}+\hat{v}_{i}$. Then we obtain the following results:

  * •

If $a_{i}+\hat{u}<0$, then $\hat{v}_{i}=-a_{i}+\hat{u}$ and $\hat{x}_{i}=0$,
following from the fact that $\hat{v}_{i}\hat{x}_{i}=0$ and $\hat{x}_{i}\leq
0$.

  * •

If $a_{i}+\hat{u}_{r}\geq 0$, then $\hat{v}_{i}=0$. This can be obtained as
follows: if $\hat{x}_{i}=0$, then $\hat{v}_{i}=-(a_{i}+\hat{u})\leq 0$. Since
$\hat{v}_{i}\geq 0$, one has $\hat{v}_{i}=0$; On the other hand, if
$\hat{x}_{i}\neq 0$, then $\hat{v}_{i}=0$, following from the fact that
$\hat{v}_{i}\hat{x}_{i}=0$. As a result, we get $\hat{x}_{i}=a_{i}+\hat{u}$.

Overall, we obtain that

$\hat{x}_{j}=0,\quad\mathrm{and}\quad\hat{x}_{i}=\max(a_{i}+\hat{u},\,0),\quad\forall\,i\neq
j.$ (49)

On the other hand, $\hat{\mathbf{x}}$ and $\hat{u}$ satisfy that
$\hat{\mathbf{x}}^{T}\mathbf{1}\geq 1$, $\hat{u}\geq 0$, and
$\hat{u}\big{(}\hat{\mathbf{x}}^{T}\mathbf{1}-1\big{)}=0$. To this end, we can
obtain the following results:

  * •

If $\sum_{i\neq j}\max(a_{i},\,0)\geq 1$, then $\hat{u}=0$, indicating that
$\hat{x}_{i}=\max(a_{i},\,0)$, for any $i\neq j$.

  * •

If $\sum_{i\neq j}\max(a_{i},\,0)<1$, then $\hat{u}\neq 0$. This is because
$\hat{u}=0$ will result in $\hat{\mathbf{x}}^{\top}\mathbf{1}<1$, which does
not satisfy the KKT system. Together with the KKT condition that
$\hat{u}\big{(}\hat{\mathbf{x}}^{T}\mathbf{1}-1\big{)}=0$, one has
$\hat{\mathbf{x}}^{\top}\mathbf{1}=1$. Therefore, one obtains that, for any
$i\neq j$, $\hat{x}_{i}=\max(a_{i}+\hat{u},\,0)$, where $\hat{u}$ is chosen
such that $\sum_{i\neq j}\max(a_{i}+\hat{u},\,0)=1$.

We note that the $\phi$ in Proposition 1 is exactly the dual optimal point
$\hat{u}$. ∎

##  Appendix C Proof of Theorem 1

###### Proof.

The convergence result stated in Theorem 1 is based on the framework
established by Theorem 2.3 in [38]. To demonstrate the validity of Theorem 1,
it suffices to establish that the conditions and assumptions of Theorem 2.3
are satisfied in our context. Our approach to block updates aligns with the
procedure delineated in equation (1.3a) of [38].

We first verify the conditions the requisite conditions of Theorem 2.3 in [38]
are met. Specifically, Theorem 2.3 stipulates that the objective function,
along with the feasible set of the optimization problem, should exhibit block
multiconvexity. For Problem (8), the objective function $f$ is convex with
respect to each of the blocks $\boldsymbol{\Theta}$ and ${\mathbf{S}}$ when
the other block is fixed, a property that defines block multiconvexity as per
[38]. Moreover, the function $f$ is strongly convex with respect to both
$\boldsymbol{\Theta}$ and ${\mathbf{S}}$.

The feasibility constraints of Problem (8) form a set $\mathcal{X}$ that
satisfies the criteria for block multiconvexity as defined in [38]. This is
due to the convexity of the individual set maps
$\mathcal{X}_{\boldsymbol{\Theta}}$ and $\mathcal{X}_{{\mathbf{S}}}$. The set
map $\mathcal{X}_{\boldsymbol{\Theta}}$ is defined as

$\mathcal{X}_{\boldsymbol{\Theta}}=\\{\boldsymbol{\Theta}\in\mathbb{R}^{N\times
N}\,|\boldsymbol{\Theta}\succeq\mathbf{0},\|\,\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{S}}\boldsymbol{\Theta}\|_{F}\leq\delta\\}$
(50)

for some given ${\mathbf{S}}$, and similarly, the set map
$\mathcal{X}_{{\mathbf{S}}}$ is defined as

$\mathcal{X}_{{\mathbf{S}}}=\\{{\mathbf{S}}\in\mathbb{R}^{N\times
N}\,|\,{\mathbf{S}}\in\mathcal{S},\|\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{S}}\boldsymbol{\Theta}\|_{F}\leq\delta\\}$
(51)

for some given $\boldsymbol{\Theta}$. Consequently, the optimization
subproblems with respect to $\boldsymbol{\Theta}$ and ${\mathbf{S}}$ in
Problem (8) are convex.

We now validate the assumptions required by Theorem 2.3 in [38] within the
context of our setting. Specifically, it is required that the objective
function $f$ is bounded below over the feasible set $\mathcal{X}$, that is,
$\inf_{(\boldsymbol{\Theta},{\mathbf{S}})\in\mathcal{X}}f(\boldsymbol{\Theta},{\mathbf{S}})>-\infty$.
This is indeed the case here, because the term
$\rho\|{\mathbf{S}}\|_{1}+\frac{\eta}{2}\|{\mathbf{S}}\|_{F}^{2}$ is
nonnegative. Additionally, the function
$-\log(\det(\boldsymbol{\Theta}))+\text{tr}(\hat{\boldsymbol{\Sigma}}\boldsymbol{\Theta})$
attains a finite infimum when $\hat{\boldsymbol{\Sigma}}$ is positive definite
under Assumption 2. To see this, first observe that
$\det(\boldsymbol{\Theta})\leq\|\boldsymbol{\Theta}\|_{2}^{N}$, where
$\|\boldsymbol{\Theta}\|_{2}$ is the largest eigenvalue of
$\boldsymbol{\Theta}$. Consequently,
$-\log\det(\boldsymbol{\Theta})\geq-N\log(\|\boldsymbol{\Theta}\|_{2})$.
Further, since
$\text{tr}(\hat{\boldsymbol{\Sigma}}\boldsymbol{\Theta})\geq\gamma\text{tr}(\boldsymbol{\Theta})\geq\gamma\|\boldsymbol{\Theta}\|_{2}$,
with $\gamma>0$ being the smallest eigenvalue of $\hat{\boldsymbol{\Sigma}}$,
we obtain

$-\log\det(\boldsymbol{\Theta})+\text{tr}(\hat{\boldsymbol{\Sigma}}\boldsymbol{\Theta})\geq-N\log(\|\boldsymbol{\Theta}\|_{2})+\gamma\|\boldsymbol{\Theta}\|_{2},$
(52)

which indeed has a finite minimum. Thus, we conclude that
$\inf_{(\boldsymbol{\Theta},{\mathbf{S}})\in\mathcal{X}}-\log\det(\boldsymbol{\Theta})+\text{tr}(\hat{\boldsymbol{\Sigma}}\boldsymbol{\Theta})>-\infty$.

Furthermore, the existence of block coordinatewise minimizers is assured by
the compactness of the feasible set. Moreover, Theorem 2.3 in [38] stipulates
that set maps change continuously during iterations. Referring to (50) and
(51), it is clear that the only constraint that changes through iterations is
$\|\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{S}}\boldsymbol{\Theta}\|_{F}\leq\delta$,
while the other constraints, $\boldsymbol{\Theta}\succeq\mathbf{0}$ and
${\mathbf{S}}\in\mathcal{S}$, remain constant. Given that
$\|\boldsymbol{\Theta}{\mathbf{S}}-{\mathbf{S}}\boldsymbol{\Theta}\|_{F}$ is a
continuous function with respect to both $\boldsymbol{\Theta}$ and
${\mathbf{S}}$, the set maps indeed change continuously, satisfying the
theorem’s conditions.

These verifications above have demonstrated that the conditions and
assumptions of Theorem 2.3 are satisfied in our context, completing the proof.
∎

## References

  * [1] A. Buciulea and A. G. Marques,  “Graph learning from Gaussian and stationary graph signals,”  in IEEE Intl. Conf. Acoust., Speech and Signal Process. (ICASSP), 2023, pp. 1–5.
  * [2] E. D. Kolaczyk,  Statistical Analysis of Network Data: Methods and Models,  Springer, New York, NY, 2009.
  * [3] O. Sporns,  Discovering the Human Connectome,  MIT Press, Boston, MA, 2012.
  * [4] D.I Shuman, S.K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst,  “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,”  IEEE Signal Process. Mag., vol. 30, no. 3, pp. 83–98, May 2013\.
  * [5] A. Ortega, P. Frossard, J. Kovacevic, J. M. F. Moura, and P. Vandergheynst,  “Graph signal processing: Overview, challenges, and applications,”  Proc. IEEE, vol. 106, no. 5, pp. 808–828, 2018.
  * [6] J. V. D. M. Cardoso, J. Ying, and D. P Palomar,  “Algorithms for learning graphs in financial markets,”  arXiv preprint arXiv:2012.15410, 2020.
  * [7] G. Mateos, S. Segarra, A. G. Marques, and A. Ribeiro,  “Connecting the dots: Identifying network structure via graph signal processing,”  IEEE Signal Process. Mag., vol. 36, no. 3, pp. 16–43, 2019.
  * [8] S. Sardellitti, S. Barbarossa, and P. Di Lorenzo,  “Graph topology inference based on sparsifying transform learning,”  IEEE Trans. Signal Process., vol. 67, no. 7, pp. 1712–1727, 2019\.
  * [9] N. Meinshausen and P. Buhlmann,  “High-dimensional graphs and variable selection with the lasso,”  Ann. Statist., vol. 34, pp. 1436–1462, 2006.
  * [10] G. V. Karanikolas, G. B. Giannakis, K. Slavakis, and R. M. Leahy,  “Multi-kernel based nonlinear models for connectivity identification of brain networks,”  in IEEE Intl. Conf. Acoust., Speech and Signal Process. (ICASSP), Shanghai, China, Mar. 20-25, 2016.
  * [11] J. Mei and J.M.F. Moura,  “Signal processing on graphs: Estimating the structure of a graph,”  in IEEE Intl. Conf. Acoust., Speech and Signal Process. (ICASSP), 2015, pp. 5495–5499.
  * [12] H. E. Egilmez, E. Pavez, and A. Ortega,  “Graph learning from data under Laplacian and structural constraints,”  IEEE J. Sel. Topics Signal Process., vol. 11, no. 6, pp. 825–841, 2017.
  * [13] S. Segarra, A. G. Marques, G. Mateos, and A. Ribeiro,  “Network topology inference from spectral templates,”  IEEE Trans. Signal Info. Process. Networks, vol. 3, no. 3, pp. 467–483, Sep. 2017.
  * [14] V. Chandrasekaran, P. A. Parrilo, and A. S. Willsky,  “Latent variable graphical model selection via convex optimization,”  Ann. of Statist., vol. 40, no. 4, pp. 1935–1967, 2012.
  * [15] A. Chang, T. Yao, and G. I. Allen,  “Graphical models and dynamic latent factors for modeling functional brain connectivity,”  in 2019 IEEE Data Science Wrksp. (DSW). IEEE, 2019, pp. 57–63.
  * [16] A. Buciulea, S. Rey, and A. G. Marques,  “Learning graphs from smooth and graph-stationary signals with hidden variables,”  IEEE Trans. Signal Process., vol. 8, pp. 273–287, 2022.
  * [17] P. Danaher, P. Wang, and D. M. Witten,  “The joint graphical lasso for inverse covariance estimation across multiple classes,”  J. of the Roy. Statistical Soc.: Ser. B (Statistical Methodology), vol. 76, no. 2, pp. 373–397, 2014.
  * [18] M. Navarro, Y. Wang, A. G. Marques, C. Uhler, and S. Segarra,  “Joint inference of multiple graphs from matrix polynomials,”  J. of Machine Learning Research (JMLR), vol. 23, no. 1, pp. 3302–3336, 2022.
  * [19] A. G. Marques, S. Segarra, G. Leus, and A. Ribeiro,  “Stationary graph processes and spectral estimation,”  IEEE Trans. Signal Process., vol. 65, no. 22, pp. 5911–5926, 2017\.
  * [20] A. Sandryhaila and J.M.F. Moura,  “Discrete signal processing on graphs,”  IEEE Trans. Signal Process., vol. 61, no. 7, pp. 1644–1656, Apr. 2013.
  * [21] N. Perraudin and P. Vandergheynst,  “Stationary signal processing on graphs,”  IEEE Trans. Signal Process., vol. 65, no. 13, pp. 3462–3477, 2017\.
  * [22] B. Girault,  “Stationary graph signals using an isometric graph translation,”  in European Signal Process. Conf. (EUSIPCO), Aug 2015, pp. 1516–1520.
  * [23] J. Friedman, T. Hastie, and R. Tibshirani,  “Sparse inverse covariance estimation with the graphical lasso,”  Biostatistics, vol. 9, no. 3, pp. 432–441, 2008.
  * [24] X. Dong, D. Thanou, M. Rabbat, and P. Frossard,  “Learning graphs from data: A signal representation perspective,”  IEEE Signal Process. Mag., vol. 36, no. 3, pp. 44–63, 2019.
  * [25] S. Segarra M. Sevilla, A. G. Marques,  “Estimation of partially known gaussian graphical models with score-based structural priors,”  in 27th Intl. Conf. Artificial Intelligence and Statistics (AISTATS), 2024.
  * [26] H. Zou and T. Hastie,  “Regularization and variable selection via the elastic net,”  J. of the Roy. Statistical Soc.: Ser. B (Statistical Methodology), vol. 67, no. 2, pp. 301–320, 2005.
  * [27] R. Tibshirani,  “Regression shrinkage and selection via the lasso,”  J. of the Roy. Statistical Soc.: Ser. B (Statistical Methodology), vol. 58, no. 1, pp. 267–288, 1996.
  * [28] E. J Candes, M. B. Wakin, and S. P. Boyd,  “Enhancing sparsity by reweighted $\ell_{1}$ minimization,”  J. of Fourier Anal. and Appl., vol. 14, no. 5-6, pp. 877–905, 2008\.
  * [29] Y. Sun, P. Babu, and D. P. Palomar,  “Majorization-minimization algorithms in signal processing, communications, and machine learning,”  IEEE Trans. on Signal Processing, vol. 65, no. 3, pp. 794–816, 2017\.
  * [30] J. P. Boyle and R. L. Dykstra,  “A method for finding projections onto the intersection of convex sets in Hilbert spaces,”  in Advances in order restricted statistical inference, pp. 28–47. Springer, 1986.
  * [31] M. Li, D. Sun, and K.-C. Toh,  “A majorized ADMM with indefinite proximal terms for linearly constrained convex composite optimization,”  SIAM J. on Optimization, vol. 26, no. 2, pp. 922–950, 2016.
  * [32] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein,  “Distributed optimization and statistical learning via the alternating direction method of multipliers,”  Foundations and Trends® in Machine learning, vol. 3, no. 1, pp. 1–122, 2011.
  * [33] J. V. de M. Cardoso, J. Ying, and D. Palomar,  “Learning bipartite graphs: Heavy tails and multiple components,”  Adv. in Neural Inf. Process. Syst., vol. 35, pp. 14044–14057, 2022\.
  * [34] X. Cai, J. A. Bazerque, and G. B. Giannakis,  “Sparse structural equation modeling for inference of gene regulatory networks exploiting genetic perturbations,”  PLoS, Computational Biology, June 2013.
  * [35] M. Grant and S. Boyd,  “CVX: Matlab software for disciplined convex programming, version 2.1,” https://cvxr.com/cvx, Mar. 2014.
  * [36] J. V. de M. Cardoso and D. Palomar,  “Learning undirected graphs in financial markets,”  in Asilomar Conf. Signals, Systems, and Computers, 2020, pp. 741–745.
  * [37] L. Condat,  “Fast projection onto the simplex and the $\boldsymbol{l}_{\mathbf{1}}$ ball,”  Mathematical Programming, vol. 158, no. 1-2, pp. 575–585, 2016\.
  * [38] Y. Xu and W. Yin,  “A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion,”  SIAM J. on imaging sciences, vol. 6, no. 3, pp. 1758–1789, 2013\.