1

Graph Learning from Filtered Signals: Graph
System and Diffusion Kernel Identification

arXiv:1803.02553v1 [cs.LG] 7 Mar 2018

Hilmi E. Egilmez, Student Member, IEEE, Eduardo Pavez, Student Member, IEEE,
and Antonio Ortega, Fellow, IEEE

Abstract—This paper introduces a novel graph signal processing framework for building graph-based models from classes
of filtered signals. In our framework, graph-based modeling
is formulated as a graph system identification problem, where
the goal is to learn a weighted graph (a graph Laplacian
matrix) and a graph-based filter (a function of graph Laplacian
matrices). In order to solve the proposed problem, an algorithm
is developed to jointly identify a graph and a graph-based filter
(GBF) from multiple signal/data observations. Our algorithm is
valid under the assumption that GBFs are one-to-one functions.
The proposed approach can be applied to learn diffusion (heat)
kernels, which are popular in various fields for modeling diffusion
processes. In addition, for specific choices of graph-based filters,
the proposed problem reduces to a graph Laplacian estimation
problem. Our experimental results demonstrate that the proposed
algorithm outperforms the current state-of-the-art methods. We
also implement our framework on a real climate dataset for
modeling of temperature signals.
Index Terms—Graph learning, graph signal processing, graphbased filtering, graph system identification, diffusion kernels, heat
kernels.

G

I. I NTRODUCTION

RAPHS are fundamental mathematical structures used in
various fields to characterize data, signals and processes.
Particularly in signal processing, machine learning and statistics, structured modeling of signals/data by means of graphs
is essential in numerous problems including clustering [1]–
[3], regularized regression and denoising [4]–[6], where graphs
provide concise (sparse) representations for effective modeling
and analysis of signals/data [7]. Graph signal processing [6]
offers a new general paradigm for processing and analyzing
signals/data on graphs, referred as graph signals, by using
graph Laplacian matrices to extend basic signal processing
operations1 such as filtering [9], [10], transformation [11], [12]
and sampling [13] on graph signals. However, in practice,
datasets typically consist of an unstructured list of samples,
where the graph information (representing the structural relations between samples/features) is latent. For analysis, learning, processing and algorithmic purposes, it is often useful to
build graph-based models that provide a concise explanation
for datasets and also reduce the dimension of the problem
[14], [15] by exploiting the available prior knowledge and assumptions about the desired graph (e.g., structural information
Authors are with the Department of Electrical Engineering, University of
Southern California, Los Angeles, CA, 90089 USA. Contact author e-mail:
hegilmez@usc.edu. This paper was funded in part by NSF funding under
grant: CCF-1410009.
1 In [8], adjacency matrices are used to define basic processing operations.
In this work, we adopt the graph Laplacian based construction in [6].

xin

h(L)

xout

L graph Laplacian
h graph-based filter

Fig. 1. Input-output relation of a graph system defined by a graph Laplacian
(L) and a graph-based filter (h). In this work, we focus on joint identification
of L and h, given the observed data xout .

including connectivity and sparsity level) and observed data
(e.g., signal smoothness).
The focus of this paper is to build graph-based models from
signals/data, where the models of interest are defined by graph
systems consisting of a graph (graph Laplacian matrix L) and
a graph-based filter (GBF h) as depicted in Fig. 1. For graphbased modeling, graph systems provide a general abstraction,
in which graphs represent pairwise relations between the signal
samples and GBFs model smoothness of signals, defined on
a designated graph. In particular, graph systems (with specific
choices of GBFs) have been introduced to represent diffusion processes. Prominent examples include diffusion kernels
defined on graphs [16] and polynomials of graph Laplacian
matrices used to define localized diffusion operators on graphs
[6]. Models based on graph systems can be useful in a broad
number of applications including signal/data processing on
computer, social, sensor, energy, transportation and biological
networks, when signals are observed without the knowledge
of a graph and GBF associated with the underlying network.
For example, signals can be diffused on an unknown network
(L) with an unknown rate (β) as the smoothness parameter of
a GBF (hβ ).
In the literature, there has been a lot of recent interest
in the problem of learning graph-based models from data.
Several papers, including [17]–[19], consider the problem
of learning a graph, represented by a matrix such as graph
Laplacian, where the goal is to find the graph that provides
a best fit to the training data according to different criteria.
A few studies (see Section III for more details) consider
that the training data are the output of an unknown graph
system and propose methods under different assumptions on
GBFs. In [20], [21], the authors aim to learn a graph (i.e.,
a graph Laplacian or an adjacency matrix) from eigenvectors
of the sample covariance of training data, where the GBF is
an arbitrary function. In [22], the GBF is assumed to be a
polynomial function of a specific order, and the goal is to learn
a graph (an adjacency matrix) and coefficients of a polynomial
GBF. However, with such general assumptions on GBFs, it
is hard to provide guarantees about the performance of the
proposed methods. Specifically in [22], the solutions require

2

explicit sparsity conditions both on graphs and polynomials of
adjacency matrices, which are restrictive and may not hold in
some practical cases. Moreover, the approaches in [20], [21]
are based on the assumption that the sample covariance of the
observed data and the graph Laplacian have the same set of
eigenvectors. Thus, their performances strictly depend on the
accuracy of eigenvectors obtained from the sample covariance
of data, which is not a good estimator for the true covariance
of a model, especially when the number of data samples is
small [23], [24]. Indeed, better eigenvector estimates can be
obtained using inverse covariance estimation methods [24],
[25] or graph learning methods, such as those introduced in
our prior work [17], which minimize a regularized maximumlikelihood (ML) criterion.
In this paper, we propose an extension of our prior work
[17] to estimate the parameters of a graph system (a graph and
a GBF) from data. For this purpose, we formulate a graph
system identification (GSI) problem with a regularized ML
criterion, where our goal is to jointly identify a combinatorial
graph Laplacian (CGL) matrix L and a GBF h from multiple
signal/data observations under the main assumption that the
GBF h is a one-to-one function. While this assumption is
not as general as the assumptions on GBFs made in existing
studies [20]–[22], several useful GBFs have this property, as
shown in Table I. The key novelty of our work is to propose
new methods that rely on the one-to-one GBF property to learn
both graph and GBF with stronger optimality guarantees as
compared to the approaches in [20]–[22]. In Section IV-D,
we will show that this assumption is one of the sufficient
conditions required for graph system identifiability, and it
also allows our algorithm to perform a prefiltering step that
significantly improves the estimation accuracy. Although our
algorithm can be extended for any one-to-one GBFs, we focus
on methods to estimate the parametric GBFs listed in Table
I, which are one-to-one functions that depend on a single
parameter β and have useful applications discussed in Section
IV-C. The first three GBF types in the table define basic
scaling and shifting operations, and the exponential decay and
β-hop localized GBFs provide diffusion-based models. Since
all these GBFs yield larger filter responses (hβ ) in lower graph
frequencies (λ) as illustrated in Fig. 2, they can be used for
modeling different classes of smooth graph signals satisfying
that most of the signal energy is concentrated in the low graph
frequencies.
In order to solve the GSI problem, we propose an alternating
optimization algorithm that first optimizes the graph by fixing
the GBF and then designs the GBF by fixing the graph.
Basically, the proposed algorithm involves three main steps,
which are repeated until convergence is achieved:
• The graph-based filter (GBF) identification step designs
a GBF h for current estimate of graph Laplacian L so
that its inverse (i.e., h−1 ) will be used for the prefiltering
step. Note that h has to be a one-to-one function for its
inverse h−1 to exist.
−1
• The prefiltering step filters the observed signals using h
to compute the covariance matrix that will be used in the
graph Laplacian estimation step. We show that this step
significantly improves the accuracy of graph estimation.

