1

Characterization and Inference of Graph Diffusion
Processes from Observations of Stationary Signals

arXiv:1605.02569v4 [] 6 Jun 2017

Bastien Pasdeloup, Vincent Gripon, Grégoire Mercier, Dominique Pastor, and Michael G. Rabbat

Abstract—Many tools from the field of graph signal processing
exploit knowledge of the underlying graph’s structure (e.g., as
encoded in the Laplacian matrix) to process signals on the
graph. Therefore, in the case when no graph is available, graph
signal processing tools cannot be used anymore. Researchers have
proposed approaches to infer a graph topology from observations
of signals on its nodes. Since the problem is ill-posed, these
approaches make assumptions, such as smoothness of the signals
on the graph, or sparsity priors. In this paper, we propose a
characterization of the space of valid graphs, in the sense that
they can explain stationary signals. To simplify the exposition in
this paper, we focus here on the case where signals were i.i.d.
at some point back in time and were observed after diffusion
on a graph. We show that the set of graphs verifying this
assumption has a strong connection with the eigenvectors of
the covariance matrix, and forms a convex set. Along with a
theoretical study in which these eigenvectors are assumed to be
known, we consider the practical case when the observations
are noisy, and experimentally observe how fast the set of valid
graphs converges to the set obtained when the exact eigenvectors
are known, as the number of observations grows. To illustrate
how this characterization can be used for graph recovery, we
present two methods for selecting a particular point in this set
under chosen criteria, namely graph simplicity and sparsity.
Additionally, we introduce a measure to evaluate how much
a graph is adapted to signals under a stationarity assumption.
Finally, we evaluate how state-of-the-art methods relate to this
framework through experiments on a dataset of temperatures.

I. I NTRODUCTION
In many applications, such as brain imaging [1] and hyperspectral imaging [2], it is convenient to model the relationships
among the entries of the signals studied using a graph. Tools
such as graph signal processing can then be used to help understand the studied signals, providing a spectral view of them.
However, there are many cases where a graph structure is not
readily available, making such tools not directly appliable.
Graph topology inference from only the knowledge of
signals observed on the vertices is a field that has received a lot
of interest recently. Classical methods to obtain such a graph
are generally based on estimators of the covariance matrix
using tools such as covariance selection [3] or thresholding of
This was supported by the European Research Council under the European
Union’s Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement n◦ 290901, by the Labex CominLabs Neural Communications, and by
the Natural Sciences and Engineering Research Council of Canada through
grant RGPAS 429296-12.
B. Pasdeloup, V. Gripon, G. Mercier, and D. Pastor are with UMR CNRS
Lab-STICC, Télécom Bretagne, 655 Avenue du Technopole, 29280, Plouzané,
France. Email: {name.surname}@telecom-bretagne.eu.
M.G. Rabbat is with the Department of Electrical and Computer Engineering, McGill University, 3480 University Street, Montréal, H3A 0E9, Canada.
Email: michael.rabbat@mcgill.ca.

the empirical covariance matrix [4]. More recent approaches
make assumptions on the graph, and enforce properties such as
sparsity of the graph and/or smoothness of the signals [5]–[7].
A common aspect of all these techniques is that they
propose graph inference strategies that directly find a particular
topology from the signals based on some priors. Rather than
performing a direct graph inference, we explore an approach
that proceeds in two steps. First, we characterize the matrices
that may explain the relationships among signal entries. Then,
we introduce criteria to select a matrix from this set.
In this paper, we consider the case of stationary signals
[8]–[10]. These signals are such that their covariance matrix
has the same eigenvectors as the graph Fourier transform
operator. To simplify the exposition in this paper, we focus
here on the case of diffusion matrices, but the same ideas
and methods could work for general observations of stationary
signals on graphs. We assume that the signals were i.i.d. at
some point back in time. The relationships among entries of
the signals were then introduced by a diffusion matrix applied
a variable number of times on each signal. This matrix has
non-null entries only when a corresponding edge exists in
the underlying graph, and therefore is compliant with the
underlying graph structure, modeling a diffusion process on
it. Such matrices are referred to as graph shift operators [11],
[12], examples of which are the adjacency matrix or the graph
Laplacian. Under these settings, we address in this paper the
following question: How can one characterize an adapted
diffusion matrix from a set of observed signals?
To answer this question, we choose to focus in this paper on
a particular family of matrices to model the diffusion process
for the signals. Similar results can be obtained with other graph
shift operators by following the same development.
We show that retrieving a diffusion matrix from signals can
be done in two steps, by first characterizing the set of admissible candidate matrices, and then by introducing a selection
criterion to encourage desirable properties such as sparsity
or simplicity of the matrix to retrieve. This particular set of
admissible matrices is defined by a set of linear inequality
constraints. A consequence is that it is a convex polytope, in
which one can select a point by defining a criterion over the
set of admissible diffusion matrices and then maximizing or
minimizing the criterion.
We show that all candidate matrices share the same set of
eigenvectors, namely those of the covariance matrix. Along
with a theoretical study in which these eigenvectors are
assumed to be known, we consider the practical case when
only noisy observations of them are available, and observe
the speed of convergence of the approximate set of solutions

2

to the limit one, as the number of observed signals increases.
Two criteria for selecting a particular point in this set are
proposed. The first one aims to recover a graph that is simple,
and the second one encourages sparsity of the solution in the
sense of the L1,1 norm. Additionally, we propose a method to
obtain a diffusion matrix adapted to stationary signals, given
a graph inferred with other methods based on other priors.
This paper is organized as follows. First, Section II introduces the problem addressed in this article, and presents the
notions and vocabulary that are necessary for a full understanding of our work. Then, Section III reviews the work that has
been done in graph recovering from the observation of signals.
Section IV studies the desired properties that characterize
the admissible diffusion matrices, both in the ideal case and
in the approximate one. Section V introduces methods to
select an admissible diffusion matrix in the polytope based
on a chosen criterion. Then, in Section VI, these methods
are evaluated on synthetic data. Finally, in Section VII, a
dataset of temperatures in Brittany is studied, and additional
experiments are performed to establish whether current stateof-the-art methods can be used to infer a valid diffusion matrix.
II. P ROBLEM FORMULATION
A. Definitions
We consider a set of N random variables (vertices) of interest. Our objective is, given a set of M realizations (signals)
of these variables, to infer a diffusion matrix adapted to the
underlying graph topology on which unknown i.i.d. signals
could have evolved to generate the given M observations.
Definition 1 (Graph): A graph G is a pair (V, E) in which
V = {1, . . . , N } is a set of N vertices and E ⊆ V × V is a
set of edges. In the remainder of this document, we consider
positively weighted undirected graphs. Therefore, we make no
distinction between edges (u, v) and (v, u). We denote such
an edge by {u, v}. A convenient way to represent G is through
its adjacency matrix W:

αuv if {u, v} ∈ E
W(u, v) ,
; αuv ∈ R+ ; ∀u, v ∈ V .
0
otherwise
Graph shift operators are defined by Sandryhaila et al. [11],
[12] as local operations that replace a signal value at each
vertex of a graph with the linear combination of the signal
values at the neighbors of that vertex. The adjacency matrix
is an example of graph shift operator, as its entries are non-null
if and only if there exists a corresponding edge in E.
Graphs may have numerous properties that can be used as
priors when inferring an unknown graph. In this paper, we
are particularly interested in the sparsity of the graph, that
measures the density of its edge, and in its simplicity. A graph
is said to be simple if no vertex is connected to itself, i.e., if the
diagonal entries of the associated adjacency matrix are null.
These two properties are often desired in application domains.
In the general case, we consider graphs that can have selfloops, i.e., non-null elements on the diagonal of W.
A signal on a graph can be seen as a vector that attaches a
value to every vertex in V.
Definition 2 (Signal): A signal x on a graph G of N vertices
is a function on V. For convenience, signals on graphs are

represented by vectors in RN , in which x(i) is the signal
component associated with the ith vertex of V.
A widely-considered matrix that allows the study of signals
on a graph G is the normalized Laplacian of G.
Definition 3 (Normalized Laplacian): The normalized
Laplacian Ł of a graph G with adjacency matrix W is a
1
1
differential operator on G, defined by Ł , I − D− 2 WD− 2 ;
where D isPthe diagonal matrix of degrees of the vertices:
D(u, u) , v∈V W(u, v); ∀u ∈ V, and I is the identity matrix
of size N . Note that for Ł to be defined, D must contain only
non-null entries on its diagonal, which is the case when every
vertex has at least one neighbor.
An interpretation of this matrix is obtained by considering
the propagation of a signal x on G using Ł. By definition,
1
1
Łx = Ix − (D− 2 WD− 2 )x. Therefore, the normalized Laplacian models the variation of a signal x when diffused through
one step of a diffusion process represented by the graph shift
1
1
operator TŁ , D− 2 WD− 2 . More generally, in this paper, we
define diffusion matrices as follows:
Definition 4 (Diffusion matrix): A diffusion matrix T is a
symmetric matrix such that
• ∀u, v ∈ V : T(u, v) ≥ 0;
• λ1 = 1;
• ∀i, ∈ {2, . . . , N } : |λi | ≤ 1,
where λi are the eigenvalues of T, in descending order.
The idea behind these constraints is that we want to model a
diffusion process by a matrix. Such process propagates signal
components from vertex to vertex and consequently consists
of positive entries indicating what quantity of signal is sent
to the neighboring vertices. Enforcing all eigenvalues to have
their modulus be at most 1 imposes a scale factor, and has
the interesting consequence to cause the series Ti x i to be
bounded, for any signal x.
Note that by construction, the largest eigenvalue of TŁ is 1.
In our experiments, we will use this particular matrix TŁ to
diffuse signals on the graph. Other popular matrices could be
used instead to diffuse signals. For example, any polynome of
TŁ could be used [13]–[15].
B. Graph Fourier transform
One of the cornerstones of signal processing on graphs
is the analogy between the notion of frequency in classical
signal processing and the eigenvalues of the Laplacian. The
eigenvectors of the Laplacian of a binary ring graph correspond to the classical Fourier modes (see e.g., [16] for a
detailed explanation). The lowest eigenvalues are analogous
to low frequencies, while higher ones correspond to higher
frequencies. Using this analogy, researchers have successfully
been able to use graph signal processing techniques on nonring graphs (e.g., [17], [18]).
To be able to do so, the Laplacian matrix of the studied
graph must be diagonalizable. Although it is a sufficient,
but not necessary, condition for diagonalization, we only
consider undirected graphs in this article (see Definition 1),
for which the normalized Laplacian as defined in Definition 3
is symmetric. Note that there also exist definitions of the
Laplacian matrix when the graphs are directed [19].

3