The graph Laplacian estimation step learns a graph
from the covariance of prefiltered signals obtained in the
previous step by using the combinatorial graph Laplacian
(CGL) estimation algorithm introduced in our prior work
[17]. Although the present paper focuses on graphs associated with CGL matrices, our solution can be easily
extended to other types of graph Laplacian matrices, e.g.,
generalized graph Laplacians [26].
In order to accommodate the GBFs in Table I in our algorithm, we propose specific methods (for the GBF identification
step) to find the filter parameter β fully characterizing the
selected GBF. Our proposed algorithm guarantees optimal
identification of β and L in an `1 -regularized ML sense
for frequency shifting, variance shifting and β-hop localized
GBFs. However, for frequency scaling and exponential decay
GBFs, our algorithm cannot find the optimal β in general,
but it guarantees that the estimated L is optimal up to a
constant factor. In practice, the type of GBF and its parameter
β can also be selected based on the prior knowledge available
about the problem and the application. In this case, different
GBFs and their parameters can be tested until the estimated
graph and GBF pair achieves the desired performance (e.g.,
in terms of mean square error, likelihood or sparsity) where
the parameter β serves as a regularization parameter, which
can be used for tuning smoothness of the signal, for example.
Moreover, our algorithm can be extended to support GBFs
beyond the filters in Table I (including one-to-one functions
with more than one parameter) by developing specific methods
for the GBF identification step. As long as a specified GBF
(h) has an inverse function (h−1 ), the proposed prefiltering
and graph Laplacian estimation steps can be directly utilized
to learn graphs from signals/data.
To the best of our knowledge, this is the first work that (i)
formulates the graph-based modeling problem as identification
of graph systems with the types of GBFs in Table I under
a regularized ML criterion and (ii) proposes a prefilteringbased algorithm to jointly identify a graph and a GBF. Existing
related studies (see Section III) consider optimization of different criteria (see Table III), and do not use prefiltering, which
can be shown to be optimal in some cases (see Section V).
The proposed approach can significantly improve the accuracy
of graph learning from filtered signals, as compared to the
existing methods that estimate a graph directly from observed
signals without prefiltering. Particularly, if observed signals are
diffused/filtered on an unknown network/graph to be learned,
applying a graph learning algorithm (e.g., CGL algorithm in
[17]) on diffused signals potentially results in a dense graph
due to diffusion, even when the underlying graph is sparse. On
the other hand, our proposed algorithm reverses the effect of
diffusion via prefiltering (h−1 ), whose corresponding GBF (h)
is jointly estimated with the graph (L). Thus, the latent graph
can be learned more accurately. Note that the diffusion (heat)
kernels [16] used in a number of applications [4], [5], [27] are
special cases of graph systems with exponential decay GBFs,
and thus our proposed algorithm can be applied to learn their
parameters.
The rest of the paper is organized as follows. Section
II presents the notation and some concepts used throughout
•

3

TABLE I
GBF S WITH PARAMETER β AND CORRESPONDING INVERSE FUNCTIONS
†
WHERE λ DENOTES THE PSEUDOINVERSE OF A SCALAR λ, THAT IS ,
λ† = 1/λ FOR λ 6= 0 AND λ† = 0 FOR λ = 0.
Filter Name

(

Frequency
scaling

hβ (λ)
1
βλ

0

λ>0
λ=0

Frequency
shifting

(λ + β)†

Variance
shifting

λ† + β

Exponential
decay

h−1
β (s)
1
βs

s>0
s=0

0
1
s

exp(−βλ)
(
1
λ>0
λβ
0
λ=0

β-hop
localized

(

(

−β

Range of β
β ∈ R, β > 0
β ∈ R, β ≥ 0

(s − β)†

β ∈ R, β ≥ 0

−log(s)/β

β ∈ R, β > 0

1
1
s

0

β

s>0
s=0

β ∈ N, β ≥ 1

TABLE II
L IST OF S YMBOLS AND T HEIR M EANING
Symbols

Meaning

L | L, Lc

graph Laplacian or CGL matrix | set of CGLs

λi , λi (L)

i-th eigenvalue of L in ascending order

h, hβ | H
n|k

GBF function | set of GBFs

number of vertices | number of data samples

I|W|D

identity matrix | adjacency matrix | degree matrix

Θ−1 | Θ†

inverse of Θ | pseudo-inverse of Θ

1|0

ΘT | θ T

det(Θ) | |Θ|
(Θ)ij

column vector of ones | column vector of zeros
transpose of Θ | transpose of θ

determinant of Θ | pseudo-determinant of Θ
entry of Θ at i-th row and j-th column

(θ)i

i-th entry of θ

≥ (≤)

element-wise greater (less) than or equal to operator
Θ is a positive semidefinite matrix

Θ0

Tr | logdet(Θ)

trace operator | natural logarithm of det(Θ)

p(x)

probability density function of random vector x

x ∼ N(0, Σ)

zero-mean multivariate Gaussian with covariance Σ

kθk22 , kΘk2F

sum of squared values of all elements

kθk1 , kΘk1

sum of absolute values of all elements (`1 -norm)

the paper. Section III discusses the prior related work. In
Section IV, the GSI problem and its variations are formulated.
Additionally, the graph system identifiability conditions are
introduced. In Section V, we derive optimality conditions and
develop methods to solve the GSI problem. The experimental
results are presented in Section VI, and Section VII draws
some conclusions.
II. N OTATION AND P RELIMINARIES
Throughout the paper, lowercase normal (e.g., a and θ),
lowercase bold (e.g., a and θ) and uppercase bold (e.g., A and
Θ) letters denote scalars, vectors and matrices, respectively.
Unless otherwise stated, calligraphic capital letters (e.g., E and
S) represent sets. Notation is summarized in Table II.
The graph-based models considered in this paper are defined
based on undirected, simple weighted graphs with nonnegative
edge weights and no self-loops. Let G = (V, E, fw ) be a simple
weighted graph with n vertices in the set V = {v1 , . . . , vn },

where the edge set E = { e | fw (e) 6= 0, ∀ e ∈ Pu } is a subset
of Pu , the set of all unordered pairs of distinct vertices, and
fw ((vi , vj )) ≥ 0 for i 6= j is a real-valued edge weight function.
The adjacency matrix of G is an n × n symmetric matrix, W,
such that (W)ij = (W)ji = fw ((vi , vj )) for (vi , vj ) ∈ Pu .
The degree matrix
G is an n × n diagonal matrix, D, with
Pof
n
entries (D)ii = j=1 (W)ij and (D)ij = 0 for i 6= j. The
combinatorial graph Laplacian (CGL) of G is defined as L =
D − W. Alternatively, the set of CGL matrices can also be
written as
Lc = {L | L  0, (L)ij ≤ 0 for i 6= j, L1 = 0} .

(1)

By construction, CGLs are symmetric and positive semidefinite matrices, so each of them has a complete set of orthogonal eigenvectors u1 , u2 , . . . , un whose associated eigenvalues λ1 ≤ λ2 ≤ · · · ≤ λn are nonnegative real numbers2 . In
addition, the CGL of a connected graph always has a zero
eigenvalue (i.e., λ1 = 0 with
√ multiplicity one) whose associated
eigenvector is u1 = (1/ n)1.
In graph signal processing, the eigenpairs of a CGL matrix,
(λ1 , u1 ), (λ2 , u2 ), . . . , (λn , un ), provide a Fourier-like spectral interpretation for signals defined on graphs, where the
graph frequency spectrum is defined by eigenvalues of the
CGL, which are called graph frequencies, and eigenvectors
of the CGL u1 , u2 , . . . , un are the harmonics associated with
the graph frequencies. Based on the eigenvectors of a CGL
matrix L, the graph Fourier transform (GFT) is defined as the
orthogonal matrix U (i.e., satisfying UT U = I) obtained by
eigendecomposition of L = UΛUT , where Λ denotes the diagonal matrix consisting of eigenvalues λ1 , λ2 , . . . , λn (graph
frequencies).
For a given signal x = [x1 x2 · · · xn ] defined on a graph
G with n vertices, where xi is attached to vi (i-th vertex), the
b = UT x. Naturally, the variation of the GFT
GFT of x is x
basis vectors gradually increase with the graph frequencies,
and the GBF basis vectors corresponding to low frequencies
are relatively smooth.
As a specific example, the GFT basis
√
vector u1 = (1/ n)1, associated with lowest graph frequency
(λ1 = 0), is the smoothest among other GFT basis vectors.
To quantify smoothness of a graph signal x, a common
variation measure used in graph signal processing is the graph
Laplacian quadratic form, xT Lx, which can be written in
terms of edge weights of G as
X
2
xT Lx =
(W)ij (xi − xj )
(2)
T

(i,j)∈I

where (W)ij = −(L)ij , and I = {(i, j) | (vi , vj ) ∈ E} is the set
of index pairs of vertices associated with the edge set E. A
smaller Laplacian quadratic term (xT Lx) indicates a smoother
graph signal (x). A more general notion of smoothness can
be obtained by defining graph-based filters (GBFs), i.e., matrix functions h(L) = Uh(Λ)UT , where U is the GFT and
(h(Λ))ii = h((Λ)ii ) = h(λi ) ∀i. Note that a GBF h serves as
a scalar function of λ that maps graph frequencies λ1 , . . . , λn

2 Without loss of generality, we assume that the eigenvalues are ordered as
λ1 ≤ λ2 ≤ · · · ≤ λn throughout the paper.

0

2

4

λ

6

8

10

(a) Frequency scaling GBF
Fig. 2.

β = 0.25
β = 0.50
β = 0.75

0

2

4

λ

6

8

10

(b) Frequency shifting GBF

1.4
1.2
1
0.8
0.6
0.4
0.2
0

β = 0.25
β = 0.50
β = 0.75

0

2

4

λ

6

8

(c) Variance shifting GBF

10

1.4
1.2
1
0.8
0.6
0.4
0.2
0

β = 0.25
β = 0.50
β = 0.75
hβ (λ)

1.4
1.2
1
0.8
0.6
0.4
0.2
0

hβ (λ)

β=1
β=2
β=3

hβ (λ)

1.4
1.2
1
0.8
0.6
0.4
0.2
0

hβ (λ)

hβ (λ)

4

0

2

4

λ

6

8

10

(d) Exponential decay GBF

1.4
1.2
1
0.8
0.6
0.4
0.2
0

β=1
β=2
β=3

0

2

4

λ

6

8

10

(e) β-hop localized GBF

The illustration of graph-based filters (GBFs) in Table I with different β parameters as a function of graph frequency hβ (λ).

to filter responses3 h(λ1 ), . . . , h(λn ). Thus, the measure in (2)
can be generalized using a GBF h as
n
X
xT h(L)x = xT Uh(Λ)UT x =
h(λi )b
x2i ,
(3)
i=1

where x
bi = uTi x for i = 1, . . . , n denote GFT coefficients.
Based on a graph and a GBF, we formally define a graph
system as follows.

Definition 1 (Graph System). A graph system is defined by a
simple graph G and a GBF h such that the input-output relation
of the system is xout = h(L)xin , where L ∈ Lc is a CGL matrix
associated with G, and xin is the input signal vector.

In this paper, the graph system parameter h is selected from
a set of GBFs H determined based on the GBF types in Table
I. As a specific example, the exponential decay GBFs lead to
the set H = {h | h(λ) = exp(−βλ), β ∈ R and β > 0} and
the operator h(L) = exp(−βL), which is also known as the
diffusion (heat) kernel defined on a graph [16].
Definition 2 (Diffusion Kernels on Graphs). The diffusion
kernel over graph G is the matrix exponential of L, that is
t

βL
t ∈ N,
(4)
exp(−βL) = lim I −
t→∞
t
where L denotes a graph Laplacian associated with G and the
parameter β is a real number called diffusion bandwidth.
Appendix A presents a derivation of diffusion kernels (i.e.,
exponential GBFs) based on a basic class of random diffusion
processes defined by CGLs.
III. R ELATED W ORK
Table III summarizes prior related studies by comparing
against our present work in terms of (i) target variables,
(ii) underlying assumptions and (iii) optimization criteria. In
our previous work on graph learning from data [17], we
have proposed algorithms for estimating models based on
three different types of graph Laplacians (including CGLs).
The present paper extends our prior work by introducing
more general graph-based models based on graph systems,
defined by a CGL and a GBF, and proposing an algorithm to
learn parameters of a graph system from multiple signal/data
observations. In the literature, there are several papers on graph

3 Filter responses corresponding to the eigenvalues with multiplicity more
than one have the same value. Formally, if λi = λi+1 = · · · = λi+c−1
for some i > 1 and multiplicity c > 1, then h(λi ) = h(λi+1 ) = · · · =
h(λi+c−1 ).

learning from signals/data4 . Methods to estimate CGLs from
smooth signals are proposed in [18] and [19], but, in contrast
with our work, GBFs are not considered explicitly. Instead, the
graph estimation problem is formulated as minimization of a
regularized graph Laplacian quadratic form (i.e., a smoothness
metric for graph signals). In [20], [21], the authors focus on
learning graph shift/diffusion operators (e.g., adjacency and
graph Laplacian matrices) from a complete set of eigenvectors
that has to be computed from observed data. More specifically,
Segarra et al. [20] solve a sparse recovery problem by minimizing the `1 -norm of the target variable (i.e., minimizing
kΘk1 where Θ is the target variable) to infer the topology
of a graph, and Pasdeloup et al. [21] estimate an adjacency
matrix by minimizing its trace, Tr(Θ), as well as its `1 -norm,
kΘk1 . In fact, the problems in [20] and [21] are equivalent
if the target matrix is constrained to be a CGL, since the
problems of minimizing kLk1 and Tr(L) over the set of
CGL matrices (i.e., L ∈ Lc ) lead to the same solution. Note
that both of these methods [20], [21] only use the eigenvectors of the sample covariance, while its eigenvalues (which
also carry graph information implicitly) are not exploited in
graph estimation. Although our approach also requires a set
eigenvectors to be computed, they are not directly used to
estimate a graph. Instead, the computed eigenvectors and their
associated eigenvalues are used in the prefiltering step, then a
graph is estimated from the covariance of prefiltered signals
by minimizing a regularized ML criterion. The optimality of
our prefiltering-based approach is discussed in Section V.
Thanou et al. [28] address the estimation of a graph and
a sparse input (i.e., localized sources) from a set of observed
signals, and they propose a dictionary-based method, where a
graph estimation problem is solved to construct a dictionary
consisting of multiple diffusion kernels and the resulting
dictionary is used for identifying the localized sources in
the diffusion process. Noting that our paper focuses on the
GSI problem without locality assumptions on diffused/filtered
signals (see Table III), when no locality assumptions are
imposed and a single diffusion kernel is used in the dictionary,
the problem in [28] reduces to the graph estimation problem
in [18], formulated as minimization of a regularized graph
Laplacian quadratic form. In contrast with the work in [28],
our algorithm iteratively solves a different graph learning problem with a regularized ML criterion (i.e., the CGL estimation
problem in [17]) and also performs prefiltering on observed
signals. Since our algorithm often results in a more accurate
4 We

refer to [17] for a discussion of methods beyond CGL estimation.

5

TABLE III
A N OVERVIEW OF THE RELATED WORK ON LEARNING GRAPHS FROM SMOOTH / DIFFUSED SIGNALS .
References

Target Variable(s)

Graph

graph Laplacian

None

None

None

regularized Laplacian quadratic form

Segarra et al. [20]

graph shift operator

None

None

None

`1 -norm of graph shift operator

Pasdeloup et al. [21]

adjacency matrix

None

None

None

trace and `1 -norm of adjacency

Thanou et al. [28]

multiple heat kernels
source locations

None

Heat
kernel

Localized

regularized least squares

Mei and Moura [22]

polynomials of adjacency

Sparse

Polynomial

None

regularized least squares

Egilmez et al. [17]

graph Laplacian

None

None

None

regularized maximum likelihood

This work

graph Laplacian
GBF (Table I)

None

One-to-one
function

None

regularized maximum likelihood

IV. P ROBLEM F ORMULATION : G RAPH S YSTEM
I DENTIFICATION
Our formulation is based on the following general assumption on a GBF h, which holds for the GBFs in Table I.
Assumption 1. We assume that a graph-based filter h(λ) is a
nonnegative and one-to-one function of λ.
A. Filtered Signal Model
We formulate the GSI problem in a probabilistic setting
by assuming that the observed (filtered) signals have been
sampled from a zero mean n-variate Gaussian distribution6
practice, powers of adjacency/Laplacian matrices lead to dense matrices.
zero mean assumption is made to simplify the notation. Our model
can be trivially extended to a multivariate Gaussian distributions with nonzero
mean.
6 The

Optimization Criterion

Signal

Dong et al. [18]
Kalofolias [19]

graph estimation as compared to the method in [18] (see
Section VI for the details), it can be extended to construct
potentially better dictionaries as alternatives to the ones in
[28] used for identifying localized sources. Such an extension
is out of the scope of the present paper. In [22], Mei and
Moura address the estimation of polynomials of adjacency
matrices by solving a regularized least-squares problem. As
their counterparts, polynomials of graph Laplacians can be
used to approximate the GBFs in Table I as well as many other
types of filters such as bandlimited GBFs. Since polynomial
filters provide more degrees of freedom in designing GBFs,
they can be considered as more general compared to our
GBFs of interest. However, polynomials are not one-to-one
functions in general, so our proposed prefiltering step cannot
be applied. The algorithm in [22] provides some optimality
guarantees in a mean-square sense, but this is under a restrictive set of assumptions, which require the polynomials
of adjacency matrices to be sparse5 . Our proposed algorithm
provides stronger optimality guarantees in a regularized ML
sense without imposing any assumptions on the sparsity of
graphs (see Table III), as long as GBFs are one-to-one. The
identification of graph systems with polynomial GBFs will be
studied as part of our future work.

5 In

Assumptions
Filter

x ∼ N(0, h(L)),




1 T
†
p(x|h(L)) =
exp − x h(L) x , (5)
n/2
1/2
2
(2π) |h(L)|
1

with the covariance Σ = h(L) defined based on a CGL matrix
L and a GBF h.
For modeling smooth graph signals, it is reasonable to
choose h(λ) to be a monotonically decreasing function7 satisfying h(λ1 ) ≥ h(λ2 ) ≥ · · · ≥ h(λn−1 ) ≥ h(λn ) > 0,
where the graph frequencies (eigenvalues of L) are ordered as
0 = λ1 ≤ λ2 ≤ · · · ≤ λn . Thus, the corresponding covariance
matrix h(L) represents graph signals whose energy is larger in
lower graph frequencies. Note that the nonnegativity condition
in Assumption 1, that is, h(λi ) ≥ 0 for i = 1, . . . , n, ensures
that the covariance h(L) is a positive semidefinite matrix.
From the graph system perspective (see Fig. 1 and Definition 1), the above probabilistic model corresponds to the
case when the input is the n-variate white Gaussian noise
xin ∼ N(0, Σin = I). Then, the covariance of the output vector
T
xout = h(L)xin becomes Σout = h(L)Σin h(L) = h(L)2 ,
which can be translated into our model in (5) by simply
1/2
transforming Σout as Σ = Σout
= h(L), where the mapping
between h(L)2 and h(L) is one-to-one (i.e., there is no loss
of information after transformation), since eigenvalues of h(L)
are nonnegative.
B. Problem Formulation
Our goal is to estimate the graph system parameters L and
h based on k random samples, xi for i = 1, . . . , k, obtained
from the signal model in (5) by maximizing the likelihood of
h(L), that is