To understand the link between diffusion of signals on the
graph and the notion of smoothness on the graph introduced in
Section II-C, we need to introduce the graph Fourier transform
[16], [20], that transports a signal x defined on the graph into
its spectral representation b
x:
Definition 5 (Graph Fourier transform): Let Λ =
(λ1 , . . . , λN ) be the set of eigenvalues of Ł, sorted by increasing value, and X = (χ1 , . . . , χN ) be the matrix of associated
eigenvectors. The graph Fourier transform of a signal x is the
projection of x in the spectral basis defined by X : b
x , X > x.
N
b
x is a vector in R , in which b
x(i) is the spectral component
associated with χi .
This operator allows the transportation of signals into a
spectral representation defined by the graph. Note that there
exist other graph Fourier transform operators, based on the
eigenvectors of the non-normalized Laplacian L , D − W or
on those of the adjacency matrix [11].
An important property of the normalized Laplacian states
that the eigenvalues of Ł lie in the closed interval [0, 2], with
the multiplicity of eigenvalue 0 being equal to the number of
connected components in the graph, and 2 being an eigenvalue
for bipartite graphs only [20]. We obtain that the eigenvalues
of TŁ lie in the closed interval [-1, 1], with at least one of
them being equal to 1. Also, since TŁ and Ł only differ by an
identity, both matrices share the same set of eigenvectors. If
the graph is connected then TŁ has a single eigenvalue equal
to 1, being associated with a constant-sign eigenvector χ1 :
s
D(i, i)
,
(1)
∀i ∈ {1, . . . , N } : χ1 (i) =
Tr(D)
where D is the matrix of degrees introduced in Definition 3,
and all other eigenvalues of TŁ are strictly less than 1.
Therefore, diffusing a signal x using TŁ shrinks the spectral
contribution of the eigenvectors of Ł associated with high
eigenvalues more than those associated with lower ones.
It is worth noting that since one of the eigenvalues of
TŁ is equal to 1, then the contribution of the associated
eigenvector χ1 does not change after diffusion. Therefore,
after numerous diffusion steps, (b
x(i))i∈[2;N ] become close to
null and x becomes stable on any non-bipartite graph. As a
consequence, we consider in our experiments signals that are
diffused a limited number of times.
C. Smoothness of signals on the graph
A commonly desired property for signals on graphs is
smoothness. Informally, a signal is said to be smooth on the
graph if it has similar entries where the corresponding vertices
are adjacent in the graph. In more details, given a diffusion
matrix T for a graph G, smoothness of a signal x can be
measured via the following quantity:
X
2
S(x) ,
T(u, v) (x(u) − x(v)) .
(2)
{u,v}∈E

From this equation, we can see that the lower S(x) is,
the more regular are the entries of x on the graph. When
using TŁ as a diffusion matrix, signals that are low-frequency,

i.e., that mostly have a spectral contribution of the lower
eigenvectors of the Laplacian, have a low value of S(x) and are
then smooth on the graph. As mentioned above, diffusion of
signals using TŁ shrinks the contribution of eigenvectors of Ł
associated with higher eigenvalues more than the contribution
of the ones associated with lower eigenvalues. Thence the
property that diffused signals become low-frequency after
some diffusion steps, and hence smooth on the graph. In
addition to seeing diffusion as a link between graphs and
signals naturally defined on them, this interesting property
justifies the assumption, made in many papers, that signals
should be smooth on a graph modeling their support [5]–[7].
D. Stationarity of signals on the graph
Considering stationary signals is a very classical framework
in traditional signal processing that facilitates the analysis of
signals. Analogously, stationary processes on graphs have been
recently defined to ease this analysis in the context of signal
processing on graphs [8]–[10].
A random process on a graph is said to be (wide-sense)
stationary if its first moment is constant over the vertex
set and its covariance matrix is invariant with respect to
the localization operator [10]. In particular, white noise is
stationary for any graph, and any number of applications
of a graph shift operator on such noise leaves the process
stationary. This implies that the covariance matrix of stationary
signals shares the same eigenvectors as this particular operator
(see Section III-D for details).
Diffusion of signals is a particular case of stationary processing. The example we develop in this article when studying
diffusion of signals through a matrix T can be generalized to
any stationary process and any graph shift operator, with only
few adaptations.
E. Problem formulation
Using the previously introduced notions, we can formulate
the problem we address in this paper as follows. Let X =
(x1 , . . . , xM ), xi ∈ RN , be a N ×M matrix of M observations,
one per column. Let Y = (y1 , . . . , yM ), yi ∈ RN , be a N ×M
unknown matrix of M i.i.d. signals; i.e., the entries Y(i, j) are
zero-mean, independent random variables. Let k ∈ RM
+ be an
unknown vector of M positive numbers, corresponding to the
number of times each signal is diffused before observation.
Given X, we aim to characterize the set of all diffusion
e such that there exist Y and k with:
matrices1 T
e k(i) y .
∀i ∈ {1, . . . , M } : xi = T
i

(3)

This framework can be seen as a particular case of graph
filters [12], containing only a monomial of the diffusion
matrix. From a practical point of view, this corresponds to
the setup where all signals are observed at a given time t,
but have been initialized at various instants t − k(i). More
generally, all polynomials of the diffusion matrix share the
same eigenvectors. The key underlying assumption in our work
1 Throughout this article we will denote recovered/estimated quantities using
a tilde.

4

is that each observation is the result of passing white noise
through a graph filter whose eigenvectors are the same as those
of the normalized Laplacian. Consequently, our approach can
be applied to any graph filter
∀i ∈ {1, . . . , M } : xi =

∞
X
ej y ,
(Ki )j T
i

(4)

j=0

for M sequences K1 , . . . , KM .
To summarize the following sections, we infer a diffusion
matrix in two steps. First, we characterize the convex set of
solutions using the method in Section IV. Then, we select a
point from this set using some criteria on the matrix we want
to infer. The strategies we propose are given in Section V.
III. R ELATED WORK
While much effort has gone into inferring graphs from
signals, the problem of characterizing the set of admissible
graphs under diffusion priors is relatively new, and forms the
core of our work. In this section we review related work on
reconstructing graphs from the observation of diffused signals
and make connections to the approach we consider. Additional
approaches exist but consider different signal models such as
time series [13], [21], band-limited signals [22] or combinations of localized functions [14], [15].
A. Estimation of the covariance matrix
As stated in the introduction, obtaining the eigenvectors of
the covariance matrix is a cornerstone of our approach. They
allow us to define a polytope limiting the set of matrices that
can be used to model a diffusion process.
h
i
Since the covariance matrix Σ , E XX> is not obtainable
in practical cases, a common approach involve estimating Σ
e
using the sample covariance matrix Σ:
e ,
Σ
where M(i, j) ,

1
M

1
(X − M)(X − M)> ,
M −1
M
P
k=1

(5)

X(i, k) is an N ×M matrix with each

row containing the mean signal value for the associated vertex.
An interesting property of this matrix is that its eigenvectors
converge to those of the covariance matrix as the number of
signals increases (see Section IV-B).
Other methods exist to infer a covariance matrix [23]–[26]
and may be interesting to consider in place of the sample
covariance matrix. Methods for retrieving a sparse covariance
matrix based on properties of its spectral norm are described
in [23] and [24]. However, these works do not provide any
information on the convergence rate of the eigenvectors of
their solutions to the eigenvectors of the covariance matrix,
as the number of signals increases. Similarly, [25] and [26]
retrieve covariance matrices that converge in operator norm
or in distribution. An intensive study of covariance estimation
methods could be interesting to find techniques that improve
the convergence of eigenvectors. This paper focuses on the use
of the sample covariance matrix.

B. Graphical lasso for graph inference
A widely-used approach to provide a graph is the graphical
lasso [27], which recovers a sparse precision matrix (i.e., ine under the assumption that the data
verse covariance matrix) Θ
are observations from a multivariate Gaussian distribution. The
core of this method consists in solving the following problem,


e = argmin Tr(ΣΘ)
e
Θ
− log det(Θ) + λkΘk1 ,
(6)
Θ≥0

e is the sample covariance matrix and λ is a regularwhere Σ
ization parameter controlling sparsity.
Numerous variations of this technique have been developed
[28]–[31], and several applications have been using graphical
lasso-based methods for inferring a sparse graph. Examples
can be found for instance in the fields of neuroimaging [32],
[33] or traffic modeling [34].
What makes this method interesting, in addition to its fast
convergence to a sparse solution, is a previous result from
Dempster [3]. In the covariance selection model, Dempster
proposes that the inverse covariance matrix should have numerous null off-diagonal entries. An additional result from
Wermuth [35] states that the non-null entries in the precision
matrix correspond to existing edges in a graph that is representative of the studied data.
Therefore, in our experiments, we evaluate whether considering the result of the graphical lasso as a graph makes it
admissible or not to model a diffusion process. However, when
considering (6), we can see that the method does not impose
any similarity between the eigenvectors of the covariance
matrix and those of the inferred solution. For this reason,
we do not expect this method to provide a solution that is
admissible in our settings.
Close to the graphical lasso, [36] and [37] propose an
algorithm to infer a precision matrix by adding generalized
Laplacian constraints. While this allows for good recovery of
the precision matrix, it proceeds in an iterative way by following a block descent algorithm that updates one row/column
per iteration. As for the graphical lasso, it does not force
the eigenvectors of the retrieved matrix to match those of the
covariance matrix, and therefore does not match our stationarity assumption. Interestingly, these methods could also be
mentioned in the next section, dedicated to smoothness-based
methods. In particular, [36] has pointed out that minimizing
e
the quantity Tr(ΣΘ)
promotes smoothness of the solution
when Θ is a graph Laplacian. Additionally, [38] promotes
sparsity of the inferred graph by applying a soft threshold to
the precision matrix, and shows that the solution matches a
smoothness assumption on signals.
C. Smoothness-based methods for graph inference
Another approach to recover a graph is to assume that the
signal components should be similar when the vertices on
which they are defined are linked with a strong weight in
W, thus enforcing natural signals on this graph to be lowfrequency (smooth). Using the definition of smoothness of
signals on a graph in (2), we can see that the smaller S(x),
the more regular the components of x on the graph.

5

A first work taking this approach has been proposed by
Lake and Tenenbaum [5], in which they solve a convex
optimization problem to recover a sparse graph from data
to learn the structure best representing some concepts. More
recently, Dong et al. [6] have proposed a similar method that
outperforms the one by Lake and Tenenbaum. In order to find a
graph Laplacian that minimizes S in (2) for a set of signals, the
authors propose an iterative algorithm that converges to a local
solution, based on the resolution of the following problem:
L∗ = arg min kX − Yk2F + α Tr(Y> LY) + βkLk2F
 L,Y
 Tr(L) = N
L(i, j) = L(j, i) ≤ 0, i 6= j
s. t.
,
PN

∀i ∈ {1, . . . , N } : j=1 L(i, j) = 0

(7)