k
k
Y
Y
k
1
p(xi |h(L)) = c |h(L)|− 2
exp − xi T h(L)† xi
(6)
2
i=1
i=1
− kn

with a constant c = (2π) 2 . The maximization of (6) can be
equivalently stated as minimizing its negative log-likelihood,
7 All of the GBFs in Table I are monotonically decreasing functions on the
interval λ ∈ (0, ∞). The exponential decay and frequency shifting GBFs
further satisfies h(λ1 ) ≥ h(λ2 ) at the zero frequency (i.e., 0 = λ1 ≤ λ2 ),
while for the other GBFs, we have h(λ2 ) ≥ h(λ1 ).

6

which leads to the following problem for estimating the graph
system parameters L and h,
)
( k
 k
1X  T
†
†
b
b
(L, h) = argmin
Tr xi h(L) xi − log|h(L) |
2 i=1
2
L,h
n 

o
†
†
= argmin Tr h(L) S − log|h(L) |
L,h

(7)
where S denotes the sample covariance calculated using k
samples, xi for i = 1, 2, . . . , k. In our problem formulation,
we additionally introduce a sparsity promoting weighted `1 regularization, kL Hk1 for robust graph estimation, where
H is symmetric regularization matrix and denotes elementwise multiplication. Thus, the proposed GSI problem is formulated as follows:
†
†
minimize Tr(hβ (L) S)−log|hβ (L) |+kL Hk1
L0,β
(8)
subject to L 1 = 0, (L)ij ≤ 0 i 6= j
where β is the (unknown) parameter for a given type of GBF
hβ from Table I, and H denotes the regularization matrix. The
constraints ensure that L is a CGL matrix8 . In this work, we
particularly choose H = α(2I−11T ) so that the regularization
term in (8) reduces to the standard `1 -regularization αkLk1
with parameter α.
C. Special Cases of Graph System Identification Problem and
Applications of Graph-based Filters
Graph Learning. The problem in (8) can be reduced to the
CGL problem proposed in our previous work [17] by choosing
the filter hβ to be the frequency scaling or β-hop localized
GBFs with β = 1, so that we have
(
1/λ λ > 0
†
hβ (λ) = λ =
that is hβ (L) = L† .
(9)
0
λ=0
Since hβ (L) = L, we can rewrite the objective function in (8)
as
Tr(LS) − log|L|+kL Hk1
(10)
†

which is the objective function of graph learning problems
originally introduced in [17].
Graph Learning from Noisy Signals. The variance shifting
filter corresponds to a noisy signal model with variance β =
σ 2 . Specifically, assuming that the noisy signal is modeled as
b = x + e, where x ∼ N(0, L† ) denotes the original signal,
x
and e ∼ N(0, σ 2 I) is the additive white Gaussian noise vector
b
independent of x with variance σ 2 , then the noisy signal x
b = L† +σ 2 I). The same model is
b ∼ N(0, Σ
is distributed as x
obtained by using a variance shifting filter with β = σ 2 so that
hβ (λ) = λ† + β and hβ (L) = L† + βI = L† + σ 2 I.

(11)

With this type of GBFs, our GSI problem in (8) can be solved
to learn graphs from noisy signals by identifying L and the
noise parameter β = σ 2 . This problem can also be viewed
as an extension of our formulations in [17] derived based on
signal models that are assumed to be noise-free (i.e., β = 0), so
8 The formulation can be trivially extended for different types of graph
Laplacians (e.g., generalized graph Laplacian) discussed in [17].

it can be solved to improve the performance of graph learning
from data in the presence of noise [29].
Graph Learning from Frequency Shifted Signals. For the
shifted frequency filter with parameter β, we have
hβ (λ) = (λ + β)† and hβ (L) = (L + βI)† .

(12)

By substituting hβ (L) with (L+βI) , the problem in (8) can
be written as
e
e
e Hk
minimize Tr(LS)−log|
L|+k
L
1
e
L0,β≥0
(13)
e = L+βI, L1 = 0, (L)ij ≤ 0 i 6= j
subject to L
which is the shifted combinatorial Laplacian (SCGL) estimation problem that is originally proposed in [30].
Diffusion (Heat) Kernel Learning. To learn diffusion kernels,
formally defined in Definition 2, the proposed GSI problem
in (8) can be modified by choosing hβ (λ) as the exponential
decay filter, so we get
†

hβ (L) = Uhβ (Λ)UT = Uexp(−βΛ)UT = exp(−βL). (14)
The resulting problem can be solved to learn diffusion kernels
from signals/data, which are popular in many applications [4],
[5], [27].
Graph Learning from β-hop Localized Signals. To learn
graphs from β-hop localized signals, where β is a positive
integer, the GBF can be selected as h(λ) = (λ† )β so that the
corresponding inverse covariance (i.e., precision) matrix in (8)
is h(L)† = Lβ , which is generally not a graph Laplacian matrix
due to the exponent β. However, it defines diffusion operators
(on a graph associated with L) that can be used as alternatives
to heat kernels, and the corresponding signal model leads to βhop localized/diffused signals, in which each sample depends
on samples located within the neighborhood at most β-hops
away.
D. Graph System Identifiability
In statistical learning theory, a probabilistic model is identifiable if it is possible to learn the model parameters exactly from
infinite number of data samples. This requires different model
parameters to generate different probability distributions [31],
which is formally stated in the following definition.
Definition 3 (Model Identifiability [31]). Let M = {p(x|Θ) :
Θ ∈ PΘ } be the family of probabilistic models defined by the
parameter space PΘ . M is identifiable if the mapping from
PΘ to M is one-to-one, that is, if Θ1 6= Θ2 then p(x|Θ1 ) 6=
p(x|Θ2 ) for arbitrary Θ1 , Θ2 ∈ PΘ .
Based on the above definition, we introduce three different
notions of identifiability for graph systems characterizing the
probabilistic model discussed in Section IV-A.

Definition 4 (Notions of Identifiability). Let the parameter set
of the signal model in (5) be PΘ = {(L, h) | L ∈ L, h ∈ H}
such that a class of graph systems is defined by a specific
choice of sets L and H, denoting admissible graphs and filters
in the class. For any such class of graph systems, we define
the following types of identifiability:
• A class is filter identifiable if h1 6= h2 implies h1 (L) 6=
h2 (L) for arbitrary L ∈ L and h1 , h2 ∈ H.

7

•
•

A class is graph identifiable if L1 6= L2 implies h(L1 ) 6=
h(L2 ) for arbitrary h ∈ H and L1 , L2 ∈ L.
A class is jointly identifiable if any two pairs satisfying
(L1 , h1 ) 6= (L2 , h2 ) lead to h1 (L1 ) 6= h2 (L2 ) for
arbitrary L1 , L2 ∈ L and h1 , h2 ∈ H.

The following proposition categorizes the classes of systems
defined based on parametric GBFs in Table I in terms of their
identifiability.
Proposition 1. Let Lc and H define a class of graph systems.
All classes corresponding to the parametric GBFs in Table I
are both filter and graph identifiable. In addition, the classes
with β-hop localized, frequency shifting or variance shifting
GBFs in Table I are jointly identifiable.
Proof: The proof follows from Definition 4 by checking types
of identifiability for each GBF with different β parameters.
For filter identifiability the proof is straightforward, because
two different parameters β1 6= β2 imply hβ1 (L) 6= hβ2 (L).
Since GBFs of interest are one-to-one functions, corresponding
classes are also graph identifiable by Proposition 2. The classes
with frequency scaling and exponential decay GBFs are not
jointly identifiable, since we have hβ1 (L1 ) = hβ2 (L2 ) when
L1 = (β2 /β1 )L2 . However, no such construction leading
to hβ1 (L1 ) = hβ2 (L2 ) exists for the other GBF types as
shown in (24)–(30), thus the corresponding classes are jointly
identifiable.

For a general choice of H, the following proposition states
the sufficient condition for graph identifiability.
Proposition 2. A class of graph systems defined by L and H
is graph identifiable if L is a set of CGLs and H is a set of
one-to-one functions.
Proof: Assuming that all h ∈ H are one-to-one functions, we
need to show that for L1 6= L2 in L, we have h(L1 ) 6= h(L2 ).
Specifically, if L1 = UΛ1 UT and L2 = UΛ2 UT , which have
the same GFT, Uh(Λ1 )UT 6= Uh(Λ2 )UT is satisfied, because
h is one-to-one. The proof is obvious when GFTs of L1 and
L2 are different.

To derive necessary conditions for graph identifiability, additional set of assumptions are needed. For example, if graphs
(i.e., graph Laplacians in L) are assumed to be sparse, the oneto-one requirement on H may not be necessary. Particularly in
[32], identifiability of graphs from their eigenspaces (i.e., from
GFTs as in [20]) are investigated, and necessary conditions are
introduced for a specific type of identifiability, called diagonal
identifiability. To the best of our knowledge, there are no other
analysis on the identifiability of graphs in the literature. As
part of our future work, we will focus on characterizations
of necessary and sufficient conditions for joint identifiability
under more general choices of H.
V. P ROPOSED S OLUTION AND O PTIMALITY C ONDITIONS
Algorithm 1 is proposed to solve the GSI problem in (8)
for a given sample covariance S and a selected type of GBF
hβ . After obtaining U and Λs via eigendecomposition of S
and initialization of the parameter β (see lines 1 and 2), the

Algorithm 1 Graph System Identification (GSI)
Input: Sample covariance S, graph-based filter type hβ
Output: Graph Laplacian L and β filter parameter
1: Obtain U and Λs via eigendecomposition S = UΛs UT
b For variance or frequency shifting
2: Initialize parameter β:
GBFs, apply the initialization methods in Section V-C. For
the other GBF types, apply random initialization.
3: repeat
4:
Prefilter the sample covariance S:
β

5:

6:

7:
8:

(15)

−1
†
†
T
Spf = (h−1
b (S)) = U(h b (Λs )) U

β

Estimate L from prefiltered data (Spf ):
b ← Run the CGL algorithm in [17] to solve (23)
L
Update filter parameter βb or skip if βb is optimal:
βb ← Apply GBF-specific update in Section V-C

until convergence has been achieved
b and βb
return Graph system parameters L

algorithm performs three main steps to find the optimal graph
Laplacian L and the filter parameter β:
1) Prefiltering step applies an inverse filtering on the sample covariance S to reverse the effect of filter h in
the graph system. Without the prefiltering step, it may
be impossible to effectively recover L from S, since
Σ† = h(L)† is generally not a graph Laplacian. The
proposed prefiltering allows us to effectively estimate the
original eigenvalues of L (i.e., Λλ ) from the prefiltered
covariance Spf in (15).
2) Graph Laplacian estimation step uses the CGL estimab from the
tion algorithm introduced in [17] to learn L
prefiltered covariance Spf .
3) Filter parameter selection step finds the best matching
β for the graph system. Depending on the type of GBF,
we propose different methods for parameter selection.
In the rest of this section, we derive the optimality conditions and discuss optimal prefiltering, graph Laplacian estimation and filter parameter selection in Sections V-A, V-B and
V-C, respectively.
A. Optimal Prefiltering
Let CGL L and GBF h be the parameters of a graph system
so that the covariance matrix of the model in (5) is Σ = h(L).
There is a GFT matrix U jointly diagonalizing L and Σ such
that L = UΛλ UT and Σ = UΛσ UT . Under the ideal case that
S = Σ is obtained from an asymptotically large number of
samples (k), by change of variables, the objective function in
(7) becomes
†

†

J (U, h(Λλ )) = Tr(Uh(Λλ ) Λσ UT ) − log|Uh(Λλ ) UT |
(16)
which is simplified using properties of Tr(·) and |·| as
†

†

J (h(Λλ )) = Tr(h(Λλ ) Λσ ) − log|h(Λλ ) |

(17)

where the GBF h and the diagonal matrix of graph frequencies
Λλ are unknown, and the diagonal matrix Λσ is known from

8

data. By letting φi = h(λi ) = (h(Λλ ) )ii and σi2 = (Λσ )ii for
i = 1, 2, . . . , n, we can write (17) as
n
X

J (φ) =
φi σi2 − log(φi ) ,
(18)
†

†

i=1

where φ = [φ1 φ2 · · · φn ] . In minimization of the convex
function (18), the optimal solution satisfies the following
necessary and sufficient conditions [33] obtained by taking
the derivative of (18) with respect to φi and equating to zero,
1
∂J (φ)
=
− σi2 = 0.
(19)
∂φi
φi
By change of variables, the optimality conditions can be stated
as
h(λi ) = (h(Λλ ))ii = (Λσ )ii = σi2 ∀i.
(20)
T

Initialization for Variance/Frequency Shifting Filters. For
both variance and frequency shifting GBFs, the optimal β is
found by calculating σ12 = u1T Σu1 where u1 is the eigenvector associated with the zero eigenvalue of the Laplacian
((Λλ )11 = λ1 = 0). Specifically, by using the optimality condition in (21),
•

hβ (λ1 ) = σ12 = u1T Σu1 = λ1 + β,
†

since λ1 = λ1 = 0, the optimal β satisfies
β = u1T Σu1 = σ12 .
•

h−1 (Σ) = Uh−1 (Λσ )UT = UΛλ UT = L.

(22)

This condition shows that we can find the optimal Laplacian
L via inverse filtering (inverse prefiltering) h−1 (Σ). However,
in practice, we can only have access to a sample estimate of
Σ (i.e., S) obtained from a limited number of data samples
(k), which is not a good estimator especially when k is small
[24]. Thus, computing h−1 (S) generally does not lead to a
CGL matrix. In order to address this problem, the proposed
Algorithm 1 first estimates the prefiltered sample covariance
Spf as in (15), then employs the CGL estimation algorithm [17]
to find the best CGL from Spf by minimizing the criterion in
(10).
B. Optimal Graph Laplacian Estimation
For a GBF h (or hβ ) satisfying the optimal prefiltering
condition in (21), the GSI problem in (8) can be rewritten
as the following graph learning problem,
minimize Tr(LSpf ) − log|L|+kL Hk1
L0
(23)
subject to L 1 = 0, (L)ij ≤ 0 i 6= j
where Spf = (h−1 (S))† is obtained by prefiltering the empirical
covariance S. This problem is discussed in detail in our previous work [17] where we have derived the optimality conditions
for (23) and developed the CGL estimation algorithm to solve
it.
C. Filter Parameter Selection
Based on the optimality condition in (21), we propose different methods to identify the parameter β for GBFs in Table I.
Specifically, optimal parameter initializations for variance and
frequency shifting GBFs are derived, and a line search method
is proposed for β-hop localized GBF. Since graph systems with
exponential decay and frequency scaling GBFs are not jointly
identifiable (as discussed in Section IV-D), we cannot identify
β optimally. Yet, we show that graph Laplacian matrices can
be identified up to a constant factor, which depends on β.

(25)

if hβ (λ) is a frequency shifting filter, we obtain
hβ (λ1 ) = σ12 = u1T Σu1 = 1/(λ1 + β) = 1/β,

(26)

so the optimal β satisfies

(21)

where h−1 is the inverse function of h. By using the matrix
notation, we can state (21) more compactly as

(24)

†

Based on Assumption 1, we can also write (20) as
h−1 (h(λi )) = λi = h−1 (σi2 ) ∀i,

if hβ (λ) is a variance shifting filter, we get

β = 1/(u1T Σu1 ) = 1/σ12 .

(27)

Since the optimal β can be directly estimated from the sample
covariance S as calculated from Σ in (25) and (27), Algorithm
1 uses the optimized initial β (in line 2) for prefiltering, and
then estimates L, so that the filter parameter update (line 6) is
skipped for graph systems with frequency and variance shifting
filters.
Exponential Decay and Frequency Scaling Filters. Suppose
b that is
b is obtained by inverse prefiltering of Σ with β,
that L
−1
b
formally, L = h b (Σ).
β

If hβ (λ) is a frequency scaling filter, we have
1
β
2
2
h−1
b (σi ) : λ1 = σ1 = 0, b λi = b 2 i = 2, . . . , n.
β
β
βσi
(28)
• If hβ (λ) is an exponential decay filter, we have
log(σi2 )
log(exp(−βλi )) β
h−1
(σi2 ) = −
=−
= λi . (29)
b
β
b
β
βb
βb
b
Based on (28) and (29), for any selected β, the resulting
b where L and β denote
b = (β/β)L
matrix is a CGL satisfying L
the original graph system parameters. Since the inverse of a
b
different GBF (i.e., h−1
b with β 6= β) leads to a CGL for any
β
β, Algorithm 1 can only find the optimal CGL matrix L up
b In practice, the parameter βb can be
to a constant factor β/β.
tuned so that the desired normalization (scaling factor) for L
is achieved.
•

β-hop Localized Filter. For estimation of the optimal hop
count β in Algorithm 1, inverse prefiltering with βb gives,
 1/βb