where L∗ is the non-normalized Laplacian recovered, k · kF
is the Frobenius norm, Y is a matrix in RN ×M that can be
considered as a noiseless version of signals X, and α and β
are regularization parameters controlling the distance between
X and Y, and the sparsity of the solution.
Kalofolias [7] proposes a unifying framework to improve
the previous solutions of Lake and Tenenbaum, and Dong et
al., by proposing a better prior and reformulating the problem
to optimize over entries of the (weighted) adjacency matrix
rather than the Laplacian. An efficient implementation of his
work is provided in the Graph Signal Processing Toolbox [39].
His approach consists in rewriting the problem as an `1 minimization, that leads to naturally sparse solutions. Moreover,
the author has shown that the method from Dong et al. could
be encoded in his framework.
Graph inference with smoothness priors continues to receive
a lot of interest. Recently, Chepuri et al. [40] have proposed
to infer a sparse graph on which signals are smooth, using
an edge selection strategy. Finally, enforcing the smoothness
property for signals defined on a graph has also been considered by Shivaswamy and Jebara [41], where a method is
proposed to jointly learn the kernel of an SVM classifier and
optimize the spectrum of the Laplacian to improve this classification. Contrary to our approach, Shivaswamy and Jebara [41]
study a semi-supervised case, in which the spectrum of the
Laplacian is learned based on a set of labeled examples.
D. Diffusion based methods for graph inference
Recently we proposed a third approach to recover a graph
from diffused signals. In [42], we study a particular case of the
problem we consider here, namely when k is a known constant
vector. Let K denote the value in every entry of this vector. We
show in [42] that the covariance matrix of signals diffused K
times on the graph is equal to T2K . This implies that we need
to recover a particular root of the covariance matrix to obtain
T. In more details, if Y is a matrix of mutually independent
signals with independent entries, X = TK Y, and Σ is the
covariance matrix of X, we have:
i
h
i
h
>
Σ = E XX> = E TK YY> TK
= T2K ,
(8)
using the independence of Y and the symmetry of T.

Thanks to K being known, one could then retrieve a matrix
e by diagonalizing Σ, taking the 2K-square root of the obT
tained eigenvalues, and solving a linear optimization problem
to recover their missing signs. This reconstruction process was
illustrated on synthetic cases, where a graph G is generated,
and M i.i.d. signals are diffused on it using the associated
matrix TK to obtain X [42]. Experiments demonstrate that
e = T2K (which is the limit case when M grows
when using Σ
e = T.
to infinity), we can successfully recover T
However, this previous work has two principal limitations:
1) The number of diffusion steps k is constant and known,
which is a limiting assumption since in practical applications signals may be obtained after a variable,
unknown number of diffusion steps. In this work, we
remove this assumption. Taking the 2K-square root of
the eigenvalues of Σ is therefore no longer possible.
2) The number of observations M is assumed to be infinite
so that we have a perfect characterization of the eigenvectors of the covariance matrix. We also address this
assumption in this paper and show that the higher M ,
the closer the recovered graph to the ground truth.
Ongoing work by Segarra et al. [43], [44], initiated in
[45], takes a similar direction. The authors propose a two-step
approach, where they first retrieve the eigenvectors of a graph
shift operator, and then infer the missing eigenvalues based
on some criteria. They also study the case of stationary graph
processes, for which the covariance matrix shares the same
eigenbasis as the graph Fourier transform operator, and use
this information to infer a graph based on additional criteria.
However, while the characterization of the set of solutions
is identical to ours, our works differ in the matrix selection
strategy. Segarra et al. [43] focus on adjacency and Laplacian
inference, while we aim at recovering a matrix modeling
a diffusion process. Still, note that both of our works can
be easily extended to any graph shift operator, by setting
up the correct set of constraints. The authors of [43] solve
a slightly different problem, where they minimize the `1
norm of the inferred matrix under more constraints than ours,
which describe a valid Laplacian matrix. In particular, they
enforce the diagonal elements of the solution to be null,
thus considering graphs that do not admit self-loops. In more
details, they solve the following optimization problem:

PN
S = i=1 λi χi χ>
∗
i
S = arg min
kSk1 s. t.
,
S,λ1 ,...,λN
S∈S
(9)
where S∗ is the inferred graph shift operator, S is the set
of admissible solutions delimited by their constraints, and
χ1 , . . . , χN are the eigenvectors of the covariance matrix.
Contrary to their approach, we aim at inferring a matrix that
can be simple (see Section V-A) or sparse (see Section V-B),
rather than selecting a sparse matrix from the set of simple
matrices. Among other differences, we propose in Section V-C
a method to approximate the solution of any graph inference
strategy to make it match our stationary assumption on signals.
Our work also explores how the polytope of solutions can be
used to evaluate which graph, among a set of given graphs, is
the most adapted to given signals.

6