b
1
β/β
−1
2
2
i = 2, . . . , n.
h b (σi ) : λ1 = σ1 = 0, λi =
β
σi2
(30)
b
Since this requires the graph learning step to estimate Lβ/β ,
which is not a graph Laplacian in general, Algorithm 1 cannot
b In
guarantee optimal graph system identification unless β = β.
order to find the optimal β, we perform a line search for given
range of integers optimizing the following:
b − S||F .
βb = argmin||h b(L)
(31)
b
β∈N

β

9

VI. R ESULTS
A. Graph Learning from Diffusion Signals/Data
We evaluate the performance of our proposed graph system
identification method (GSI or Algorithm 1) by benchmarking
against the current state-of-the-art approaches proposed for
learning graph from smooth signals (GLS) [18], [19] as well
as the graph topology inference (GTI) in [20]. The proposed
GSI is also compared against the CGL estimation algorithm,
CGL [17] (i.e., using the CGL estimation algorithm without
prefiltering), and the inverse prefiltering (IPF) approach, which
estimates a graph Laplacian matrix by inverting the prefiltered
b = h−1 (S) = S† . For this
covariance, Spf in (15), so that L
pf
β
purpose, we generate several artificial datasets based on the
signal model in (5), defined by a graph Laplacian (L) and a
GBF (hβ ) where the dataset entries are generated by random
sampling from the distribution N(0, hβ (L)). Then, the generated data is used in the proposed and benchmark algorithms
to recover the corresponding graph Laplacian L. We repeat
our experiments for different L and hβ where graphs are
constructed by using three different graph connectivity models:
•
•
•

(n)
Grid graph, Ggrid
, consisting n vertices attached to their
four nearest neighbors (except the vertices at boundaries).
(n,p)
Random Erdos-Renyi graph, GER
, with n vertices attached to other vertices with probability p = 0.2.
Random modular graph (also known as stochastic block
(n,p1 ,p2 )
model), GM
with n vertices and four modules (subgraphs) where the vertex attachment probabilities across
modules and within modules are p1 = 0.1 and p2 = 0.2,
respectively.

Then, the edge weights of a graph are randomly selected from
the uniform distribution U(0.1, 3), on the interval [0.1, 3]. For
each L and hβ pair, we perform Monte-Carlo simulations to
test average performance of proposed and benchmark methods
with varying number of data samples (k) and fixed number of
vertices (n = 36)9 . To measure the estimation performance,
we use the following two metrics:
∗
b
b L∗ ) = kL − L kF
(32)
RE(L,
kL∗ kF
which is the relative error between the ground truth graph (L∗ )
b and
and estimated graph parameters (L),
2 tp
b L∗ ) =
(33)
FS(L,
2 tp + fn + fp
is the F-score metric (commonly used metric to evaluate binary
classification performance) calculated based on true-positive
(tp), false-positive (fp) and false-negative (fn) detection of
b with respect to the ground truth
graph edges in estimated L
∗
edges in L . F-score takes values between 0 and 1, where the
value 1 means perfect classification.
In our experiments, for the proposed GSI, the regularization
parameter α in (23) is selected from the following set:
p
{0} ∪ {0.75r (smax log(n)/k) | r = 1, 2, 3, . . . , 14}, (34)
9 Methods are evaluated on small graphs (with n = 36 vertices), since GTI
[20] is implemented using CVX [34] and do not currently have an efficient
and scalable implementation. The proposed and the other benchmark methods
are more efficient and can support larger graphs.

where smax = maxi6=j |(S)ij | is the maximum off-diagonal
enp
try of S in absolute value, and the scaling term log(n)/k
is used for adjusting the regularization according to k and
n as suggested in [24], [35]. Monte-Carlo simulations are
performed for each proposed/baseline method, by successively
solving the associated problem with different regularization
parameters to find the best regularization that minimizes
RE. The corresponding graph estimate is also used to calculate FS. For all baseline methods [18]–[20], the required
parameters are selected by fine tuning. Since CGL [17],
GLS [18], [19] and GTI [20] approaches generally result in
severely biased solutions with respect to the ground truth
L∗ (based on our observations from the experiments), RE
values are calculated after normalizing the estimated solution
b = (Tr(L∗ )/Tr(L))L. Note that, this normalization also
L as L
resolves the ambiguity in identification of graph systems with
exponential decay and frequency scaling filters up to a scale
factor (discussed in Section V-C).
Figures 3 and 4 depict the performance of different methods
applied for estimating graphs from signals modeled based
on exponential decay filters (diffusion kernels) and β-hop
localized filters. As shown in Figs. 3 and 4, the proposed GSI
significantly outperforms all baseline methods, including the
state-of-the-art GLS [18], [19] and GTI [20], in terms of average RE and FS metrics. The performance difference between
GSI and CGL [17] demonstrates the impact of the prefiltering
step, which substantially improves the graph learning accuracy.
Similarly, the performance gap between GSI and IPF shows the
advantage of Algorithm 1 compared to the direct prefiltering of
input covariance (S) as in (22), where GSI provide better graph
estimation especially when number of data samples (i.e., k/n)
is small. Besides, Figs. 5 and 6 illustrate two examples from
the experiments with grid graphs for the case of k/n = 30,
where the proposed GSI constructs graphs that are the most
similar to the ground truth (L∗ ).
B. Graph Learning from Variance/Frequency Shifted Signals
In this subsection, we compare the CGL estimation performance of GSI, CGL [17] and SCGL [30] methods from signals
modeled based on variance and frequency shifting GBFs. As
discussed in Section IV-C, the covariance matrices for signals
modeled based on these GBFs with parameter β are
†
• Σ = L +βI for variance shifting,
†
• Σ = (L+βI) for frequency shifting.
where L denotes the associated combinatorial Laplacian.
In our experiments, we construct 10 random Erdos-Renyi
(n,p)
graphs (GER
) with n = 36 vertices and p = 0.2, then generate
Σ for each GBF by varying β between 0 and 1. To evaluate
the effect of β only, we use actual covariance matrices instead
of sample covariances as input to the algorithms. So, GSI, CGL
and SCGL estimate a graph Laplacian L from Σ. The average
RE results are presented in Tables IV and V for various β.
Table IV shows that the proposed GSI significantly outperforms CGL for β > 0, and the average RE difference increases
as β gets larger. This is because the variance shifting GBF
leads to the noisy signal model with the covariance in (11)
where β represents the variance of the noise (σ 2 ), and the

10

0.4

CGL
GLS
GLS
GTI

0.35
Average RE

0.3

1

GSI
IPF
[17]
[18]
[19]
[20]

0.9
0.8
Average FS

0.45

0.25
0.2
0.15

0.6
0.5

0.1

CGL
GLS
GLS
GTI

0.4

0.05
0

0.7

5

10

30

50

100

250

500

0.3

1000

5

10

30

50

k/n

100

GSI
IPF
[17]
[18]
[19]
[20]

250

500

1000

k/n

Fig. 3. Average RE and FS results for graph estimation from signals modeled based on exponential decay GBFs tested with β = {0.5, 0.75} on 10 different
grid, Erdos-Renyi and modular graphs (30 graphs in total). The proposed GSI outperforms all baseline methods in terms of both RE and FS.

CGL
GLS
GLS
GTI

0.4

Average RE

0.35
0.3

1

GSI
IPF
[17]
[18]
[19]
[20]

0.9
0.8
Average FS

0.5
0.45

0.25
0.2

0.6
0.5

0.15
0.1

CGL
GLS
GLS
GTI

0.4

0.05
0

0.7

0.3
5

10

30

50

100

250

500

1000

k/n

5

10

30

50

100

GSI
IPF
[17]
[18]
[19]
[20]

250

500

1000

k/n

Fig. 4. Average RE and FS results for graph estimation from signals modeled based on β-hop localized GBFs tested with β = {2, 3} on 10 different grid,
Erdos-Renyi and modular graphs (30 graphs in total). The proposed GSI outperforms all baseline methods in terms of both RE and FS.
TABLE IV
AVERAGE R ELATIVE E RRORS FOR VARIANCE S HIFTING GBF
Filter parameter (β)
Method

0

0.1

CGL [17]

2×10−4

0.60

0.3

0.5

0.7

0.9

0.79

0.85

0.88

0.89

2 × 10−4

GSI

TABLE V
AVERAGE R ELATIVE E RRORS FOR F REQUENCY S HIFTING GBF
Filter parameter (β)
Method
SCGL [30]
GSI

0

0.1

0.5

0.9

0.2354

6.7×10−4

6.6×10−4

6.2×10−4

1.6 ×

10−4

prefiltering step allows GSI to perfectly estimate the parameter
β from Σ by using (25) so that the covariance is prefiltered
as in (15) based on the optimal β. The prefiltering step can
also be considered as a denoising operation (reversing the
effect of variance shifting GBFs) on the signal covariance
before the graph estimation step, while CGL works with noisy
(i.e., shifted) covariances, which diminish the CGL estimation
performance. For β = 0 (i.e., Σ is noise-free), the problem
(8) reduces to the CGL estimation problem in [17], so both
GSI and CGL lead to the same average RE.

For the frequency shifting GBFs with β > 0, GSI performs
slightly better than SCGL, since SCGL is implemented using
a general purpose solver CVX [34], which generally produces
less accurate solutions compared to our algorithm. Moreover,
our algorithm is approximately 90 times faster than SCGL on
average, and significantly outperforms SCGL for β = 0, since
SCGL method is developed for shifted covariance matrices
(i.e., L + βI) where β needs to be strictly positive.
C. Illustrative Results on Temperature Data
In this experiment, we apply our proposed method on a real
(climate) dataset10 consisting of air temperature measurements
[36]. We specifically use the average daily temperature measurements collected from 45 states in the US over 16 years
(2000-2015), so that in total there are k = 5844 samples
for each of the n = 45 states. Fig. 7 shows samples of
average temperature signals, which are spatially smooth across
different states. Also, the Rocky Mountains region has lower
average temperature values as compared to the other regions
in the western US11 .
10 NCEP Reanalysis data provided by the NOAA/OAR/ESRL PSD, Boulder,
Colorado, USA, from their website at http://www.esrl.noaa.gov/psd/
11 The Rocky Mountains cross through the states of Idaho, Montana,
Wyoming, Utah, Colorado and New Mexico. For example, in Fig. 7b, the
areas with temperature values between 0 and 10 degrees Celsius (colored in
green) correspond to the Rocky Mountains region approximately.

11

0

(a) The ground truth CGL

2.9746

0

(L∗ )

1.8932

0

(b) GLS [18]
(RE,FS) = (0.4019,0.6172)

2.6617

0

(c) GTI [20]
(RE,FS) = (0.1201,0.7229)

3.2115

(d) Proposed GSI
(RE,FS) = (0.0340,0.9917)

Fig. 5. A sample illustration of graph estimation results (for k/n = 30) from signals modeled using the exponential decay GBF with β = 0.75 and L∗ is
derived from the grid graph in (a). The edge weights are color coded where darker colors indicate larger weights. The proposed GSI leads to the graph that
is the most similar to the ground truth.

0

2.9307

(a) The ground truth CGL (L∗ )

0

2.9115

(b) GLS [18]
(RE,FS) = (0.4045,0.6890)

0

2.0173

(c) GTI [20]
(RE,FS) = (0.2118,0.6030)

0

3.0194

(d) Proposed GSI
(RE,FS) = (0.0274,0.9833)

Fig. 6. A sample illustration of graph estimation results (for k/n = 30) from signals modeled using the β-hop localized GBF with β = 2 and L∗ is derived
from the grid graph in (a). The edge weights are color coded where darker colors indicate larger weights. The proposed GSI leads to the graph that is the
most similar to the ground truth.

The goal of this experiment is to learn graphs (with 45 vertices) associated with exponential decay and β-hop localized
filters from temperature data and investigate the effect of filter
type and parameter β on the resulting graphs, representing
the similarity of temperature conditions between the 45 states.
For modeling temperature signals, a diffusion kernel is a
good candidate, because it is a fundamental solution of the
heat equation12 , which describes the distribution of heat in a
physical environment [37].
Figure 8 illustrates the graphs estimated using the GSI
method without `1 -regularization (i.e., H in (8) is set to the
zero matrix). As shown in the figure, larger edge weights are
assigned between vertices (i.e., states) that are closer to each
other in general, since temperature values are mostly similar
between states within close proximity. However, the distance
between states is obviously not the only factor effecting the
similarity of temperature values. For example, in Figs. 8b–8f,
the weights are considerably small between the states in the
Rocky Mountains region and their neighboring states in the
east (e.g., Nebraska and Kansas) due to the large differences
in altitude. Note also that different choices of GBFs can lead
to substantially different similarity graphs. Especially for the
β-hop localized GBF with β = 1 (corresponding to the CGL
12 This is the reason that diffusion kernels are also known as heat kernels
in the literature.

method), the resulting graph is significantly different than the
results in Figs. 8b and 8c, since the 1-hop localized filter
does not lead to a diffusion model. The graphs associated
with exponential decay GBF leads to sparser graphs, better
revealing the structure of the signal, compared to the graphs
in Figs. 8a–8c. The structure of the graphs in Figs. 8d–8f
are similar for different β because of the relation in (29) for
the exponential decay filter. For example, increasing β from
0.25 to 0.5 approximately halves edge weights, as shown in
Figs. 8d and 8e. Besides, the distribution of edge weights for
β-hop localized GBFs becomes more similar to the ones in
Figs. 8d–8f as β gets larger.
VII. C ONCLUSION
In this paper, we have introduced a novel graph-based
modeling framework that (i) formulates the modeling problem as the graph system identification from signals/data and
(ii) proposes an alternating optimization algorithm iteratively
solving for a graph and a GBF. At each iteration of the
algorithm, a prefiltering operation, defined by the estimated
GBF, is applied on the observed signals, and then a graph
is estimated from prefiltered signals. The experimental results
have demonstrated that the proposed algorithm outperforms
the existing methods on modeling smooth signals and learning
diffusion-based models [18]–[20] in terms of graph estimation
accuracy.

12

(a) 45th day (Winter)
Fig. 7.
states.

(b) 135th day (Spring)

(c) 225th day (Summer)

(d) 315th day (Autumn)

Average air temperatures (in degree Celsius) for (a) 45th, (b) 135th, (c) 225th and (d) 315th days over 16 years (2000-2015). Black dots denote 45

0

83.3921

0

4.7463

(a) β-hop localized GBF with β = 1

(b) β-hop localized GBF with β = 5

0

0

9.5083

(d) Exponential decay with β = 0.25

4.7722

(e) Exponential decay with β = 0.5

0

0.80833

(c) β-hop localized GBF with β = 10

0

0.81572

(f) Exponential decay with β = 3

Fig. 8. Graphs learned from temperature data using the GSI method with exponential decay and β-hop localized GBFs for fixed β parameters where
no regularization is applied (i.e., H is set to the zero matrix). The edge weights are color coded so that darker colors represent larger weights (i.e., more
similarities). The graphs associated with exponential decay GBF leads to sparser structures compared to the graphs corresponding to β-hop localized GBFs.

The proposed framework supports various types of GBFs
(in Table I) including the diffusion (heat) kernels as special
cases which are widely used in many applications [4], [5], [27].
Our future work focuses on the extensions of our algorithm
for joint identification of graphs and polynomial filters (i.e.,
estimation of polynomials of graph Laplacians), which can
provide more degrees of freedom in designing filters than
the GBFs in Table I. Also, data-oriented applications of the
proposed modeling framework is considered as part of our
future work.
A PPENDIX A
D ERIVATION OF D IFFUSION K ERNELS
Suppose that the random process is obtained by diffusion
of the initial random vector x(0) whose entries are independent, zero-mean random variables, denoted as (x(0))i for
i = 1, 2, . . . , n, with variance σ 2 , the random diffusion over a
graph G is formulated recursively as
x(t+1) = x(t)−rLx(t) = (I−rL)x(t) t ∈ {0, 1, . . . } (35)

where t represents the time index, L is the graph Laplacian
of G and the real number r ∈ (0, 1) is the diffusion rate of
the process. On i-th vertex of the graph G, we can write the
random process in (35) as
X
(x(t + 1))i = (x(t))i + r
(W)ij ((x(t))j −(x(t))i ) (36)
j∈Ni

for all i where W is the adjacency matrix of the graph G, and
Ni is the set of vertex indices neighboring i-th vertex (vi ).
More compactly, (35) can be written as
t

x(t) = G(t)x(0) where G(t) = (I − rL) .

(37)

The covariance of x(t) is
(Σ(t))ij = E [(G(t)x(0))i (G(t)x(0))j ]
" n
! n
!#
X
X
=E
(G(t))ik (x(0))k
(G(t))jl (x(0))l ,
k=1

l=1

(38)
which simplifies by using independence of x(0) as
n
X
(Σ(t))ij = σ 2
((G(t))ik (G(t))kj ) = σ 2 (G(t)2 )ij . (39)
k=1

13

Therefore, the covariance matrix leads to
2

2

2t

Σ(t) = G(t) = σ (I − rL) .

(40)

To adjust the time resolution of the process, we replace t with
t/∆t and r with r∆t in G(t) so that
t/∆t

rL
.
(41)
G(t) = σ I −
1/∆t
For arbitrarily small ∆t, G(t) converges to a matrix exponential function of L,

2t/∆t
rL
2
e
Σ(t) = σ lim I −
= σ 2 exp(−2rtL) (42)
∆t→0
1/∆t
which is equivalent to exponential decay filtering operation in
Table I with β(t) = −2rt. For arbitrarily large t (i.e., when
the diffusion reaches steady-state), the covariance is
n
X
σ2 T
e
11 , (43)
lim Σ(t)
= lim
exp(−2rtλi )ui ui T =
t→∞
t→∞
n
i=1
whose all entries have the same value, since we have λ1 = 0
for CGLs. Thus, the statistical properties eventually become
spatially flat across all vertices as intuitively expected for a
diffusion process. Yet, for nonsingular graph Laplacians (e.g.,
generalized graph Laplacians [17], [26]), the process behaves
differently so that the signal energy vanishes, since the above
limit converges to the n×n zero matrix as t gets larger.
R EFERENCES
[1] A. Y. Ng, M. I. Jordan, and Y. Weiss, “On spectral clustering: Analysis
and an algorithm,” in Proceedings of the 14th International Conference
on Neural Information Processing Systems: Natural and Synthetic, ser.
NIPS’01. Cambridge, MA, USA: MIT Press, 2001, pp. 849–856.
[2] M. Soltanolkotabi, E. Elhamifar, and E. J. Candès, “Robust subspace
clustering,” Ann. Statist., vol. 42, no. 2, pp. 669–699, 04 2014.
[3] U. Luxburg, “A tutorial on spectral clustering,” Statistics and Computing,
vol. 17, no. 4, pp. 395–416, Dec. 2007.
[4] A. J. Smola and I. R. Kondor, “Kernels and regularization on graphs.”
in Proceedings of the Annual Conference on Computational Learning
Theory, 2003.
[5] A. D. Szlam, M. Maggioni, and R. R. Coifman, “Regularization on
graphs with function-adapted diffusion processes,” J. Mach. Learn. Res.,
vol. 9, pp. 1711–1739, Jun. 2008.
[6] D. Shuman, S. Narang, P. Frossard, A. Ortega, and P. Vandergheynst,
“The emerging field of signal processing on graphs: Extending highdimensional data analysis to networks and other irregular domains,”
IEEE Signal Processing Magazine, vol. 30, no. 3, pp. 83–98, 2013.
[7] P. Buhlmann and S. van de Geer, Statistics for High-Dimensional Data:
Methods, Theory and Applications. Springer Publishing Co., Inc., 2011.
[8] A. Sandryhaila and J. M. F. Moura, “Discrete signal processing on
graphs: Frequency analysis,” IEEE Transactions on Signal Processing,
vol. 62, no. 12, pp. 3042–3054, June 2014.
[9] P. Milanfar, “A tour of modern image filtering: New insights and methods, both practical and theoretical,” IEEE Signal Processing Magazine,
vol. 30, no. 1, pp. 106–128, Jan 2013.
[10] H. E. Egilmez and A. Ortega, “Spectral anomaly detection using graphbased filtering for wireless sensor networks,” in 2014 IEEE International
Conference on Acoustics, Speech and Signal Processing (ICASSP), May
2014, pp. 1085–1089.
[11] H. E. Egilmez, Y. H. Chao, A. Ortega, B. Lee, and S. Yea, “GBST:
Separable transforms based on line graphs for predictive video coding,”
in 2016 IEEE International Conference on Image Processing (ICIP),
Sept 2016, pp. 2375–2379.
[12] S. K. Narang and A. Ortega, “Perfect reconstruction two-channel wavelet
filter banks for graph structured data,” IEEE Transactions on Signal
Processing, vol. 60, no. 6, pp. 2786–2799, June 2012.
[13] A. Anis, A. Gadde, and A. Ortega, “Efficient sampling set selection
for bandlimited graph signals using graph spectral proxies,” IEEE
Transactions on Signal Processing, vol. 64, no. 14, pp. 3775–3789, July
2016.

[14] T. Hastie, R. Tibshirani, and J. Friedman, The elements of statistical
learning: data mining, inference and prediction, 2nd ed. Springer,
2008.
[15] Y. S. Abu-Mostafa, M. Magdon-Ismail, and H.-T. Lin, Learning From
Data. AMLBook, 2012.
[16] R. I. Kondor and J. D. Lafferty, “Diffusion kernels on graphs and other
discrete input spaces,” in Proceedings of the Nineteenth International
Conference on Machine Learning, ser. ICML ’02. San Francisco, CA,
USA: Morgan Kaufmann Publishers Inc., 2002, pp. 315–322.
[17] H. E. Egilmez, E. Pavez, and A. Ortega, “Graph learning from data
under Laplacian and structural constraints,” IEEE Journal of Selected
Topics in Signal Processing, vol. 11, no. 6, pp. 825–841, Sept 2017.
[18] X. Dong, D. Thanou, P. Frossard, and P. Vandergheynst, “Learning
Laplacian matrix in smooth graph signal representations,” IEEE Transactions on Signal Processing, vol. 64, no. 23, pp. 6160–6173, Dec 2016.
[19] V. Kalofolias, “How to learn a graph from smooth signals,” in Proceedings of the 19th International Conference on Artificial Intelligence and
Statistics (AISTATS), May 2016, pp. 920–929.
[20] S. Segarra, A. G. Marques, G. Mateos, and A. Ribeiro, “Network
topology inference from spectral templates,” IEEE Transactions on
Signal and Information Processing over Networks, vol. 3, no. 3, pp.
467–483, Sept 2017.
[21] B. Pasdeloup, V. Gripon, G. Mercier, D. Pastor, and M. G. Rabbat,
“Characterization and inference of graph diffusion processes from
observations of stationary signals,” IEEE Transactions on Signal and
Information Processing over Networks, vol. PP, no. 99, pp. 1–1, 2017.
[22] J. Mei and J. M. F. Moura, “Signal processing on graphs: Causal modeling of unstructured data,” IEEE Transactions on Signal Processing,
vol. 65, no. 8, pp. 2077–2092, April 2017.
[23] I. M. Johnstone and A. Y. Lu, “On consistency and sparsity for principal
components analysis in high dimensions,” Journal of the American
Statistical Association, vol. 104, no. 486, pp. 682–693, 2009.
[24] P. Ravikumar, M. Wainwright, B. Yu, and G. Raskutti, “High dimensional covariance estimation by minimizing l1-penalized log-determinant
divergence,” Electronic Journal of Statistics (EJS), vol. 5, pp. 935–980,
2011.
[25] J. Friedman, T. Hastie, and R. Tibshirani, “Sparse inverse covariance
estimation with the graphical lasso,” Biostatistics, vol. 9, no. 3, pp. 432–
441, Jul. 2008.
[26] T. Bıyıkoglu, J. Leydold, and P. F. Stadler, “Laplacian eigenvectors of
graphs,” Lecture notes in mathematics, vol. 1915, 2007.
[27] J. Lafferty and G. Lebanon, “Diffusion kernels on statistical manifolds,”
J. Mach. Learn. Res., vol. 6, pp. 129–163, Dec. 2005.
[28] D. Thanou, X. Dong, D. Kressner, and P. Frossard, “Learning heat diffusion graphs,” IEEE Transactions on Signal and Information Processing
over Networks, vol. 3, no. 3, pp. 484–499, Sept 2017.
[29] M. G. Rabbat, “Inferring sparse graphs from smooth signals with
theoretical guarantees,” in 2017 IEEE International Conference on
Acoustics, Speech and Signal Processing (ICASSP), March 2017, pp.
6533–6537.
[30] B. M. Lake and J. B. Tenenbaum, “Discovering structure by learning
sparse graph,” in Proceedings of the 33rd Annual Cognitive Science
Conference, 2010, pp. 778–783.
[31] E. L. Lehmann and G. Casella, Theory of Point Estimation, 2nd ed.
New York, NY, USA: Springer-Verlag, 1998.
[32] Y. D. Castro, T. Espinasse, and P. Rochet, “Reconstructing undirected
graphs from eigenspaces,” Journal of Machine Learning Research,
vol. 18, no. 51, pp. 1–24, 2017.
[33] D. P. Bertsekas, Nonlinear Programming.
Belmont, MA: Athena
Scientific, 1999.
[34] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex
programming, version 2.1,” http://cvxr.com/cvx, Mar. 2014.
[35] S. Zhou, P. Rutimann, M. Xu, and P. Buhlmann, “High-dimensional
covariance estimation based on Gaussian graphical models,” J. Mach.
Learn. Res., vol. 12, pp. 2975–3026, Nov. 2011.
[36] E. Kalnay, M. Kanamitsu, R. Kistler, W. Collins, D. Deaven, L. Gandin,
M. Iredell, S. Saha, G. White, J. Woollen, Y. Zhu, A. Leetmaa,
R. Reynolds, M. Chelliah, W. Ebisuzaki, W. Higgins, J. Janowiak, K. C.
Mo, C. Ropelewski, J. Wang, R. Jenne, and D. Joseph, “The ncep/ncar
40-year reanalysis project,” Bulletin of the American Meteorological
Society, vol. 77, no. 3, pp. 437–471, 1996.
[37] L. C. Evans, Partial differential equations. Providence, R.I.: American
Mathematical Society, 2010.