E. Other related work
Shahrampour and Preciado [46], [47] study the context of
network inference from stimulation of its vertices with noise.
However, their method implies a series of node knockout
operations that need to individually intervene on the vertices.
Also, we note that there exist methods that aim to recover
a graph from the knowledge of its Laplacian spectrum [48].
However, we do not assume that such information is available.
Finally, a recent work by Shafipour et al. [49] has started
to explore the problem of graph inference from non-stationary
graph signals, which is a direct continuation of the work
presented in this article and of the work by Segarra et al..
f2 , λ
f3 ) ∈ [−1, 1] × [−1, 1] (using a step
Figure 1: All pairs (λ
1 0 0 
−2
e
of 10 ) for which T = X 0 λf2 0 X > is an admissible
IV. C HARACTERIZATION OF THE SET OF ADMISSIBLE
DIFFUSION MATRICES

In this section, we show that the set of diffusion matrices
verifying the properties in Definition 4 is a convex polytope
delimited by linear constraints depending on the eigenvectors
of the covariance matrix of signals diffused on the graph. Then,
we study the impact of a limited number of observations on
the deformation of this polytope, due to imprecisions in the
obtention of these eigenvectors.
A. Characterization of the polytope of solutions
In the asymptotic case when M is infinite, the covariance
matrix Σ of the given signals X is equal to a (fixed) power K
of the diffusion matrix. Thus, under these asymptotic settings,
X can be obtained using Principal Component Analysis on X
[50]. In the more global case when k is a vector, the covariance
matrix of the signals is a linear combination of multiple powers
of T, and has therefore the same set of eigenvectors, since all
powers of a matrix share the same eigenvectors. This is also
the case when considering graph filters as in (4).
In more details, if we consider signals xi = Tk(i) yi , we have
the following development. We denote by X(i) the signal at
the ith column of X, and drop the constant factor and signals
mean from (5) for readability:
e =
Σ

M
X

>

X(i)X(i)

i=1

=

X X
k∈k

>

Tk Y(i)Y(i) Tk

>

i s.t.
k(i)=k



=

X
k∈k


X
>

>
Tk 
Y(i)Y(i)  Tk

(10)

i s.t.
k(i)=k



Σ =

X
k∈k

=

X


X
>

>
Tk EY 
Y(i)Y(i) Tk
i s.t.
k(i)=k

|{i, k(i) = k}| T2k ,

k∈k

which is a linear combination of various powers of T, all
having the same eigenvectors X .

f
0 0 λ
3

diffusion matrix (in red). The exact eigenvalues of the matrix
TŁ associated with W are located using a green dot.
Let us first consider the limit case when the eigenvectors
X of Σ are available. Given the remarks in Section II-B,
e we must find
to recover an acceptable diffusion matrix T,
f
f
e
eigenvalues Λ = (λ1 , . . . , λN ) such that:
e j) ≥ 0;
• ∀i, j ∈ {1, . . . , N }; j ≥ i : T(i,
f1 = 1;
• Let χ1 be the constant-sign eigenvector in X : λ
ei ∈ [−1, 1].
• ∀i ∈ {1, . . . , N } : λ
Note that these constraints are driven by the will to recover
a diffusion matrix as defined in Definition 4. If we were
considering the case of other graph shift operators, these
constraints would be different and would yield the definition
of different constraints in (13). As an example, diffusion using
a Laplacian matrix would imply the definition of constraints
that enforce the diagonal entries to be positive and off-diagonal
ones to be negative (see [43]). Similarly, aiming to recover
1
1
the diffusion matrix TŁ , D− 2 WD− 2 associated with the
normalized Laplacian would imply additional constraints.
To illustrate how these properties translate into a set
of admissible diffusion matrices, let us consider the randomly
3 × 3 symmetric adjacency matrix W =
 0.417 generated
0.302 0.186
0.302 0.147 0.346 . We compute its associated matrix TŁ
0.186 0.346 0.397
and corresponding eigenvectors X . This simulates a perfect
retrieval of the eigenvectors of the covariance matrix of signals
f2 , λ
f3 ) ∈ [−1, 1]×
diffused by TŁ on the graph. For all pairs (λ
[−1, 1] (using a step of 10−2 ), Fig. 1 depicts those that allow
the reconstruction of a diffusion matrix.
As we can see, the set of admissible matrices is convex
(but non-strictly convex), delimited by affine equations. To
characterize these equations, let us consider any entry of index
e we want to
(i, j) in the upper triangular part of the matrix T
recover. Since X is assumed
to
be
known,
by
developing
the
f

λ1 0 0
>
e , X 0 ... 0 X , we can write every
matrix product T
0

0 λf
N

f1 , . . . , λf
e j) as a linear combination of variables λ
entry T(i,
N
by developing the scalars in X . Let αij1 , . . . , αijN be the
factors associated with λ1 , . . . , λN for the equation associated
e j), i.e.:
with the entry T(i,
f1 + · · · + αijN λf
e j) = αij1 λ
T(i,
N .

(11)

7

As an example, let us consider a 3×3 matrix T with known
eigenvectors X . Using the decomposition of T, we can write
T(2, 3) as follows:
T(2, 3) = X (2, 1)X (3, 1) λ1 + X (2, 2)X (3, 2) λ2 + X (2, 3)X (3, 3) λ3 .
|
{z
}
{z
}
{z
}
|
|
α231
α232
α233

(12)
Enforcing the value of all entries of the matrix to be positive
thus defines the following set of N (N2+1) inequalities:
f1 +· · ·+αijN λf
∀i, j ∈ {1, . . . , N }; j ≥ i : αij1 λ
N ≥ 0 , (13)
where j ≥ i comes from the symmetry property.
Our problem of recovering the correct set of eigenvalues to
reconstruct the diffusion matrix thus becomes a problem of
selecting a vector of dimension N − 1 (since one eigenvalue
is equal to 1 due to the imposed scale) in the convex polytope
delimited by (13). Since the number of possible solutions
is infinite, it is an ill-posed problem. To cope with this
issue, one then needs to incorporate additional information
or a selection criterion to enforce desired properties on the
reconstructed matrix. To illustrate the selection of a point in
the polytope, Section V presents strategies based on different
criteria, namely sparsity and simplicity.
Note that the polytope is in most cases not a singleton. The
covariance matrix belongs to the set of admissible matrices,
along with all its powers, which are different if the covariance
matrix is not the identity. When additional constraints delimit
the polytope, for example when enforcing the matrices to be
simple, there exist situations when the solution is unique [43].

B. Impact of the use of the sample covariance matrix on the
polytope definition
The results discussed above use the eigenvectors X of the
limit covariance matrix Σ as M tends to infinity, directly obtained by diagonalizing the diffusion matrix of a ground truth
graph under study. Next, we study the impact of the use of the
e of controlled signals respecting
sample covariance matrix Σ
our assumptions on the estimation of these eigenvectors.
e instead of Σ, let us
To understand the impact of using Σ
again consider an example 3 × 3 matrix and the associated
polytope, as we did in Section IV-A. We generate a random
graph with N = 3 vertices by drawing the entries of its
adjacency matrix uniformly, and compute its matrix TŁ . Using
TŁ , we diffuse M i.i.d. signals (entries are drawn uniformly)
a variable number of times (chosen uniformly in the interval
[2, 5]) to obtain X. Then, we compute the sample covariance
e and its matrix of eigenvectors X
e.
matrix of X, Σ,
Fig. 2 depicts in white the ground truth polytope (i.e., the
one delimited by equations (13) using the eigenvectors of
Σ) and the polytopes associated with 10 different sample
e obtained from different realizations
covariance matrices Σ,
of X. From each of the eigenvector sets of these empirical
f2 , λ
f3 ) that
covariance matrices, we determine the pairs (λ
satisfy the criteria in Section IV-A. Then, we plot a histogram
of the number of occurrences of these valid pairs.
As we can see, the recovered polytope more accurately
reflects the true one as M increases. This coincides with the

fact that the empirical covariance matrix converges to the real
one as M tends to infinity.
In more details, we are interested in the convergence of
e =
the eigenvectors of the empirical covariance matrix X
e N } to those of the actual covariance matrix X .
{f
χ1 , . . . , χ
Asymptotic results on this convergence are provided by Anderson [51], which extends earlier related results by Girshick
fi be the vector of cosine
[52] and Lawley [53]. Let ei , X > χ
fi and all eigenvectors of the actual
similarities between χ
covariance matrix. Anderson [51] states that, as the number
of observations tends to infinity, entries in ei have a Gaussian
distribution with a known variance. In particular, when all
eigenvalues are distinct, the inner product between the ith (for
all i) eigenvector of the covariance matrix, χi , and the j th
fj , is asymptotically
(for all j) eigenvector of its estimate, χ
Gaussian with zero mean and variance
fj
λi λ
fj )2
(M − 1)(λi − λ

, λi 6= λj ,

(14)

fj is the
where λi is the eigenvalue associated with χi , and λ
fj . As a consequence, the variance
eigenvalue associated with χ
1
decreases like M
, and it also depends on the squared difference
fj . Additionally, [51] shows that the maximum
between λi and λ
likelihood estimate λei of λi (for all i) is
1 M −1 X
λj ,
(15)
λei =
Qi M
j∈Li

where Qi is the multiplicity of eigenvalue λi , and Li is the
set of integers {Q1 + · · · + Qi−1 + 1, . . . , Q1 + · · · + Qi },
containing all indices of equal eigenvalues. In the simple case
when all eigenvalues are distinct, (15) simplifies to
M −1
λei =
λi .
(16)
M
The eigenvalues of the empirical covariance matrix thus converge to those of the actual covariance matrix as M increases.
As M tends to infinity, ei thus tends to the ith canonical vector,
fi and χi . Additionally, [51]
indicating collinearity between χ
provides a similar result for the more general case when
eigenvalues may be repeated.
As we can see from (14), the convergence of the eigenvectors of the sample covariance matrix to the eigenvectors of the
true covariance matrix is impacted by the eigenvalues of the
matrix used to diffuse the signals. We know that diffusing a
signal K times using T is equivalent to TK x. When rewriting
this equation in the spectral basis, we obtain (X ΛK X > )x,
where X and Λ are the eigenvectors and eigenvalues of T.
As we can see, the power distributes on the eigenvalues, and
due to their location in the interval ]−1, 1] (with the noticeable
fj )2 in (14) gets
exception of bipartite graphs), the term (λi − λ
smaller, and M must grow to achieve the same precision.
To illustrate the impact of the number of diffusions on the
convergence of the polytope, let us consider the following
experiment. We generate 104 occurrences of random adjacency matrices of N = 10 vertices, by drawing their entries
uniformly in [0, 1], and enforcing symmetry. Then, for each
adjacency matrix, we compute the associated matrix TŁ , and

8

(a)

(b)

(c)

f2 , λ
f3 ) is valid, in the sense of the criteria in Section IV-A,
Figure 2: Histogram representing the number of times a pair (λ
e
when used jointly with X to recover a diffusion matrix. The ground truth polytope is represented by the inequality constraints
in white. Results obtained for 10 instances of X on the same graph, for M = 10 (a), M = 100 (b) and M = 1000 (c).

Inclusion ratio

1
M
M
M
M
M

0.8
0.6

= 105
= 104
= 103
= 102
= 10

0.4
5

10
K

15

20

Figure 3: Ratio of cases when the eigenvalues of the ground
truth matrix belong to the approximate polytope, as a function
of the number of diffusions K, for various quantities of signals
M . Tests were performed for 104 occurrences of random
adjacency matrices of N = 10 vertices.
for various values of K and of M , we diffuse M randomly
generated signals K times using TŁ . From the diffused signals,
we compute the eigenvectors of the sample covariance matrix,
and check if the eigenvalues of TŁ are located in the polytope
defined by these eigenvectors. Fig. 3 depicts the ratio of times
it is the case, for each combination of K and M .
The figure demonstrates that the number of diffusions has
some imporance in the process. Too small values of K encode
too little information on the diffusion matrix in the signals,
but values that are too high concentrate the eigenvalues too
much around 0. This corroborates that, as the eigenvalues
concentrate around 0, it is necessary to have a higher value of
M to achieve the same precision.
V. S TRATEGIES FOR SELECTING A DIFFUSION MATRIX
As stated in Section IV-A, inferring a valid diffusion matrix,
in the sense that it can explain the relationships among signal
entries through a diffusion process, reduces to selecting a point
in the polytope. Since it contains an infinite number of possible
solutions, one needs to introduce additional selection criteria
in order to favor desired properties of the retrieved solution.

Note that the polytope describes a set of diffusion matrices
e
as introduced in Definition 4. Given a diffusion matrix T
selected from the polytope, unless the degrees of the vertices
are known, there is no possibility in the general case to retrieve
the corresponding adjacency matrix. However, in the particular
case when the associated adjacency matrix is binary, one can
e at 0, setting its non-null entries to 1.
just threshold T
In this section, we first propose to illustrate the selection
of points in the polytope, using two criteria: simplicity of the
solution, and sparsity. In the first case, we aim at retrieving
a diffusion matrix that has an empty diagonal. In the second
case, we aim at recovering a sparse diffusion matrix. Additionally, we introduce a third method that performs differently
from the two other methods. Numerous graph inference techniques have been developed to obtain a graph from signals,
with various priors. While most of them do not require the
retrieved matrices to share the eigenvectors of the covariance
matrix, it may still be interesting to evaluate whether these
matrices are close enough to the polytope. If one can select
a point in the polytope that is close to the solution of a
given method, while keeping the properties enforced by the
associated priors, we obtain a new selection strategy. For this
reason, we introduce in this section a method to adapt the
solutions of other methods to stationary signals.
A. Selecting a diffusion matrix under a simplicity criterion
The first criterion we consider to select a point in the
polytope is simplicity of the solution. In other words, we want
e of eigenvalues that, jointly
to encourage the retrieval of a set Λ
with the eigenvectors X of the covariance matrix, produce a
diffusion matrix that has an empty diagonal. Such a matrix
represents a process that maximizes the diffusion of a signal
evolving on it, and does not retain any of its energy.
As shown in Section IV-A, since we are considering a
diffusion matrix as defined in Definition 4, the polytope of
solutions is defined by inequality constraints (13) that each
enforce the positivity of an entry in the matrix to recover. A
consequence is that if the matrix to be retrieved contains any

9

null entry, then the point we want to select lies on an edge or
a face of the polytope, since at least one inequality constraint
holds with equality. Enforcing simplicity of the solution is
therefore equivalent to selecting a point in the polytope that is
located at the intersection of at least N constraints. Using this
observation and the fact that the trace of a matrix is equal to the
sum of its eigenvalues, retrieving the eigenvalues that enforce
simplicity of the corresponding matrix reduces to solving a
linear programming problem, stated as follows:
N
X
f1 , . . . , λf
λ
=
arg
min
λi
N
λ1 ,...,λN
i=1

 (13)
∀i ∈ {1, . . . , N } : λi ∈ [−1, 1]
s. t.

λ1 = 1

(17)
,

where the two last constraints impose a scale factor.
Equation (17) is a linear program for which it is known that
polynomial-time algorithms exist. The main bottleneck of this
method is the definition of the N (N2+1) linear constraints in
(13), that are computed in O(N 3 ) time and space.
B. Selecting a diffusion matrix under a sparsity criterion
In many applications one may believe the graph underlying
the observations is sparse. Similar to the case when trying to
recover a simple graph, finding a sparse admissible solution
can be formulated as finding a point at the intersection of multiple linear constraints. To find a sparse solution, we seek the
set of admissible eigenvalues for which the maximum number
of constraints in (13) are null. This reduces to minimizing the
`0 norm of the solution, which is an NP-hard problem [54].
A common approach to circumvent this problem is to
approximate the minimizer of the `0 norm by minimizing
the `1 norm instead [55]–[57]. In our case, we use the L1,1
matrix norm, which is the sum of all entries, since they are all
positive. In this section, we adopt this approach and consider
again a linear programming problem as follows:
 λ1 0 0 
>
f1 , . . . , λf
0 ... 0
λ
=
arg
min
1
X
X > 1N
N
N
0 0 λN
λ1 ,...,λN

 (13)
(18)
∀i ∈ {1, . . . , N } : λi ∈ [−1, 1] ,
s. t.

λ1 = 1

normalized Laplacian). Let X Tm be the eigenvectors of Tm ,
and let X Σ be the eigenvectors of the covariance matrix.
The idea here is to consider Tm as if it were expressed in
the eigenbasis of X Σ , to check whether it belongs or not to
the polytope of admissible matrices. In other words, we want
to find a matrix Am such that Tm = X Σ Am X >
Σ . Using the
fact that X Σ forms an orthonormal basis, we have Am =
X>
Σ Tm X Σ . Unless X Tm and X Σ are the same, Am is not
necessarily a diagonal matrix. Therefore, Am lies in a space of
dimension N 2 , while the polytope is defined by N variables.
Let us call Λm = (λm1 , . . . , λmN ) the vector of elements
on the diagonal of Am . Since the polytope of admissible
diffusion matrices is defined in RN , Λm is the point of this
2
set that forms the best estimate for Am , defined in RN , after
dimensionality reduction. In other words, Λm is the orthogonal
projection of Am in the polytope. If Λm does not belong to
the polytope of admissible diffusion matrices characterized by
X Σ , then the method m provides a solution that does not
satisfy the conditions to be a diffusion process. To find the
point in the polytope that is the closest to Λm in the sense of
the Euclidean norm, we solve the following problem:
g
Λ
m = arg min kΛ − Λm k2
 Λ∈RN
 (13)
∀i ∈ {1, . . . , N } : Λ(i) ∈ [−1, 1]
s. t.

Λ(1) = 1

(19)
.

g
The solution to (19) gives us a set of eigenvalues Λ
m =
g
]
(λm1 , . . . , λmN ) that represent the best approximation of Am
when restricting the search to admissible diffusion
! matrices.
Therefore, the matrix Tf
m = XΣ

]
λ
0
m1 0
0 ...
0
0
0 λ]
mN

X>
Σ is the

adaptation of the solution of method m to stationary signals.
We measure the distance between Am , the projection of Tm
g
in the space defined by X Σ , and Λ
m , the closest point in the
polytope, as follows:


λg
0
0
m1
g
 0
...
0 
d(Am , Λ
,
(20)
m ) , Am −
]
0
0 λ
mN
F

where 1N is the vector of N entries all equal to one.

where k · kF is the Frobenius norm.
Fig. 4 summarizes provides a graphical illustration of the
various steps to compute this distance.

C. Adaptation of other strategies to stationary signals

VI. N UMERICAL EXPERIMENTS

The two methods introduced before consist in selecting a
point in the polytope, given simplicity or sparsity priors. In
this section, we take a different point of view. Many graph
inference techniques exist in the literature (see Section III),
all enforcing different properties of the graph that is retrieved.
However, most of them do not impose the eigenvectors of the
inferred solution to match those of the covariance matrix. The
idea here is to adapt these solutions to stationary signals.
To do so, let us consider an inference method m providing
an adjacency or a Laplacian matrix from a set of signals X =
(x1 , . . . , xM ). Let Tm be a diffusion matrix associated with
the inferred matrix (for example, the one derived from the

To be able to evaluate reconstruction performance of the
methods presented in Section V, we first need to design
experimental settings. This section introduces a generative
model for graphs and signals, and evaluates the performance
of the methods Simple and Sparse. The regularization method
introduced in Section V-C is evaluated as a means to select a
matrix representing best some given signals, among a set of
possible matrices. Section VII presents additional experiments
and comparisons with other methods from the literature on a
non-synthetic dataset of temperatures in Brittany.
Our experiments show that the Simple method succeeds in
recovering the ground truth matrix from signals diffused on

10

Definition 7 (Erdős-Rényi graph): An Erdős-Rényi graph of
parameter P is a graph where each edge exists with probability
P independently from each other. The adjacency matrix W of
such a graph is defined by:

1 with probability P
W(i, j) ,
.
(22)
0 with probability 1 − P

Figure 4: Correction of the result of an inference method m
to match the stationarity hypothesis on the observed signals.
The eigenvalues of the result of m are expressed in the space
defined by X Σ as a matrix Am . Then, Am is approximated
by Λm , its orthogonal projection in the space of the polytope.
The closest point in the polytope (in the sense of the Euclidean
g
norm), Λ
m , is then found by solving (19). Finally, the distance
between Am and its estimate in the polytope is given by the
measure in (20). This corresponds to the norm of the vector
in green.

it, provided that the number of signals is high enough. The
Sparse method, while not being able to retrieve the ground
truth matrix, infers a matrix that has a lower L1,1 norm than
the matrix yielding the polytope. Additionally, we show that
the regularization method allows the selection of the ground
truth diffusion matrix from a set of candidate matrices, even
for a small number of signals. Finally, comparison with other
inference methods on a dataset of temperatures show that
the methods Simple and Sparse return the best solutions with
respect to their objectives, and that the regularization strategy
applied to methods favorizing smoothness of signals yields a
diffusion matrix on which signals are relatively smooth.

Once a graph is generated using the model presented above,
and given a number M of signals to produce, signals verifying
our settings are created as follows:
1) Create Y a N × M matrix with i.i.d. entries. We denote
by Y(i) the ith column of Y. In these experiments,
entries of Y are drawn uniformly.
2) Create k a vector of M i.i.d. integer entries, comprised in
the interval {1, . . . , 10}. These values are chosen not too
high in reaction to a remark in Section II-B, not to obtain
signals that are already stable. In these experiments,
entries of k are also drawn uniformly.
3) Compute the diffusion matrix TŁ associated with G.
4) Create X a N × M matrix of signal as follows: ∀i ∈
k(i)
{1, . . . , M } : X(i) , TŁ Y(i).
Once these four steps are performed, the objective becomes:
given X and some criteria on the graph to retrieve, infer
e for the diffusion matrix of the signals. To
an estimate T
summarize the previous sections, we proceed as follows:
e , an estimate for the eigenvectors of the diffusion
1) Find X
matrix. Here, this is done by computing the eigenvectors
of the sample covariance matrix.
e , compute the constraints in (13) that define the
2) Using X
polytope of solutions.
3) Select a point from the polytope, using one of the
strategies in Section V.

A. Generative model for graphs and signals
In our experiments, we consider randomly generated graphs,
produced by a random geometric model. These are frequently
used to model connectivity in wireless networks [58].
Definition 6 (Random geometric graph): A random geometric graph of parameter R is a graph built from a set of N
uniformly distributed random points on the surface of a unit
2-dimensional torus, by adding an edge between those being
closer than R according to the geodesic distance d(i, j) on
the torus. We then add a weight on the existing edges that
is inversely proportional to the distance separating the points.
Here, we choose to use the inverse of d(i, j). The adjacency
matrix W of such a graph is defined by:
 1
if d(i, j) < R and i 6= j
d(i,j)
.
(21)
W(i, j) ,
0
otherwise
Note that, by construction, such graphs are simple and
relatively sparse. Therefore, we expect methods using such
selection criteria to be able to retrieve them.
In some of our experiments, we also consider random graphs
generated by an Erdős-Rényi model [59], in which two vertices
are linked with a given probability independently from each
other. Such graphs are defined as follows:

B. Error metrics
To be able to evaluate the reconstruction error for our
techniques, we use multiple metrics. Let T be the ground truth
diffusion matrix, with eigenvalues Λ = (λ1 , . . . , λN ), and let
e be the one that is recovered using the assessed technique,
T
f1 , . . . , λf
e = (λ
with eigenvalues Λ
N ).
The first metric we propose is the mean error per reconstructed entry (MEPRE):
e ,
MEPRE(T, T)

1
N

e
T
T
−
e F
kTkF
kTk

.

(23)

F

This quantity measures the mean error for all entries in the
e using
reconstructed matrix, where we first normalize T and T
their Frobenius norm k · kF to avoid biases related to scale.
The second metric we propose is the reconstruction error of
the powered retrieved eigenvalues (REPRE). We define it as
the Euclidean distance between the Kth-power of the ground
e for
truth vector of eigenvalues ΛK and the recovered ones Λ,
the best value of K possible:
1
K∈R N

e , min
REPRE(Λ, Λ)

e
ΛK
Λ
−
K
e ∞
kΛ k∞
kΛk

.
2

(24)

11

Here, the normalization using k · k∞ comes from the constraint in Section IV-A that the highest eigenvalue should
be equal to 1. Therefore, it imposes a scale on the set of
eigenvalues. Also, we divide the error by N to make it
independent of the number of vertices.
Since in the limit case Σ = T2K , for an unknown K ∈ R,
then the algorithms should be able to recover at least a power
of T. Indeed, there is absolutely no way to distinguish for
example one step of diffusion of a signal x using T2K(i.e.,
2

T2K x) and two diffusion steps of x using TK (i.e., TK x).
A consequence is that, if the algorithm cannot fully retrieve
T, a power of T should also be an acceptable answer.
These two metrics provide information on the ability of
methods to infer a diffusion matrix that is close to the ground
truth one, or one of its powers. More classical metrics are also
considered to evaluate whether the most significant entries
of the inferred matrices correspond to existing edges in the
ground truth graph. Since such metrics are defined for binary,
relatively sparse graphs, we evaluate them on thresholded
versions of the inferred matrix. Entries of the inferred matrix
that are above this threshold t are set to 1, and others are set to
0. The resulting matrix is denoted in the following equations
e t . To find the optimal value, we perform an exhaustive
as T
search among all possible thesholds, and keep the one that
maximizes the F-measure, defined below.
The first metric we consider is the Precision, measuring the
fraction of relevant edges among those retrieved.
et) ,
precision(T, T

et (i,j)>0 and T(i,j)>0}
#{(i,j)|T
et k0,1
kT

,

(25)

where k · k0,1 denotes the L0,1 matrix norm, that counts the
number of non-null entries.
A second metric we consider is the Recall, that measures
the fraction of relevant edges effectively retrieved.
et) ,
recall(T, T

et (i,j)>0 and T(i,j)>0}
#{(i,j)|T
kTk0,1

,

(26)

Both metrics are often combined into a single one, called
F-measure, that can be considered as a harmonic mean of
precision and recall.
et) , 2
F-measure(T, T

et )·recall(T,T
et )
precision(T,T
et )+recall(T,T
et )
precision(T,T

.

(27)

Note that in practical cases, the optimal threshold is not
available, and depends on a desired sparsity of the inferred
matrix. In the following experiments, we show the compromise
between true positive edges and false positive edges for all
possible thresholds using ROC curves.
C. Performance of the Simple method
In the situation when the eigenvectors are not available,
due to a limited number of signals, recovery methods must
use estimated eigenvectors. Linear programming problems
introduced in Section V must then be solved on a polytope
defined by noisy eigenvectors. We have previously shown in
Section IV-B that increasing the number of signals allows
this approximate polytope to be more precise. The following

experiments evaluate the quality of the solutions retrieved by
the two methods introduced in Section V-A and Section V-B.
Fig. 5 illustrates the convergence of the Simple method to a
solution, when the number of signals increases. In this experiment, we generate 1000 random geometric graphs (N = 10,
R = 0.6) and, for each of them, we diffuse M ∈ {10i }1≤i≤6
signals using TŁ , as described in Section VI-A. For each
configuration, we retrieve a diffusion matrix by solving the
problem in (17) using the CVX [60] package for MATLAB
[61], with default parameters. Then, we compute the mean
errors, and measure the distance to the ground truth solutions
in terms of trace value, which is the objective function in (17):


e − Tr(TŁ ) ,
e , 1 Tr(T)
(28)
diff simple (TŁ , T)
N
e , 1 Tr(T)
e
which in our case simplifies to diff simple (TŁ , T)
N
since by construction Tr(TŁ ) = 0.
These results show that both error measurements decrease
— except for very low values of N due to the high variance
— as M increases. As the approximate polytope converges
to the ground truth one, the solution of (17) converges to the
ground truth matrix TŁ used to produce the signals. This is
confirmed by performance metrics as well as the ROC curves,
which indicate that the ground truth edges are recovered more
successfully as the number of signals increases.
As stated in Section VI-A, matrices generated using the
random geometric model are by construction simple, and
have therefore null traces. When defining the polytope using
inequalities in (13), we enforce the positivity of all entries of
the admissible matrices. Therefore, TŁ is a matrix for which
eigenvalues lie on a plane where λ1 + λ2 + λ3 = 0. While
the optimal solution may not be unique (see below), counterexamples to this uniqueness must respect particular constraints
and are very unlikely to happen for random matrices as well
as for non-synthetic cases. As a consequence, minimizing the
sum of the eigenvalues as an objective function enforces the
retrieval of the correct result in nearly all cases. This implies
that the error measurements will certainly converge to 0 as M
e
grows to infinity. This can easily be verified by replacing X
by X — the eigenvectors of TŁ — in (13) for the resolution
of (17).
As stated above, the solution is not necessarily unique.
Multiple matrices with a minimum trace can be found when
there exists a frontier of the polytope along which all points
have a sum that is minimal among all admissible vectors of
eigenvalues. As an example, let us consider the 8×8 Hadamard
matrix as a matrix of eigenvectors X : When defining the
polytope associated with these eigenvectors using constraints
in (13), we obtain for T(2, 6) the following constraint:
λ2 + λ3 + C ≥ 0 ,

(29)

where C is some value that does not depend on λ2 and λ3 .
As a consequence, any point located on this particular plane
corresponds to a matrix with identical sum of eigenvalues,
leading to the same trace. When considering random matrices
or non-synthetic cases, the case when there exists such a plane
that is aligned with the objective is very unlikely.

12

1

1
MEPRE
REPRE
diffsimple

0.8
Performance

Error

0.06
0.04

True positives

0.08

0.9

0.8
Recall
Precision
F-measure

0.02
0 1
10

102

103

104

105

106

0.7 1
10

102

103

M

(a)

104

105

0.6

M
M
M
M
M
M

0.4
0.2

106

0

0

M

(b)

0.2

0.4
0.6
False positives

= 106
= 105
= 104
= 103
= 102
= 10

0.8

1

(c)

Figure 5: Image (a) depicts the mean MEPRE and REPRE measurements of the solutions retrieved by the Simple method,
for M ∈ {10i }1≤i≤6 signals. Additionally, diff simple shows the distance to the ground truth solutions in terms of trace value,
which is the objective function of problem presented in (17). Image (b) shows the results in terms of edge reconstruction by
e t . Finally, image (c) depicts the ROC
studying the recall, precision and F-measure using the binarized thresholded matrix T
curves, that show the compromise between true positive edges and false positive edges when varying the threshold. All tests
were performed for 1000 occurrences of random geometric graphs with parameters N = 10 and R = 0.6.

D. Performance of the Sparse method
In this section, we perform the same experiment as for the
Simple method, but solving the problem in (18) instead of
(17). Here, since we minimize the L1,1 norm as an objective
function, we measure the mean difference between the sparsity
of the retrieved matrices and the sparsity of the ground truth
matrices TŁ , computed as follows:


e 1,1 − kTŁ k1,1 .
e , 1 kTk
(30)
diff sparse (TŁ , T)
N2
For space considerations, the results are not detailed here.
Contrary to the method for selecting a simple graph in
Section V-A, we have observed that the error and performance
measurements stay approximately constant for all values of M
(MEPRE ≈ 0.1, REPRE ≈ 6 × 10−2 , diff sparse ≈ −2 × 10−2 ,
F-measure ≈ 0.76). Similar results were obtained for additional
experiments conducted on random geometric graphs with
different values of R to assess different levels of sparsity.
The results suggest that the method fails at recovering the
matrix that was used for diffusing the signals, even when
the number of observations is high. However, the negative
difference indicates that the method does not fail at recovering
a sparse graph, in the sense of the L1,1 norm. Graphs generated using the random geometric model, while being sparse
by construction, are not necessarily the sparsest within the
associated admissible set, especially when considering the L1,1
norm. This implies that, as the number of signals increases,
the method in fact converges to the sparsest solution in the
polytope, although it is in most cases not the matrix we started
e by X — the eigenvectors of TŁ — in
from. Replacing X
(13) for the resolution of (18) confirms that there exist sparser
solutions than the ground truth graph.
Note that in their work, Segarra et al. [43] also use the L1,1
norm minimization as an objective, and suceed in retrieving
the ground truth graph. This is the case because they have
additional constraints that enforce the solution to be simple.
Therefore, the set of solutions is a lot smaller, and the solution
is most likely unique.

E. Impact of the parameters
The quality of the solutions inferred by the methods assessed above mostly depends on how close the eigenvectors
of the sample covariance matrix are from the ground truth
ones. Noticing that the density of eigenvalues in the interval
[−1, 1] increases with N for any diffusion matrix, it follows
that larger graphs have eigenvalues pairwise closer than for
smaller graphs. In this respect, Section IV-B tells us that the
number of signals necessary for a precise estimate of the
covariance matrix needs to be higher for larger graphs. This
raises the question of the scalability of this method.
For a fixed value of M = 105 signals, we study the
performance of the Simple and Sparse methods on random
geometric graphs of orders ranging from N = 10 to N = 100
vertices. For each value of N , we set R ∝ √1N so that the value
of N does not impact the average neighborhood of each vertex.
This value is chosen in accordance with the experiments
performed earlier in this section. Fig. 6 depicts the F-measure
performance measurements obtained for the Simple and Sparse
methods, for 1000 occurrences of random geometric graphs
with the previously described settings. Additionally, we plot
in this figure the results obtained for other families of graphs,
namely Erdős-Rényi graphs with P ∝ logNN , and the ring
graph of N vertices.
These measurements confirm that the inference methods
need a large quantity of signals to work properly. Additionally, it appears that the family of graphs has an impact on
the reconstruction performance. The ring graph has repeated
eigenvalues, which has a strong impact on the convergence of
the eigenvectors of the sample covariance matrix to those of
the real covariance matrix (see Section IV-B and [51]). However, this is a marginal case, since real-weighted graphs almost
surely have distinct eigenvalues, as illustrated by the studies on
random geometric and Erdős-Rényi graphs. Another parameter
that has importance on these experiments is the number of
diffusions of signals before observations, represented by the
vector k. This has been illustrated in Section IV-B.

13

1

1

F-measure

0.8
0.6

(RG)
(RG)
(ER)
(ER)
(ring)
(ring)

Success ratio

Simple
Sparse
Simple
Sparse
Simple
Sparse

0.8
0.6
0.4

50

0.4
0.2
0

20

40

60

80

100

100
M

150

200

Figure 7: Ratio of times when the diffusion matrix Ti was chosen by the algorithm as the most adapted to signals Xi among
a set of 20 possible diffusion matrices, for i ∈ {1, . . . , 20}.
Mean results for 100 iterations of the experiment.

N

Figure 6: F-measure scores obtained for the Simple and Sparse
methods, for various families of graphs and for M = 105 signals, as a function of the graph order. Tests were performed on
1000 occurrences of each family of graph, for each technique.
F. Application of regularization to graph hypothesis testing
To evaluate the practical interest of the regularization strategy introduced in Section V-C, let us consider the situation
where some signals are observed, and various diffusion matrices are provided by inference methods to explain these signals.
The objective is to determine which of the proposed solutions
matches the signals best under a stationarity assumption.
In the following experiment, we proceed as follows: let
(T1 . . . T20 ) be a set of 20 diffusion matrices corresponding
to graphs of N = 10 vertices, equally divided into random
geometric graphs (with R drawn uniformly in [0.2, 0.6]) and
Erdős-Rényi graphs (with P drawn uniformly in [0.2, 0.6]).
For each of these matrices, let us diffuse M random signals
as detailed in Section VI to obtain observations (X1 . . . X20 ),
where Xi is the set of signals obtained after diffusion by
Ti . From these sets, we can compute the eigenvectors of the
f1 . . . X
g
sample covariance matrices (X
20 ).
f
For each matrix of eigenvectors X i , i ∈ {1, . . . , 20}, and
for each diffusion matrix Tj , j ∈ {1, . . . , 20}, we compute
fi and the
the distance between the polytope yielded by X
projection of Tj in the space of the polytope (see Section V-C)
using (20). The graph that minimizes the distance is then
selected as the most appropriate. Fig. 7 depicts the ratio of
times Ti is selected as the most appropriate diffusion matrix
when considering signals Xi , for various values of M .
The results show that the regularization strategy selects the
matrix used to diffuse the signals in most cases, even when
M is low. Additional experiments were performed for larger
graphs, and similar results were observed. Also, increasing the
number of signals eventually leads to a selection of the correct
diffusion matrix in all cases. This experiment illustrates that
the regularization strategy introduced in Section V-C can be
successfully used to select the graph that is the most adapted
to given signals among a set of candidates.
An interesting direction for future work includes evaluation
of the performance of the method when considering selection
of the most adapted matrix from noisy versions of Ti .

VII. E VALUATION OF INFERENCE METHODS ON A DATASET
The two methods introduced in Section V-A and Section V-B present solutions to infer a graph from signals, while
ensuring that it is compliant with our diffusion prior. While
other methods from the literature do not clearly impose this
prior, evaluating whether they provide solutions that match
a diffusion assumption is interesting, as it would provide
additional selection strategies for admissible diffusion matrices
if it is the case. This section explores the application of
the regularization strategy in Section V-C to the method of
Kalofolias [7], and shows that it provides matrices that do not
belong to the polytope of solutions. The closest point in the
polytope is considered, and evaluation on a dataset shows that
the result has interesting similarities with the original matrix.
Throughout this section, we study an open dataset2 of
temperature observations from 37 weather stations located in
Brittany, France [62]. Our inference methods, as well as other
existing methods, are evaluated on this dataset in terms of
sparsity, trace of the solution, and smoothness.
A. Detailed evaluation of the method from Kalofolias
The method from Kalofolias has two major qualities: it recovers a graph in a very short amount of time, and encourages
smoothness of the solution, which can be a desirable property.
To evaluate whether the retrieved solution happens to match
a diffusion process, let us consider the following experiment:
1) Let G be a random geometric graph of N = 10
vertices (R = 0.6), and let TŁ the diffusion matrix
associated with its normalized Laplacian. Using this
matrix, we diffuse M = 106 i.i.d. signals as presented
in Section VI-A to obtain a matrix X. Using Principal
e , an estiComponent Analysis on X [50], we obtain X
mate for the eigenvectors of TŁ . This set of eigenvectors
yields a polytope of admissible solutions.
2) Then, we use the method from Kalofolias to infer a
graph G K from X, and compute the associated matrix
TŁK . Since the log method from Kalofolias depends
on parameters α and β, we keep the minimal distance
obtained for values of α and β ranging from 0.01 to 2,
with a step of 10−2 . Equation (20) gives us the distance
between the polytope and the inferred solution.
2 In

http://data.gouv.fr.

Occurrences

14

400

Kalofolias (dK )
Random (dR )

200
0
0

0.5
1
1.5
Distance to polytope

2

Figure 8: Number of times a distance to the ground truth polytope was observed using either the method from Kalofolias
[7] (dK ), or a random geometric graph (dR ). Distances are
grouped in bins of size 10−2 . Tests were performed for 105
occurrences of graphs per method, with M = 106 signals.
3) Additionally, we generate a random geometric graph G R
(independent from the ground truth one) from X, using
the same settings as for G (N = 10, R = 0.6). This
gives us a baseline of how close a random graph with
the same edges distribution can be to the ground truth
one, and gives information on whether the results of
Kalofolias are closer to the ground truth than a random
matrix. Again, (20) measures the distance between the
polytope and the associated matrix TŁR .
We perform these three steps for 105 occurrences of random
geometric graphs. Let dK be the vector of distances to the
polytope obtained for each ground truth graph using the
method of Kalofolias, and dR the vector of distances to the
polytope for the baseline random graphs. In Fig. 8, we plot a
histogram of the number of times each distance was observed.
From these results, a first observation is that neither the
methods from Kalofolias nor the random method ever returned
a graph that was located in the polytope of solutions. Two
direct interpretations of this result can be made: first, it implies
that the set of admissible matrices per ground truth graph
is small relatively to the set of random graphs. Second, it
implies that the method from Kalofolias does not succeed in
recovering a graph that matches diffusion priors on the signals.
Mann-Whitney U test [63] on dK and dR shows that the
distributions differ significantly (U = 9.9813×107 , P < 10−5
two-tailed). This implies that the results obtained with the
method from Kalofolias are most of the time closer to an
admissible matrix than random solutions. This observation
can be explained by the remarks in Section II-C. Diffusion
of signals on a graph tends to smoothen them, as the low
frequencies are attenuated slower than higher ones. Since the
method of Kalofolias retrieves a graph on which signals are
smooth, the observation that it provides solutions that are
closer to the polytope than random solutions is quite natural.
The question is then whether the closest point to the
retrieved solution in the polytope has interesting properties.
Let us evaluate this solution on the dataset of temperatures.
Fig. 9 depicts the 10% most significant connections in the
adjacency matrix of the graph G K retrieved by Kalofolias, as
well as those of the matrix associated with the closest point in
e=X
eΛ
eX
e > , where Λ
e is the solution of (19).
the polytope, T

(a)

(b)

Figure 9: Most significant connections in the adjacency matrix
of the graph G K retrieved by the method from Kalofolias (a),
and most significant connections from the matrix associated
e (b).
with the closest point in the polytope, T
The method from Kalofolias retrieves a matrix that has
stronger connections between stations with similar locations.
There is a strong connectivity among stations located on the
south coast of Brittany. Stations located more in the land also
tend to be linked to close inland stations. The regularized
matrix appears to keep these properties: the strong links on the
coasts still appear, and the result also still gives importance
on the coastal versus inland aspect of the stations. Still,
differences can be seen, as the regularized matrix appears to
give more importance to the relations between stations on the
north coasts. Such relations also exist in the original matrix,
but are not depicted due to the threshold.
When computing the total smoothness of the signals with
both matrices, we obtain that the solution from Kalofolias
has a higher value of smoothness (see Fig. 10) than the
closest point in the polytope. This implies that the signals are
smoother on the approximate matrix than on the one recovered
by the method of Kalofolias. This may seem counter-intuitive,
since the solutions of the method by Kalofolias are not restricted to the polytope. However, the method from Kalofolias
imposes inference of a matrix with an empty diagonal, which
is not the case of the approximate one. These measurements,
in addition to those below, suggest that inferring a graph using
the method from Kalofolias, and considering the closest point
in the polytope, is an interesting method to infer a valid graph
on which signals are smooth.
B. Evaluation of graph inference methods on the dataset
We have proposed in Section V-C a technique to find a
valid matrix in the polytope that approximates the solution
of any method. Therefore, it is possible to evaluate all methods in terms of properties, such as L1,1 sparsity, trace, or
smoothness. Since all methods do not impose the same scale
on the inferred matrices, these quantities are computed for the
inferred diffusion matrices after normalization such that their
first eigenvalue equals one, as in the constraints in Definition 4.
When applying our methods Simple (Section V-A) and
Sparse (Section V-B), as well as those of Kalofolias [7],
Segarra et al. [43] and the graphical lasso [27] on the dataset
of temperatures, we have obtained the results in Fig. 10.

15

Simple
Sparse
Kalofolias [7]
Kalofolias closest
Segarra et al. [43]
Segarra et al. closest
Graphical lasso [27]
Graphical lasso closest

polytope
X
X
0.0313
X
0.0062
X
1.3730
X

L1,1
36.9974
36.9971
36.9979
36.9974
36.9993
36.9974
35.3539
36.9984

Tr
0.0013
0.9093
0
0.0298
1.97 × 10−5
0.0046
13.4977
13.3584

S(X)
0.0551
0.0585
0.0751
0.0548
0.0245
0.0551
32.8421
0.0335

Figure 10: Sparsity, trace and smoothness obtained for the
dataset of temperatures. Elements in bold denote the method
performing best among those that return a solution located in
the polytope. If a method provides a solution that does not
belong to the polytope, the distance to the closest point is
indicated in the first column. The last column
P indicates the
total smoothness for all signals, i.e., S(X) , x∈X S(x).
First, we notice that the method from Segarra et al. [43]
returns a matrix that is at a distance of 0.0062 from the
polytope. As the polytope description is the same for their
method and for ours, we would expect this distance to be
0. This small difference comes from their implementation. In
order to keep their equality constraints enforcing the elements
in the polytope to have an empty diagonal, while coping with
the noise in the eigenvectors, they allow small deviations from
the polytope. They do not return a matrix S∗ that shares the
eigenvectors of the covariance matrix, but a matrix Ŝ∗ such
that kS∗ − Ŝ∗ kF ≤ ε. Here, experiments were performed for
ε = 10−3 . For this reason, the matrix they return is located
slightly outside of the polytope of solutions. When considering
the closest point to this result in the polytope, it appears to be
very close to the solution returned by the Simple method.
As expected, the Sparse method recovers the matrix with
the lowest L1,1 norm. It is also interesting to remark that
the projection of the solution of the graphical lasso on the
polytope is smoother than the projection of the solution
obtained by the method from Kalofolias. This echoes the
remark in Section III-B that minimization of the quantity
e
Tr(ΣΘ)
in (6) tends to promote smoothness of the signals
on the graph when Θ is a Laplacian matrix, which appears
to be encouraged by the regularization algorithm. Note that
the method from Kalofolias infers a graph which projection
on the polytope gets the second best smoothness score, while
having a small trace. On the other hand, the solution infered
by the graphical lasso appears to have most of its energy on
the diagonal entries. This is confirmed by the traces of the
matrices, both for the original solution and its approximate in
the polytope. Therefore, these two solutions provide interesting
ways to find a graph on which stationary signals are smooth,
with different simplicity assumptions.
VIII. C ONCLUSIONS
In this article, we have proposed a method for characterizing
the set of matrices that may be used to explain the relationships
among signal entries assuming a diffusion process. We have
shown that they are part of a convex polytope, and have
illustrated how one could choose a point in this particular set,
given additional selection criteria such as sparsity of simplicity

of the graph to infer. Finally, we have shown that most of
other existing methods do not infer matrices that belong to
the polytope of admissible solutions for stationarity signals,
and have introduced a method to consider the closest valid
matrix. An experiment was performed to illustrate that this
particular method can be useful for graph hypothesis testing.
Future directions based on this work are numerous. First
of all, reviewing the covariance estimation techniques is an
interesting direction, as obtention of the eigenvectors of the
covariance matrix is a cornerstone of our approach, and some
techniques may provide such information more precisely than
the sample covariance. We could also explore new strategies
to select a point in the polytope, for example by enforcing
the reconstruction of a binary matrix. Another interesting
direction would then be to propose selection strategies that
do not imply the full definition of the N (N2+1) constraints
defining the polytope. Finally, our immediate next work will
be to complement our experiments on graph hypothesis testing,
considering noisy versions of the candidate diffusion matrices.

ACKNOWLEDGEMENTS
The authors would like to thank the reviewers of previous
versions of this paper, whose remarks helped improving the
quality of our work. Also, we would like to thank Xiaowen
Dong and Santiago Segarra for kindly providing their codes, as
well as Benjamin Girault for sharing his dataset. Additionally,
we would like to thank Pierre Vandergheynst and his group at
EPFL for the inspiring discussions that led to this work.
R EFERENCES
[1] F. D. V. Fallani, J. Richiardi, M. Chavez, and S. Achard, “Graph
analysis of functional brain networks: Practical issues in translational
neuroscience,” Phil. Trans. R. Soc. B, vol. 369, no. 1653, p. 20130521,
2014.
[2] G. Camps-Valls, T. V. B. Marsheva, and D. Zhou, “Semi-supervised
graph-based hyperspectral image classification,” Geoscience and Remote
Sensing, IEEE Transactions on, vol. 45, no. 10, pp. 3044–3054, 2007.
[3] A. P. Dempster, “Covariance selection,” Biometrics, pp. 157–175, 1972.
[4] P. J. Bickel and E. Levina, “Covariance regularization by thresholding,”
The Annals of Statistics, vol. 36, no. 6, pp. 2577–2604, 2008. [Online].
Available: http://www.jstor.org/stable/25464728
[5] B. M. Lake and J. B. Tenenbaum, “Discovering structure by learning
sparse graph,” in Proceedings of the 33rd Annual Cognitive Science
Conference, 2010.
[6] X. Dong, D. Thanou, P. Frossard, and P. Vandergheynst, “Learning
laplacian matrix in smooth graph signal representations,” arXiv preprint
arXiv:1406.7842, 2014.
[7] V. Kalofolias, “How to learn a graph from smooth signals,” in Journal
of Machine Learning Research (JMLR), ser. Workshop and Conference
Proceedings, 2016.
[8] B. Girault, “Stationary Graph Signals using an Isometric Graph
Translation,” in Eusipco, Nice, France, Aug. 2015, pp. 1531–1535.
[Online]. Available: https://hal.inria.fr/hal-01155902
[9] N. Perraudin and P. Vandergheynst, “Stationary signal processing
on graphs,” CoRR, vol. abs/1601.02522, 2016. [Online]. Available:
http://arxiv.org/abs/1601.02522
[10] A. G. Marques, S. Segarra, G. Leus, and A. Ribeiro, “Stationary graph
processes and spectral estimation,” CoRR, vol. abs/1603.04667, 2016.
[Online]. Available: http://arxiv.org/abs/1603.04667
[11] A. Sandryhaila and J. Moura, “Discrete signal processing on graphs,”
Signal Processing, IEEE Transactions on, vol. 61, no. 7, pp. 1644–656,
2013.
[12] A. Sandryhaila and J. M. Moura, “Discrete signal processing on graphs:
Frequency analysis,” IEEE Transactions on Signal Processing, vol. 62,
no. 12, pp. 3042–3054, 2014.

16

[13] J. Mei and J. M. F. Moura, “Signal processing on graphs: Estimating
the structure of a graph,” in 2015 IEEE International Conference on
Acoustics, Speech and Signal Processing, ICASSP 2015, South Brisbane,
Queensland, Australia, April 19-24, 2015, 2015, pp. 5495–5499.
[Online]. Available: http://dx.doi.org/10.1109/ICASSP.2015.7179022
[14] D. Thanou, X. Dong, D. Kressner, and P. Frossard, “Learning
heat diffusion graphs,” CoRR, vol. abs/1611.01456, 2016. [Online].
Available: http://arxiv.org/abs/1611.01456
[15] H. Petric Maretic, D. Thanou, and P. Frossard, “Graph learning under
sparsity priors,” in International Conference on Acoustics, Speech and
Signal Processing (ICASSP), 2017.
[16] 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,” Signal Processing Magazine, IEEE, vol. 30, no. 3, pp. 83–98,
2013.
[17] A. Agaskar and Y. M. Lu, “A spectral graph uncertainty principle,”
Information Theory, IEEE Transactions on, vol. 59, no. 7, pp. 4338–
4356, 2013.
[18] N. Tremblay, “Réseaux et signal: des outils de traitement du signal pour
l’analyse des réseaux,” Ph.D. dissertation, Ecole Normale Supérieure de
Lyon, 2014.
[19] F. Chung, “Laplacians and the cheeger inequality for directed graphs,”
Annals of Combinatorics, vol. 9, no. 1, pp. 1–19, 2005. [Online].
Available: http://dx.doi.org/10.1007/s00026-005-0237-z
[20] F. R. Chung, Spectral Graph Theory. American Mathematical Soc.,
1997, vol. 92.
[21] Y. Shen, B. Baingana, and G. B. Giannakis, “Topology inference of
directed graphs using nonlinear structural vector autoregressive model,”
in International Conference on Acoustics, Speech and Signal Processing
(ICASSP), 2017.
[22] S. Sardellitti, S. Barbarossa, and P. Di Lorenzo, “Graph topology
inference based on transform learning,” in Signal and Information
Processing (GlobalSIP), 2016 IEEE Global Conference on. IEEE,
2016, pp. 356–360.
[23] T. Cai and W. Liu, “Adaptive thresholding for sparse covariance matrix
estimation,” Journal of the American Statistical Association, vol. 106,
no. 494, pp. 672–684, 2011.
[24] T. T. Cai, H. H. Zhou et al., “Optimal rates of convergence for sparse
covariance matrix estimation,” The Annals of Statistics, vol. 40, no. 5,
pp. 2389–2420, 2012.
[25] W. B. Wu and M. Pourahmadi, “Banding sample autocovariance matrices of stationary processes,” Statistica Sinica, pp. 1755–1768, 2009.
[26] H. Xiao, W. B. Wu et al., “Covariance matrix estimation for stationary
time series,” The Annals of Statistics, vol. 40, no. 1, pp. 466–493, 2012.
[27] J. Friedman, T. Hastie, and R. Tibshirani, “Sparse inverse covariance
estimation with the graphical lasso,” Biostatistics, vol. 9, no. 3, pp. 432–
441, 2008.
[28] A. J. Rothman, P. J. Bickel, E. Levina, J. Zhu et al., “Sparse permutation
invariant covariance estimation,” Electronic Journal of Statistics, vol. 2,
pp. 494–515, 2008.
[29] D. M. Witten, J. H. Friedman, and N. Simon, “New insights and faster
computations for the graphical lasso,” Journal of Computational and
Graphical Statistics, vol. 20, no. 4, pp. 892–900, 2011.
[30] R. Mazumder and T. Hastie, “The graphical lasso: New insights and
alternatives,” Electronic Journal of Statistics, vol. 6, p. 2125, 2012.
[31] K. M. Tan, D. Witten, and A. Shojaie, “The cluster graphical lasso
for improved estimation of gaussian graphical models,” Computational
Statistics & Data Analysis, vol. 85, pp. 23–36, 2015.
[32] S. Huang, J. Li, L. Sun, J. Liu, T. Wu, K. Chen, A. Fleisher, E. Reiman,
and J. Ye, “Learning brain connectivity of alzheimer’s disease from
neuroimaging data,” in Advances in Neural Information Processing
Systems, 2009, pp. 808–816.
[33] S. Yang, Q. Sun, S. Ji, P. Wonka, I. Davidson, and J. Ye, “Structural
graphical lasso for learning mouse brain connectivity,” in Proceedings
of the 21th ACM SIGKDD International Conference on Knowledge
Discovery and Data Mining, 2015, pp. 1385–1394.
[34] S. Sun, R. Huang, and Y. Gao, “Network-scale traffic modeling and
forecasting with graphical lasso and neural networks,” Journal of Transportation Engineering, vol. 138, no. 11, pp. 1358–1367, 2012.
[35] N. Wermuth, “Analogies between multiplicative models in contingency
tables and covariance selection,” Biometrics, pp. 95–108, 1976.
[36] E. Pavez and A. Ortega, “Generalized laplacian precision matrix estimation for graph signal processing,” in Acoustics, Speech and Signal
Processing (ICASSP), 2016 IEEE International Conference on. IEEE,
2016, pp. 6350–6354.

[37] H. E. Egilmez, E. Pavez, and A. Ortega, “Graph learning from data
under structural and laplacian constraints,” CoRR, vol. abs/1611.05181,
2016. [Online]. Available: http://arxiv.org/abs/1611.05181
[38] M. G. Rabbat, “Inferring sparse graphs from smooth signals with theoretical guarantees,” in International Conference on Acoustics, Speech
and Signal Processing (ICASSP), 2017.
[39] N. Perraudin, J. Paratte, D. Shuman, L. Martin, V. Kalofolias, P. Vandergheynst, and D. K. Hammond, “Gspbox: A toolbox for signal
processing on graphs,” ArXiv e-prints, 2014.
[40] S. P. Chepuri, S. Liu, G. Leus, and A. O. H. III, “Learning sparse graphs
under smoothness prior,” in International Conference on Acoustics,
Speech and Signal Processing (ICASSP), 2017.
[41] P. K. Shivaswamy and T. Jebara, “Laplacian spectrum learning,” in
Machine Learning and Knowledge Discovery in Databases. Springer,
2010, pp. 261–276.
[42] B. Pasdeloup, M. Rabbat, V. Gripon, D. Pastor, and G. Mercier, “Graph
reconstruction from the observation of diffused signals,” in Proceedings
of the 53rd Annual Allerton Conference on Communication, Control,
and Computing, 2015.
[43] S. Segarra, A. G. Marques, G. Mateos, and A. Ribeiro, “Network
topology inference from spectral templates,” CoRR, vol. abs/1608.03008,
2016. [Online]. Available: http://arxiv.org/abs/1608.03008
[44] ——, “Robust network topology inference,” in International Conference
on Acoustics, Speech and Signal Processing (ICASSP), 2017.
[45] ——, “Network topology identification from spectral templates,” CoRR, vol. abs/1604.02610, 2016. [Online]. Available:
http://arxiv.org/abs/1604.02610
[46] S. Shahrampour and V. M. Preciado, “Reconstruction of directed networks from consensus dynamics,” in American Control Conference
(ACC), 2013. IEEE, 2013, pp. 1685–1690.
[47] ——, “Topology identification of directed dynamical networks via power
spectral analysis,” IEEE Transactions on Automatic Control, vol. 60,
no. 8, pp. 2260–2265, 2015.
[48] M. Ipsen and A. S. Mikhailov, “Evolutionary reconstruction of networks,” Physical Review E, vol. 66, no. 4, p. 046109, 2002.
[49] R. Shafipour, A. G. Marques, G. Mateos, and A. Ribeiro, “Network
topology inference from non-stationary graph signals,” in International
Conference on Acoustics, Speech and Signal Processing (ICASSP), 2017.
[50] K. Pearson, “On lines and planes of closest fit to system of points in
space,” Philiosophical Magazine, vol. 2, pp. 559–572, 1901.
[51] T. W. Anderson, “Asymptotic theory for principal component analysis,”
Ann. Math. Statist., vol. 34, no. 1, pp. 122–148, 03 1963. [Online].
Available: http://dx.doi.org/10.1214/aoms/1177704248
[52] M. A. Girshick, “On the sampling theory of roots of determinantal
equations,” Ann. Math. Statist., vol. 10, no. 3, pp. 203–224, 09 1939.
[Online]. Available: http://dx.doi.org/10.1214/aoms/1177732180
[53] D. N. Lawley, “Tests of significance for the latent roots of covariance
and correlation matrices,” Biometrika, vol. 43, no. 1/2, pp. 128–136,
1956. [Online]. Available: http://www.jstor.org/stable/2333586
[54] E. Amaldi and V. Kann, “On the approximability of minimizing nonzero
variables or unsatisfied relations in linear systems,” Theoretical Computer Science, vol. 209, no. 1, pp. 237–260, 1998.
[55] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition
by basis pursuit,” SIAM Review, vol. 43, no. 1, pp. 129–159, 2001.
[56] R. Gribonval and M. Nielsen, “Sparse representations in unions of
bases,” Information Theory, IEEE Transactions on, vol. 49, no. 12, pp.
3320–3325, 2003.
[57] F. Rinaldi, “Mathematical programming methods for minimizing the
zero-norm over polyhedral sets,” Ph.D. dissertation, Sapienza, University
of Rome, 2009.
[58] M. Nekovee, “Worm epidemics in wireless ad hoc networks,” New
Journal of Physics, vol. 9, no. 6, p. 189, 2007.
[59] P. Erdős and A. Rényi, “On random graphs I,” Publ. Math. Debrecen,
vol. 6, pp. 290–297, 1959.
[60] M. Grant and S. Boyd, “CVX: MATLAB software for disciplined convex
programming, version 2.1,” http://cvxr.com/cvx, Mar. 2014.
[61] MATLAB, version 7.14.0 (R2012a).
Natick, Massachusetts: The
MathWorks Inc., 2012.
[62] B. Girault, “Signal processing on graphs-contributions to an emerging
field,” Ph.D. dissertation, Lyon, École normale supérieure, 2015.
[63] H. B. Mann and D. R. Whitney, “On a test of whether one of two
random variables is stochastically larger than the other,” Ann. Math.
Statist., vol. 18, no. 1, pp. 50–60, 03 1947. [Online]. Available:
http://dx.doi.org/10.1214/aoms/1177730491

