1

Kernel Regression for Signals over Graphs

arXiv:1706.02191v3 [cs.IT] 30 Jan 2018

Arun Venkitaraman, Saikat Chatterjee, Peter Händel
School of Electrical Engineering and ACCESS Linnaeus Center
KTH Royal Institute of Technology, SE-100 44 Stockholm, Sweden
arunv@kth.se, sach@kth.se, ph@kth.se

Abstract—We propose kernel regression for signals over
graphs. The optimal regression coefficients are learnt using the
assumption that the target vector is a smooth signal over an
underlying graph. The condition is imposed using a graphLaplacian based regularization. We discuss how the proposed
kernel regression exhibits a smoothing effect, simultaneously
achieving noise-reduction and graph-smoothness. We further
extend the kernel regression to simultaneously learn the underlying graph and the regression coefficients. We validate our
theory by application to various synthesized and real-world graph
signals. Our experiments show that kernel regression over graphs
outperforms conventional regression, particularly for small sized
training data and under noisy training. We also observe that
kernel regression reveals the structure of the underlying graph
even with a small number of training samples.
Index Terms—Linear model, regression, kernels, machine
learning, graph signal processing, graph-Laplacian.

EDICS−NEG-SPGR, NEG-ADLE, MLR-GRKN.
I. I NTRODUCTION
Kernel regression constitutes one of the fundamental building blocks in supervised and semi-supervised learning strategies, be it for simple regression tasks or in more advanced
settings such as support vector machines [1], Gaussian processes [2], deep neural networks [3], [4] and extreme learning
machines [5], [6]. Kernel regression is a data-driven approach
which works by projecting the input vector into a higher
dimensional space using kernel functions and performs a
linear regression in that space [7], [8]. The projection is
achieved typically by using a nonlinear transformation. The
target output y (for both scalar and vector cases) for a given
test input x is then expressible as a linear combination of the
kernels between the test input and the training inputs {xi }:
X
y=
ρ i k(x, xi ),
i

where k(·, ·) denotes the kernel function. ρ i s denote the
coefficient vectors learnt from the training input-output pairs
{(xi , yi )}. Given sufficiently large number of samples of
reliable training data, kernel regression leads to accurate
target predictions. However, in many cases, such as medical
experiments or epidemic studies, it might not be feasible
to collect large amounts of clean data. Insufficient or noisy
training data often leads to deterioration in the prediction
performance of kernel regression, just as in the case of any
data-driven approach. In such cases, incorporating additional
The authors would like to acknowledge the support received from the
Swedish Research Council.

structure or regularization during the training phase can boost
performance. In this article, we propose kernel regression for
signals over graphs, where it is assumed that target signals
are smooth over graphs: we employ the graph signal structure
to enhance prediction performance of kernel regression under
noisy and scarce data. Experiments with real-world data such
as temperature sensor network, electroencephalogram(EEG),
atmospheric tracer diffusion, and functional MRI data demonstrate the potential of our approach.
Recently, graph signal processing has emerged as a promising area for incorporating graph-structural information in the
analysis of vector signals [9], [10]. This takes particular
relevance in contemporary applications that deal with data over
networks or graphs. Several traditional signal processing methods have been extended to graphs [9]–[12] including many
conventional spectral analysis concepts such as the windowed
Fourier transforms, filterbanks, multiresolution analysis, and
wavelets, [13]–[26]. Spectral clustering approaches based on
graph signal filtering have been proposed recently [27], [28].
A number of graph-aware sampling schemes have also been
developed [29]–[36]. Principal component analysis (PCA)
techniques were proposed for graph signals by Shahid et
al. for recovery of low-rank matrices [37], [38] and were
shown to perform on par with nuclear-norm based approaches.
Parametric dictionary learning approaches have demonstrated
that graph specific information can help to outperform standard
approaches like the K-SVD at small sample sizes [39]–[43].
Statistical analysis of graph signals particularly with respect
to stationarity have also been investigated [44]–[48]. Berger
et al. [49] recently proposed a primal-dual based convex
otpimization approach for recovery of graph signals from
a subset of nodes by minimizing the total variation while
controlling the global empirical error. A number of approaches
have also been proposed for learning the graph-Laplacian
directly from data using different metrics [50]–[54].
Graph-based kernels have been investigated for labeling/coloring of graph nodes and graph clustering [55]–[62].
They deal with binary valued signals or data over the nodes.
Kernel-based reconstruction strategies for graph signals were
proposed recently by Romero et al. in the framework of
reproducing kernel Hilbert spaces [63], [64]. Using the notion
of joint space-time graphs, Romero et al. have also proposed
kernel based reconstruction of graph signals and an extension
of the Kalman filter for kernel-based learning [65], [66].
Recently, Shen et al. used kernels in structural equation models
for identification of network topologies from graph signals
[67]. Though these works employ kernel-based strategies they

2

differ from our proposed approach by formulation. The prior
works define kernels across nodes of the same graph, while
considering a subset of the nodes to be the input. We consider
that the target to be predicted is a vector signal lying over
a graph, and the input signal is not necessarily associated
with the underlying graph. Given a new input, we provide
a prediction output which is a smooth graph signal.
A. Our contributions
We assume that the output of kernel regression and the
target vector are smooth over a graph with a specified graphLaplacian matrix. By noting that the input observations enter
the regression model only in the form of inner-products, we
develop kernel regression over graphs (KRG) using a dual
representation or kernel trick. We also discuss how our approach may be applied to cases when the graph is not specified
apriori. In this case, we jointly learn the graph-Laplacian
and the regression coefficients by formulating an alternating
minimization algorithm. In other words, we jointly perform the
prediction for unseen inputs and learning of underlying graph
topology. Experiments with both synthesized and real-world
data confirm our hypothesis that a significant gain in prediction
performance vis-a-vis conventional kernel regression can be
obtained by regularizing the target vector to be a smooth
signal over an underlying graph. Our experiments also include
realistic noisy training scenarios of missing samples and large
perturbations.
B. Signal processing over graphs
We next briefly review the concepts from graph signal
processing. Let G = (V, E, A) denote a graph with M
nodes indexed by the vertex set V = {1, · · · , M }. Let E
and A denote the edge set containing pairs of nodes, and
the weighted adjacency matrix, respectively. The (i, j)th entry
of the adjacency matrix A(i, j) denotes the strength of the
edge between the ith and jth nodes. There exists an edge
between ith and jth nodes if A(i, j) > 0 and the edge pair
(i, j) ∈ E ⇐⇒ A(i, j) 6= 0. In our analysis, we consider only
undirected graphs with symmetric edge-weights or A = A> .
The graph-Laplacian matrix L of G is defined as
L = D − A,
where D is the diagonal degree matrix with ith diagonal
element given by the sum of the elements in the ith row of
A. A vector x = [x(1) x(2) · · · x(M )]> ∈ RM is said to be a
graph signal if x(i) denotes the value of the signal at the ith
node of G. The quadratic form of x with L is given by
X
x> Lx =
A(i, j)(x(i) − x(j))2 .
(i,j)∈E
>

We observe that x Lx is minimized when the signal x takes
the same value across all connected nodes, which agrees
with the intuitive concept of a smooth signal. In general, a
graph signal is said to be smooth or a low-frequency signal
if it has similar values across connected nodes in a graph,
and is said to be a high-frequency signal if it has dissimilar
values across connected nodes, x> Lx being the measure of

similarity. This motivates the use of x> Lx as a constraint in
applications where either the signal x or the graph-Laplacian
L is to be estimated [51], [68]. The eigenvectors of L are
referred to as the graph Fourier transform basis for G, and
the corresponding eigenvalues are referred to as the graph
frequencies. The smaller eigenvalues (the smallest being zero
by construction) are referred to as low frequencies since the
corresponding eigenvectors minimize the quadratic form of L
and vary smoothly over the the nodes. Similarly, the larger
eigenvalues are referred to as the high frequencies. Then,
a smooth graph signal is one which has GFT coefficients
predominantly in the low graph frequencies.
II. K ERNEL R EGRESSION OVER G RAPHS
A. Linear basis model for regression over graphs
Let {xn }N
n=1 denote the set of N input observations. Each
xn is paired with target tn ∈ RM . Our goal is to model the
target tn with yn given by:
yn = W>φ (xn ),

(1)

where φ (xn ) ∈ RK is a function of xn and W ∈ RK×M denotes the regression coefficient matrix. Equation (1) is referred
to as linear basis model for regression, often shortened to linear
regression (cf. Chapters 3 and 6 in [7]). For brevity, we herafter
follow this shortened nomenclature and refer to the outcome
of using (1) as linear regression. Our central assumption is
that target yn is a smooth signal over an underlying graph G
with M nodes. We learn the optimal parameter matrix W by
minimizing the following cost function:
C(W) =

N
X

ktn −yn k22 +αtr(W> W)+β

n=1

N
X

yn> Lyn , (2)

n=1

where regularization coefficients α, β ≥ 0, tr(·) denotes the
trace operator, and kxk2 the `2 norm of x. The cost in
(2) is convex in W, since L is positive semidefinite on
virtue of it being the graph-Laplacian matrix. The choice
of α, β depends on the problem, and in our analysis we
compute these parameters through crossvalidation. The penalty
tr(W> W) = kWk2F ensures that W remains bounded. The
penalty or regularization yn> Lyn enforces yn to be smooth
over G. We define matrices T, Y and Φ as follows:
T = [t1 t2 · · · tN ]> ∈ RN ×M ,
Y = [y1 y2 · · · yN ]> ∈ RN ×M ,
φ(x1 ) φ (x2 ) · · · φ (xN )]> ∈ RN ×K .
Φ = [φ

(3)

Using (1) and (3), the cost function (2) is expressible as (4)
where we use properties of the matrix trace operation. Since
the cost function is quadratic in W, we get the optimal and
unique solution by setting the gradient of C with respect to W
equal to zero. Using the following matrix derivative relations
∂
tr (M1 W) = M>
1,
∂W

∂
>
tr W> M1 WM2
= M>
1 WM2 + M1 WM2 ,
∂W
where M1 and M2 are matrices, and setting

∂C
∂W

= 0 we get

3

C(W)

=
=

=

=

X

X
ktn − W>φ (xn )k22 + α tr(W> W) + β
φ (xn )> WLW>φ (xn )
n
n
!
!
X
X
X
2
>
>
>
ktn k2 − 2 tr
φ (xn ) Wtn + tr
φ (xn ) WW φ (xn ) + α tr(W> W)
n
n
n
!
X
+β tr
φ (xn )> WLW>φ (xn )
n
!
!
X
X
X
2
>
>
>
φ(xn )
ktn k2 − 2 tr
tnφ (xn ) W + tr WW
φ (xn )φ
+ α tr(W> W)
n
n
n
!
X
φ(xn )>
+β tr WLW>
φ (xn )φ
n
X


ktn k22 − 2 tr T> ΦW + tr(W> Φ> ΦW) + α tr(W> W) + β tr W> Φ> ΦWL .

(4)

n

that
>

>

>

−Φ T + Φ ΦW + αW + βΦ ΦWL = 0,
or, Φ> Φ + αIK W + βΦ> ΦWL = Φ> T.

(5)

On vectorizing both sides of (5), we have that


vec(Φ> T) = (IM ⊗(Φ> Φ+αIK ))+(βL⊗Φ> Φ) vec(W),
where vec(·) denotes the standard vectorization operator and
⊗ denotes the Kronecker product operation [69]. Then, the
?

optimal W denoted by W follows the relation:
?

−1
vec(W) = (IM ⊗(Φ> Φ+αIK ))+(βL⊗Φ> Φ) vec(Φ> T).
The predicted target for a new input x is then given by
?

(7) is referred to as a dual representation of (4) in kernel regression literature (cf. Chapter 6 of [7]). Taking the derivative
of C(Ψ) with respect to Ψ and setting it to zero, we get that
(IM ⊗ (K + αIN ))vec(Ψ) + (βL ⊗ K)vec(Ψ) = vec(T), or
[(IM ⊗ (K + αIN )) + (βL ⊗ K)] vec(Ψ) = vec(T),
−1

vec(Ψ) = [(IM ⊗ (K + αIN )) + (βL ⊗ K)]

B = (IM ⊗ (K + αIN )),
C = (βL ⊗ K).

(6)

From (6), it appears that the proposed target prediction
approach requires the explicit knowledge about the function
φ (·). We next show that using the ‘kernel trick’ this explicit
requirement of φ (·) is circumvented and that the target prediction may be done using only the knowledge of inner products
φ (xm )>φ (xn ), ∀m, n, in other words, using the matrix ΦΦ> .
Towards this end, we next discuss a dual representation of the
cost in (2). We hereafter refer to (6) as the output of linear
regression over graphs (LRG).
B. Dual representation of cost using kernel trick
We now use the substitution W = Φ> Ψ and express the
cost function in terms of the parameter Ψ ∈ RN ×M . This
substitution is motivated by observing that on rearranging the
terms in (5), we get that
W = Φ>Ψ ,

where Ψ = α (T − βΦWL − Φ) . On substituting W =
Φ> Ψ in (4), where Ψ becomes the dual parameter matrix
that we wish to learn, and omitting terms that do not depend
on W we get that


C(Ψ) = − 2tr TT ΦΦ> Ψ + tr Ψ> ΦΦ> ΦΦ> Ψ

+ α tr(Ψ> ΦΦ> Ψ) + β tr Ψ> ΦΦ> ΦΦ> ΨL


= − 2tr T> KΨ + tr Ψ> KKΨ

+ α tr(Ψ> KΨ) + β tr Ψ> KKΨL ,
(7)
1

km,n = φ (xm )T φ (xn ).

or

vec(T).

We define the matrices

? >

t = W φ (x).

where K = ΦΦ> ∈ RN ×N denotes the kernel matrix for the
training samples such that its (m, n)th entry is given by

(8)

Then, we have that
vec(Ψ) = (B + C)

−1

vec(T).

(9)

Once Ψ is computed, the kernel regression output for a new
test input x is given by
y

= W> φ (x)
φ(x)
= Ψ> Φφ
= Ψ> k(x)


>
−1
=
mat (B + C) vec(T)
k(x),

(10)

where k(x) = [k1 (x), k2 (x), · · · , kN (x)]> and kn (x) =
φ (xn )>φ (x). mat(·) denotes the reshaping operation of an
argument vector into an appropriate matrix of size N × M by
concatenating subsequent N length sections as columns. In
general, a variety of valid kernel functions may be employed
in kernel regression. Any kernel function k(x, x0 ) is a valid
kernel as long as it expresses an association between x and
x0 and the kernel matrix K is positive semi-definite for
all observation sizes [7]. The Gaussian kernel is one such
popularly employed kernel which we use in our experiments.
We hereafter refer to (10) as the output of kernel regression
over graphs (KRG). We note that KRG is a generalization
over the conventional kernel regression which does not use
any knowledge of the underlying graph structure. On setting

4

β = 0, the KRG output (10) reduces to the conventional KR
output as follows


>
−1
y =
mat (B + C) vec(T)
k(x)


>
−1
=
mat (B) vec(T)
k(x)


>
−1
=
mat (IM ⊗ (K + αIN )) vec(T)
k(x)

>
−1
= mat (I−1
)vec(T)
k(x)
M ⊗ (K + αIN )

>
−1
= (K + αIN ) T k(x)
=

T> (K + αIN )−1 k(x),

(11)

where we have used the Kronecker product equality: vec((K+
−1
αIN )−1 T) = (I−1
)vec(T). Further, we note
M ⊗ (K + αIN )
that KRG reduces to LRG on setting K = ΦΦ> .
C. Interpretation of KRG – a smoothing effect
We next discuss how the output of KRG for training input
is smooth across the training samples {1, · · · , N } and over
the M nodes of the graph. Before proceeding with KRG,
we review a similar property exhibited by the KR output [8].
Using (11) and concatenating KR outputs for the N training
samples, we get that
[y1 y2 . . . yN ] = T> (K + αIN )−1 [k(x1 ) k(x2 ) . . . k(xN )]
or, Y = K(K + αIN )−1 T,
where we use K = [k(x1 ) k(x2 ) . . . k(xN )]. Assuming K
is diagonalizable, let K = UJK UT , where JK and U
denote the eigenvalue and eigenvector matrices of K, respectively. Then, the values of mth (m ∈ {1, · · · , , M })
component of y collected over all time instances ỹ(m) =
[y1 (m), y2 (m), · · · , yN (m)]> , given by the mth column of
Y is equal to
ỹ(m) =

N
X
i=1

θi
u> t̃(m),
θi + α i

where θi and ui denote the ith eigenvalue and eigenvector of
K, respectively, and t̃(m) denotes vector containing the mth
component of the target vector for all training samples (mth
column of T). Thus, we observe that the KR output performs
a shrinkage of t̃(m) along various eigenvector directions for
each m. The contribution from the eigenvectors corresponding to the smaller eigenvalues θi < α are effectively, and
only those corresponding to larger eigenvalues θi > α are
retained. Since the eigenvectors corresponding to larger eigenvalues of K represent smooth variations across observations
{1, · · · , N }, we observe that KR performs a smoothing of tn .
We next show that such is also the case with KRG: KRG
acts as a smoothing filter across both observations and graph
nodes. Using (10) and concatenating KRG outputs, we have
that
[y1 y2 . . . yN ]
or, Y

On vectorizing both sides of (12) , we get that
vec(Y)

= (IM ⊗ K)vec(Ψ)
(a)

−1

= (IM ⊗ K) (B + C)

= KΨ.

(13)

where we have used (9) in (a). Let L be diagonalizable with
the eigendecomposition:
L = VJL V> ,
where JL and V denote the eigenvalue and eigenvector matrices of L, respectively. We also assume that K is diagonalizable
as earlier. Let λi and vi denote the ith eigenvalue and
eigenvector of L, respectively. Then, we have that
=

[v1 v2 · · · vN ] and JL = diag(λ1 , λ2 , · · · , λN ),

U =

[u1 u2 · · · uM ] and JK = diag(θ1 , θ2 , · · · , θM ).

V

Now, using (8), we have
B + C = (IM ⊗ (K + αIN )) + (βL ⊗ K)
= [(VIM V> ) ⊗ (U(JK + αIN )U> )]
+[β(VJL V> ) ⊗ (UJK U> )]
(a)

= [(V ⊗ U)(IM ⊗ (JK + αIN ))(V> ⊗ U> )]
+[β(V ⊗ U)(JL ⊗ JK )(V> ⊗ U> )]
= ZJZ> ,

(14)

where Z = V ⊗ U is the eigenvector matrix and J is the
diagonal eigenvalue matrix given by
J = (IM ⊗ (JK + αIN )) + β(JL ⊗ JK ).

(15)

In (14)(a), we have used the distributivity of the Kronecker
product: (M1 ⊗ M2 )(M3 ⊗ M4 ) = M1 M3 ⊗ M2 M4 where
{Mi }4i=1 are four matrices. We note that J is a diagonal matrix
of size M N . Let J = diag(η1 , η2 , . . . , ηM N ). Then each ηi is
a function of {λj }, {θj }, α, β. On dropping the subscripts for
simplicity, we observe that any eigenvalue ηi has the form
η = (θ + α) + β(λθ).
where θ and λ are the appropriate eigenvalues of K and L,
respectively. Similarly, we have that
(IM ⊗ K)

= (VIM V> ) ⊗ (UJK U> )
= (V ⊗ U)(IM ⊗ JK )(V> ⊗ U> )
= Z(IM ⊗ JK )Z> ,

(16)

and note that (IM ⊗ JK ) is also a diagonal matrix of size
M N . Then, on substituting (15) and (16) in (13), we get that
vec(Y)

−1

= (IM ⊗ K) (B + C) vec(T)
= (Z(IM ⊗ JK )Z> )(ZJ−1 Z> )vec(T)
= (Z(IM ⊗ JK )J−1 Z> )vec(T).

(17)

We note again that (IM ⊗ JK )J−1 is a diagonal matrix with
size M N . Let (IM ⊗JK )J−1 = diag(ζ1 , ζ2 , . . . , ζM N ). Then,
on dropping the subscripts, any ηi is of the form
ζ=

θ
θ
=
.
η
(θ + α) + β(λθ)

From (17), we have that

>

= Ψ [k(x1 ) k(x2 ) . . . k(xN )] (12)

vec(T),

vec (Y) =

M
N
X
i=1

ζ i zi z>
i vec (T) ,

5

where zi are column vectors of Z. In the case when ζi  1,
the component in vec(T) along zi is effectively eliminated.
For most covariance or kernel functions k(·, ·) used in practice,
the eigenvectors corresponding to the largest eigenvalues of K
correspond to the low-frequency components. Similarly, the
eigenvectors corresponding to the smaller eigenvalues of L
are smooth over the graph [9]. We observe that the condition
ζ  1 is achieved when corresponding θ is small and/ or
λ is large. This condition in turn corresponds to effectively
retaining only components in vec(T) which vary smoothly
across the samples {1, · · · , N } as well as over the M nodes
of the graph.
D. Learning an underlying graph
In developing KRG, we have so far assumed that the underlying graph is known apriori in terms of the graph-Laplacian
matrix L. Such an assumption may not hold for many practical
applications since there is not necessarily one best graph for
given networked data. This motivates us to develop a joint
learning approach where we learn L and KRG parameter
matrix W (or its dual representation parameter Ψ). Our goal
in this section is to provide a simple means of estimating
the graph that helps enhance the prediction performance if a
graph is not known apriori. We note that a vast and expanding
body of literature exists in the domain of estimating graphs
from graph signals [50]–[54], [68], [70]–[73], and any of
these approaches may be used as an alternative to the learning
approach proposed in this Section. Nevertheless, we pursue
the approach taken in this Section due to its effectiveness
in terms of ease of implementation and reasonable prediction
performance as demonstrated by experiments.
We propose the minimization of the following cost function
to achieve our goal:
X
C(W, L) =
ktn − yn k22 + α tr(W> W)
n

+β

X

yn> Lyn + ν tr(L> L),

n

where ν ≥ 0. Since our goal is to recover an undirected
graph, we impose appropriate constaints [74]. Firstly, any
non-trivial graph has a graph-Laplacian matrix L which is
symmetric and positive semi-definite [74]. Secondly, the vector
of all ones 1 forms the eigenvector of the graph-Laplacian
corresponding to the zero eigenvalue. Since L = D − A,
we have that the Laplacian being positive semi-definite is
equivalent to the constraint that all the off-diagonal elements
of L are non-positive which is a simpler constraint than the
direct positive semi-definiteness constraint. Then, the solution
to joint estimation of W and L is obtained by solving the
following:
min C(W, L) such that

W,L

L(i, j) ≤ 0 ∀i 6= j,
L = L> , L1 = 0.

The optimization problem (18) is jointly non-convex over W
and L, but convex on W given L and vice-versa. Hence, we
adopt an alternating minimization approach and solve (18) in
two steps alternatingly as follows:

For a given L, solve min C(W) using KRG approach of
W
Section II.
• For
a
given
W,
solve
min C(L) such that L(i, j) ≤ 0 ∀i 6= j, L = L> , and
L
P
L1 = 0. Here C(L) = β n yn> Lyn + ν tr(L> L).
We start the alternating optimization using the intialization
L = 0, which corresponds to KR. In order to keep the
successive L estimates comparable, we scale L such that the
largest eigenvalue modulus is unity at each iteration.
•

III. E XPERIMENTS
We consider application of KRG to various synthesized and
real-world signal examples. We use normalized-mean-squareerror (NMSE) as the performance measure:


EkY − T0 k2F
,
NMSE = 10 log10
EkT0 k2F
where Y denotes the regression output matrix and T0 the
true value of target matrix, that means T0 does not contain
any noise. The expectation operation E(·) in performance
measure is realized as sample average. In the case of realworld examples, we compare the performance of the following
four approaches:
1) Linear regression (LR): km,n = φ (xm )>φ (xn ), where
φ (x) = x and β = 0,
2) Linear
regression
over
graphs
(LRG):
km,n = φ (xm )>φ (xn ) and β 6= 0, where φ (x) = x,
3) Kernel regression (KR): Usingradial basis function

kxm − xn k22
(RBF) kernel km,n = exp −
and
σ2
β = 0, and
4) Kernel regression
Using RBF kernel
 over graphs 2(KRG):

kxm − xn k2
km,n = exp −
and β 6= 0.
σ2
The regularization parameters α, β, and ν for different training
data sizes are found by five-fold crossvalidation on the training
set. While performing crossvalidation, we assume that clean
target vectors are available for the validation
set. The RBF
P
kernel parameter σ is set as σ 2 = 1/N m0 ,n0 kxm0 − xn0 k22
in our experiments, where the summation is made over the
training inputs. This choice of σ is not necessarily optimized
for the training datasets.
A. Regression for synthesized data
We perform regression for synthesized data where the target
vectors are known to be smooth over a specified graph, and the
various target vectors are correlated among each other. A part
of the data is used for training in presence of additive noise,
and our task is to predict for the remaining part of the data,
given the information of the correlations in form of kernels.
In order to generate synthesized data, we use random smallworld graphs from Erdős-Rényi and Barabási-Albert models
[75] with the number of nodes M = 50. We generate a
total of S target vector realizations. We adopt the following data generation strategy: We first pick M independent
vector realizations from an S-dimensional Gaussian vector

6

-2

-4
-6

-15
15

5

10

-8
-10

SNR (dB)

SNR (dB)

(a)

(b)

Fig. 1. Results for synthesized data using Erdős-Rényi graphs. NMSE at
different SNR levels with training data size of N = 50 for (a) training data,
and (b) test data.

source N (0, CS ) where CS is S-dimensional covariance
matrix drawn from inverse Wishart distribution with identity
hyperparameter matrix. We use a highly correlated covariance
matrix CS such that each S-dimensional vector has strongly
correlated components. Thus, we create a data matrix with M
columns such that each column is an S-dimensional Gaussian
vector. Each row vector of the data matrix has a size M and
the row vectors of the data matrix are correlated between each
other. We denote a row vector by r. We then select row vectors
{r} one-by-one and project them on the specified graph to
generate target vectors {t} that are smooth over the graph
while maintaining the correlation between observations, by
solving the following optimization problem:

t = arg min kr − zk22 + z> Lz .
z

We randomly divide the S data samples into training and test
sets of equal size Ntr = Nts = S/2. We define the kernel
function between the ith and jth data samples zi and zj to be
ki,j , CS (i, j),
considering the same kernel for all the graph nodes. The choice
of the kernel is motivated by the assumed generating model.
Given the training set of size Ntr , we choose a subset of N
data samples to make predictions for the Nts test data samples
using the kernel regression over graphs. The training target
vectors are corrupted with additive white Gaussian noise at
varying levels of signal-to-noise ratio (SNR). We repeat our
experiments over 100 realizations of the graphs and noise
realizations. We compare the performance of KRG with KR.
We observe from Figure 1 that for a fixed training data size
of N = 50, KRG outperforms KR by a significant margin at
low SNR levels (below 10dB). As the SNR-level increases,
the normalized MSE of KRG and KR almost coincide. This
is expected since KR is able to make better predictions due to
availablity of cleaner training observations. A similar trend is
observed for also the case of Barabási-Albert graphs, which is
not reported for brevity. In Figure 2, we show the normalized
MSE obtained with KRG and KR on both graph models as
a function of training data size at an SNR-level of 5 dB.
We observe that KRG consistently outperforms KR and that
the gap between NMSE of KRG and KR reduces as training
data size increases as expected. The trend is similar for large
dimensional graphs as well and is not reported here to avoid
repetition.

KR
KRG

-4
-6

10

15

-8

1

10

Training sample size N

1

Training sample size N

(a)

NMSE (dB)

10

-6

-12

-8

5

KR
KRG

(b)
KR
KRG

-4

-6
-8

NMSE (dB)

-10

-2

-4

NMSE (dB)

KR
KRG

NMSE (dB)

KR
KRG

NMSE

NMSE (dB)

-5

KR
KRG

-4
-5
-6

-10
10

1

10

Training sample size N

(c)

1

Training sample size N

(d)

Fig. 2. Results for synthesized data using NMSE at 5 dB SNR level: (a)
for training data, and (b) for test data for Erdős-Rényi graphs. NMSE at 5
dB SNR level: (c) for training data, and (d) for test data for Barabási-Albert
graphs.

B. Temperature prediction for cities in Sweden
We collect temperature measurements from 45 most populated cities in Sweden for a period of three months from
September to November 2017. The data is available publicly
from the Swedish Meteorological and Hydrological Institute
[76]. We consider the target vector tn to be the temperature
measurement of a particular day, and xn to be the temperature
measurements (in degree celsius units) from the previous day.
We have 90 input-target data pairs in total, divided into training
set and test set of sizes Ntr = 60 and Nts = 30, respectively.
Let dij denote the geodesic distance between cities i and j in
kilometres, ∀i, j ∈ {1, · · · , 45}. We construct the adjacency
matrix by setting
!
d2ij
A(i, j) = exp − P
.
2
i,j dij
In order to remove self loops, the diagonal of A is set to
zero. We simulate noisy training data by adding zero-mean
white Gaussian noise to the targets. For each N , we compute
the NMSE by averaging over 50 different random training
subsets of N drawn from full training set of size Ntr . In
Figure 3, we show the NMSE for test set at SNR-levels of
5 dB and 0 dB. We observe that KRG outperforms other
regression methods by a significant margin, particularly at low
sample sizes N . In order to simulate a more realistic noisy
training scenario, we further consider two different cases of
noisy training: missing data and large perturbations. In the
case of missing data experiment, we zero out the target values
at randomly chosen five nodes for each training sample and
apply the various regression approches. In the perturbation
experiment, we chose five of the nodes randomly for each
training sample and perturb the signal value to five times the
original values. The NMSE obtained for LR, LRG, KR, and
KRG as a function of N is shown in Figures 3(d)-(e). Thus, we
observe that graph-structure improves robustness of regression
in terms of prediction performance for LR and KR.

7

12

70 ° N

10

-6
65 ° N

6
4
2

60 ° N

NMSE (dB)

8

LR
LRG
KR
KRG

-10

0

-15

-2
10 ° E

(a)

LR
LRG
KR
KRG

LR
LRG
KR
KRG

-8

NMSE (dB)

NMSE (dB)

15

-10
-15

-10
-12
-14

5

15

30

45 60

LR
LRG
KR
KRG

-5

-10

-15
5

45 60

30

Training sample size N
(c)

(b)

-5

-20

5

-4

°
25 E

°
20 E

°

15 E

NMSE (dB)

55 ° N
5° E

15

30

45 60

5

Training sample size N

Training sample size N
(d)

15

30

45 60

Training sample size N

(e)

(f)

Fig. 3. Results for the temperature data of 45 cities in Sweden. (a) Map with cities, (b) True temperature measurement signal for a particular day in the
dataset shown in units of degree celsius (The graph is the same as that in subplot (a) but without the underlying map.), (c) NMSE for test data with additive
white Gaussian noise at 5dB SNR level, (d) NMSE for test data with additive white Gaussian noise at 0dB SNR level, (e) NMSE test data at with missing
samples, and (f) NMSE test data at with large perturbations.

LR
LRG-learnt L
KR
KRG-learnt L

NMSE (dB)

-8
-10

1
0.8

10

0.6

20
0.4

-14
-16

30

15

30

Training sample size N

(a)

45 60

40

0

10

20

(b)

30

40

0.6

20
0.4

0.4

30

30
0.2

5

0.8

0.8 10

10

0.6

20

-12

1

1

0.2

0.2

40

0

10

20

(c)

30

40

40

0

10

20

(d)

30

40

Fig. 4. Results with graph learning for temperature data of cities in Sweden, for training data samples with additive white Gaussian noise at 5dB SNR-level
for N = 45. (a) Normalized MSE for test data, (b) Geodesic Laplacian, (c) Estimated Laplacian from LRG, (d) Estimated Laplacian from KRG.

We also consider the joint learning of the graph from
training data using the approach described in Section II-D.
We perform the training at an SNR level of 5 dB. We
find experimentally that the alternating optimization based
graph learning algorithm converges typically after five to ten
iterations. The NMSE for training and test data obtained using
geodesic L and the iteratively learnt L for both LRG and KRG
are shown in Figure 4 for N = 45. The results show that
our alternating optimization based graph learning algorithm is
useful in the case when we do not know the underlying graph
a-priori. We also observe that the estimated Laplacian matrices
from both LRG and KRG are close to each other and are in
reasonable agreement with that of the geodesic graph. This
validates our intuition that the graph signal holds sufficient
information to both infer the underlying graph structure and
perform predictions.

C. Prediction on EEG signals
We repeat the prediction experiment on 64-channel EEG
signal data taken from the Physionet database [77]. The signals
are obtained by placing electrodes at various locations of a
subject’s head, as shown in the schematic in Figure 5 (a).
At each time instant, we consider a subset of 16 nodes as
the input and the target graph signal to be the M = 48dimensional vector with components given by signal values at
the remaining electrodes. Our goal is then to predict the signal
at the M electrodes for the rest of the electrodes given those
at the input electrodes. Since the graph is not given apriori,
we construct the graph for the M electrodes from a part of
noiseless data using Gaussian kernels as follows:
!
krk − rl k22
A(i, j) = exp − P
∀k, l ∈ {1, · · · , M },
2
k,l krk − rl k2

8

-1.25
23

24

AF3

25

AFZ
27

26

F7

F5

30

31

FT7

FC5

39

29

28

F3

F1

FZ

F2

F4

32

33

34

35

36

FC3

1

AF8

AF4

FC1

FCZ

38

37

FT8

FC6

FC4

FC2

F8

F6

40

5

6

7

4

T9

T7

C5

C3

C1

CZ

C2

C4

C6

T8

T10

43

41

8

9

10

11

12

13

14

42

44

CPZ

CP4

CP6

20

21

TP7
45

P7

2

3

CP5

CP3

CP1

15

16

17

P5
48

P3

P1

49

50

18

CP2
19

PZ

P2

51

52

P4
53

P6
54

47

TP8

NMSE (dB)

22

AF7

PO3
57

56

POZ
58

O1
61

PO4
59

-2

-4

-2
-3
-4
-5

-6

P8

PO8
60

O2
OZ

LR
LRG
KR
KRG

46

55

PO7

-0.5
-1

LR
LRG
KR
KRG

FP2

NMSE (dB)

FPZ

FP1

4

63

62

8

16

32

64

4

128

8

(a)

16

32

64

128

Training sample size N

Training sample size N

IZ
64

(b)

(c)

Fig. 5. Results for EEG data. (a) Schematic showing positions of 64 electrode channels, (b) NMSE for test data at 5dB SNR-level. (c) NMSE for test data
at 0dB SNR-level.

(a)
-1

-4
-6
-8
5

10

15

Training sample size N

(b)

25 30

LR
LRG
KR
KRG

-2
-4
-6
-8

-2

NMSE (dB)

LR
LRG
KR
KRG

NMSE (dB)

NMSE (dB)

-2.5

LR
LRG with learnt L
KR
KRG with learnt L

-4
-6
-8

5

10

15

Training sample size N
(c)

25 30

5

10

15

25 30

Training sample size N
(d)

Fig. 6. Results for ETEX data. (a) Schematic showing positions of the ground stations, (b) NMSE for test data at 5dB SNR-level, (c) NMSE for test data at
0dB SNR-level, (d) NMSE with estimated Laplacian matrices at 5dB SNR.

where rl denotes the vector with first 1000 time-samples for
the lth electrode. Thus, we use the first 1000 time-samples
to generate the diffusion or Gaussian kernel graph for the
electrodes. Since EEG signals are highly nonstationary [78],
we consider the dataset corresponding to 20 seconds consisting
of approximately 256 input-target data samples, divided randomly into training and test data of 128 samples each. We then
consider training with noisy targets corrupted with additive
white Gaussian noise at SNR levels of 0 dB and 5 dB. The
NMSE is computed by averaging over 50 different random
training subsets of size N drawn from the full training set of
size Ntr = 128. In Figure 5, we show the NMSE for LR, LRG,
KR, and KRG for different training data sizes. We observe
once again that LRG and KRG consistently outperform LR
and KR, respectively, particularly at small training data sizes.

D. Prediction for atmospheric tracer diffusion data
We next perform prediction experiments on the European
Tracer Experiment data (ETEX) which consists of measurements from diffusion of perfluorocarbon tracers released into
the atmosphere starting from a fixed location (Rennes in
France) [79]. This is related to environmental pollution. The
observations are made at 168 ground stations over Europe
and consist of two experiments carried out over 72 hours
giving 30 measurements each. We consider the problem of
predicting the tracer concentrations on one half (84) of the
total ground stations, given the concentrations at the remaining
locations. The target signal tn is then a graph signal over a
graph of M = 84 nodes, the input vector xn being also of the
same dimension, where n denotes the measurement index. A
schematic diagram of the ground station locations is shown in
Figure 6 (a). The red markers with edges correspond to the
output nodes, whereas the rest of the markers denote the input.

9

-15
-20
-25
4

16

32

64

Training sample size N
(a)

145

LR
LRG
KR
KRG

-5
-10
-15
-20
-25

4

16

32

64

Training sample size N
(b)

145

NMSE (dB)

-10

NMSE (dB)

NMSE (dB)

-5
LR
LRG
KR
KRG

LR
LRG with learnt L
KR
KRG with learnt L

-10
-15
-20
-25
-30

4

16

32

64

145

Training sample size N
(c)

Fig. 7. Results for Cerebellum data. (a) NMSE for test data at 10dB SNR-level, (b) NMSE for test data at 0dB SNR-level, (c) NMSE for test data at with
learnt Laplacian at 10dB SNR-level.

We generate noisy training samples by adding white Gaussian
noise at different SNR levels. As with the case of temperature
data considered in Section III-B, we consider the geodesic
distance based graph. We randomly divide the total dataset of
60 samples equally into training and test datasets. We compute
the NMSE for the different approaches by averaging over 50
different randomly drawn training subsets of size N from
the full training set of size Nts = 30. We report NMSE as
a function of N in Figures 6(b)-(c) which show that graph
structure enhances the prediction performance under noisy and
low sample size conditions. We also perform LRG and KRG
with graph-learning and the NMSE is shown in Figure 6(d).
The results again show that our alternating optimization based
graph learning algorithm is useful in the absence of a-priori
knowledge about underlying graph.

is analytically tractable. The development of kernel regression
over graphs is useful for several practical signals and systems.
Our experiments show that it is particularly beneficial for
applications where available training data is limited in size
and noisy. We have also seen that it is possible to learn
underlying graph from data. A practical limitation of the
proposed approach is that which it shares with all kernel-based
approaches that it involves inversion of a large system matrix.

E. Prediction for fMRI data of the cerebellum region
Finally, we consider prediction of voxel intensities in the
functional magnetic resonance imaging (fMRI) data obtained
for the cerebellum region of the brain. The data is available
publicly at https://openfmri.org/dataset/ds000102. Each pixel
is considered as a node of the graph and the intensity to
be the signal. The original data graph is of dimension 4000
obtained by mapping the 4000 cerebellum voxels anatomically
following the atlas template [80], [81]. In our analysis, we
consider the first 100 voxels from 4000 voxels due to ease of
experiments in the sense of computational complexity. We use
the intensity values at the first 10 voxels as input x ∈ R10 to
make predictions for the output t ∈ R90 comprising the voxel
instensities at the remaining 90 voxels. The dataset consists
of 290 input-target data points. We use one half of the data
for training and the other half for testing. We construct noisy
training data at a signal-to-noise (SNR) levels of 10 dB and 0
dB. For test data, the NMSE averaged over 50 different random
selections of training data of size N from the full training
dataset is reported in Figure 7(a)-(c). We observe that LRG
and KRG have superior performance over their conventional
counterparts particularly for small N and low SNR-level. We
also see that algorithm for estimation of underlying graph does
provide a reasonable performance.

R EFERENCES

IV. C ONCLUSIONS
We conclude that extension of a standard kernel regression with inclusion of graph-Laplacian based regularization

V. R EPRODUCIBLE RESEARCH
In the spirit of reproducible research, all the codes
relevant to the experiments in this article are made available
at
https://www.researchgate.net/profile/Arun Venkitaraman
and
https://www.kth.se/ise/research/reproducibleresearch1.433797.

[1] C. Cortes and V. Vapnik, “Support-vector networks,” Machine Learning,
vol. 20, no. 3, pp. 273–297, 1995.
[2] M. Seeger, “Gaussian processes for machine learning,” Int. J. Neural
Syst., vol. 14, no. 02, pp. 69–106, 2004.
[3] Y. B. L. Yann and G. Hinton., “Deep learning,” Nature, vol. 521, no.
7553, pp. 436–444, 2015.
[4] Y. Cho and L. K. Saul, “Kernel methods for deep learning,” in Adv.
Neural Inform. Process. Syst. Curran Associates, Inc., 2009, pp. 342–
350.
[5] G.-B. Huang, “An insight into extreme learning machines: Random
neurons, random features and kernels,” Cognitive Computation, vol. 6,
pp. 376–390, 2014.
[6] A. Iosifidis, A. Tefas, and I. Pitas, “On the kernel extreme learning
machine classifier,” Patt. Recognition Lett., vol. 54, pp. 11 – 17, 2015.
[7] C. M. Bishop, Pattern Recognition and Machine Learning (Information
Science and Statistics). Secaucus, NJ, USA: Springer-Verlag New York,
Inc., 2006.
[8] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine
Learning (Adaptive Computation and Machine Learning). The MIT
Press, 2005.
[9] D. I. 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 Process. Mag., vol. 30, no. 3, pp. 83–98, 2013.
[10] A. Sandryhaila and J. M. F. Moura, “Discrete signal processing on
graphs,” IEEE Trans. Signal Process., vol. 61, no. 7, pp. 1644–1656,
2013.
[11] ——, “Big data analysis with signal processing on graphs: Representation and processing of massive data sets with irregular structure,” IEEE
Signal Process. Mag., vol. 31, no. 5, pp. 80–90, 2014.
[12] ——, “Discrete signal processing on graphs: Frequency analysis,” IEEE
Trans. Signal Process., vol. 62, no. 12, pp. 3042–3054, 2014.
[13] D. I. Shuman, B. Ricaud, and P. Vandergheynst, “A windowed graph
Fourier transform,” IEEE Statist. Signal Process. Workshop (SSP), pp.
133–136, Aug 2012.

10

[14] S. K. Narang and A. Ortega, “Local two-channel critically sampled filterbanks on graphs,” Proc. IEEE Int. Conf. Image Process. (ICIP), pp.
333–336, 2010.
[15] ——, “Perfect reconstruction two-channel wavelet filter banks for graph
structured data,” IEEE Trans. Signal Process., vol. 60, no. 6, pp. 2786–
2799, 2012.
[16] ——, “Compact support biorthogonal wavelet filterbanks for arbitrary
undirected graphs,” IEEE Trans. Signal Process., vol. 61, no. 19, pp.
4673–4685, 2013.
[17] R. R. Coifman and M. Maggioni, “Diffusion wavelets,” Appl. Computat.
Harmonic Anal., vol. 21, no. 1, pp. 53–94, 2006.
[18] D. Ganesan, B. Greenstein, D. Estrin, J. Heidemann, and R. Govindan,
“Multiresolution storage and search in sensor networks,” ACM Trans.
Storage, vol. 1, no. 3, pp. 277–315, 2005.
[19] D. K. Hammond, P. Vandergheynst, and R. Gribonval, “Wavelets on
graphs via spectral graph theory,” Appl. Computat. Harmonic Anal.,
vol. 30, no. 2, pp. 129–150, 2011.
[20] R. Wagner, V. Delouille, and R. Baraniuk, “Distributed wavelet denoising for sensor networks,” Proc. 45th IEEE Conf. Decision Control,
pp. 373–379, 2006.
[21] D. I. Shuman, B. Ricaud, and P. Vandergheynst, “Vertex-frequency
analysis on graphs,” Applied and Computational Harmonic Analysis,
vol. 40, no. 2, pp. 260 – 291, 2016.
[22] D. I. Shuman, M. J. Faraji, and P. Vandergheynst, “A multiscale pyramid
transform for graph signals,” IEEE Trans. Signal Process., vol. 64, no. 8,
pp. 2119–2134, April 2016.
[23] O. Teke and P. P. Vaidyanathan, “Extending classical multirate signal
processing theory to graphs-part i: Fundamentals,” IEEE Trans. Signal
Process., vol. 65, no. 2, pp. 409–422, Jan 2017.
[24] ——, “Extending classical multirate signal processing theory to graphspart ii: M-channel filter banks,” IEEE Trans. Signal Process., vol. 65,
no. 2, pp. 423–437, Jan 2017.
[25] A. Venkitaraman, S. Chatterjee, and P. Handel, “On Hilbert
transform of signals on graphs,” Proc. Sampling Theory Appl.,
2015. [Online]. Available: http://w.american.edu/cas/sampta/papers/
a13-venkitaraman.pdf
[26] A. Venkitaraman, S. Chatterjee, and P. Händel, “Hilbert transform,
analytic signal, and modulation analysis for graph signal processing,”
CoRR, vol. abs/1611.05269, 2016. [Online]. Available: http://arxiv.org/
abs/1611.05269
[27] N. Tremblay, G. Puy, P. Borgnat, R. Gribonval, and P. Vandergheynst,
“Accelerated spectral clustering using graph filtering of random signals,”
IEEE Int. Conf. Acoust. Speech Signal Process., 2016.
[28] N. Tremblay and P. Borgnat, “Joint filtering of graph and graph-signals,”
Asilomar Conference on Signals, Systems and Computers, pp. 1824–
1828, Nov 2015.
[29] S. Chen, R. Varma, A. Sandryhaila, and J. Kovacevic, “Discrete signal
processing on graphs: Sampling theory,” IEEE Transactions on Signal
Processing, vol. 63, no. 24, pp. 6510–6523, Dec 2015.
[30] F. Gama, A. G. Marques, G. Mateos, and A. Ribeiro, “Rethinking
sketching as sampling: Linear transforms of graph signals,” Proc.
Asilomar Conf. Signals, Systems Comput., pp. 522–526, Nov 2016.
[31] S. P. Chepuri and G. Leus, “Subsampling for graph power spectrum
estimation,” Proc. IEEE Sensor Array Multichannel Signal Process.
Workshop (SAM), pp. 1–5, July 2016.
[32] S. Chen, R. Varma, A. Singh, and J. Kovačević, “Signal recovery on
graphs: Fundamental limits of sampling strategies,” IEEE Trans. Signal
Inform. Process. over Networks, vol. 2, no. 4, pp. 539–554, Dec 2016.
[33] A. G. Marques, S. Segarra, G. Leus, and A. Ribeiro, “Sampling of graph
signals with successive local aggregations,” IEEE Trans. Signal Process.,
vol. 64, no. 7, pp. 1832–1843, April 2016.
[34] H. Q. Nguyen and M. N. Do, “Downsampling of signals on graphs via
maximum spanning trees,” IEEE Trans. Signal Process., vol. 63, no. 1,
pp. 182–191, Jan 2015.
[35] M. Tsitsvero, S. Barbarossa, and P. D. Lorenzo, “Signals on graphs:
Uncertainty principle and sampling,” IEEE Trans. Signal Processing,
vol. 64, no. 18, pp. 4845–4860, 2016.
[36] A. Anis, A. Gadde, and A. Ortega, “Efficient sampling set selection for
bandlimited graph signals using graph spectral proxies,” IEEE Trans.
Signal Process., vol. 64, no. 14, pp. 3775–3789, July 2016.
[37] N. Shahid, V. Kalofolias, X. Bresson, M. Bronstein, and P. Vandergheynst, “Robust principal component analysis on graphs,” IEEE Int.
Conf. Comput. Vision (ICCV), pp. 2812 – 2820, 2015.
[38] N. Shahid, N. Perraudin, V. Kalofolias, G. Puy, and P. Vandergheynst,
“Fast robust pca on graphs,” IEEE Journal of Selected Topics in Signal
Processing, vol. 10, no. 4, pp. 740–756, June 2016.

[39] D. Thanou, D. I Shuman, and P. Frossard, “Learning parametric dictionaries for signals on graphs,” IEEE Trans. Signal Process., vol. 62,
no. 15, pp. 3849–3862, 2014.
[40] D. Thanou and P. Frossard, “Multi-graph learning of spectral graph
dictionaries,” in IEEE Intl. Conf. Acoust. Speech Signal Processing,
April 2015, pp. 3397–3401.
[41] ——, “Multi-graph learning of spectral graph dictionaries,” IEEE Int.
Conf. Acoust. Speech Signal Process. (ICASSP), pp. 3397–3401, April
2015.
[42] D. Thanou, D. I. Shuman, and P. Frossard, “Learning parametric
dictionaries for signals on graphs,” IEEE Trans. Signal Process., vol. 62,
no. 15, pp. 3849–3862, Aug 2014.
[43] Y. Yankelevsky and M. Elad, “Dual graph regularized dictionary learning,” IEEE Trans. Signal Info. Process. Networks, vol. 2, no. 4, pp.
611–624, Dec 2016.
[44] B. Girault, “Stationary graph signals using an isometric graph translation,” Proc. Eur. Signal Process. Conf. (EUSIPCO), pp. 1516–1520, Aug
2015.
[45] S. Segarra, A. G. Marques, G. Leus, and A. Ribeiro, “Stationary graph
processes: Nonparametric spectral estimation,” Proc. IEEE Sensor Array
Multichannel Signal Process. Workshop (SAM), pp. 1–5, July 2016.
[46] N. Perraudin and P. Vandergheynst, “Stationary signal processing
on graphs,” CoRR, vol. abs/1601.02522, 2016. [Online]. Available:
http://arxiv.org/abs/1601.02522
[47] B. Girault, P. Goncalves, and E. Fleury, “Translation on graphs: An
isometric shift operator,” IEEE Signal Process. Lett., vol. 22, no. 12,
pp. 2416–2420, Dec. 2015.
[48] B. Girault, “Stationary graph signals using an isometric graph translation,” Proc. Eur. Signal Process. Conf. (EUSIPCO), 2015.
[49] P. Berger, G. Hannak, and G. Matz, “Graph signal recovery via primaldual algorithms for total variation minimization,” IEEE J. Selected
Topics Signal Process., vol. 11, no. 6, pp. 842–855, Sept 2017.
[50] X. Dong, D. Thanou, P. Frossard, and P. Vandergheynst, “Learning
laplacian matrix in smooth graph signal representations,” IEEE Trans.
Signal Process., vol. 64, no. 23, pp. 6160–6173, Dec 2016.
[51] C. Hu, L. Cheng, J. Sepulcre, G. E. Fakhri, Y. M. Lu, and Q. Li, “A
graph theoretical regression model for brain connectivity learning of
alzheimer’s disease,” in 2013 IEEE 10th International Symposium on
Biomedical Imaging, April 2013, pp. 616–619.
[52] S. I. Daitch, J. A. Kelner, and D. A. Spielman, “Fitting a graph to vector
data,” in Proceedings of the 26th Annual International Conference on
Machine Learning, ser. ICML ’09. New York, NY, USA: ACM, 2009,
pp. 201–208.
[53] N. Leonardi and D. V. D. Ville, “Wavelet frames on graphs defined
by fMRI functional connectivity,” in IEEE International Symposium on
Biomedical Imaging, March 2011, pp. 2136–2139.
[54] N. Leonardi, J. Richiardi, M. Gschwind, S. Simioni, J.-M. Annoni,
M. Schluep, P. Vuilleumier, and D. V. D. Ville, “Principal components
of functional connectivity: A new approach to study dynamic brain
connectivity during rest,” NeuroImage, vol. 83, pp. 937–950, 2013.
[55] R. I. Kondor and J. Lafferty, “Diffusion kernels on graphs and other
discrete structures,” Proc.of the ICML, pp. 315–322, 2002.
[56] A. J. Smola and R. Kondor, Kernels and Regularization on Graphs.
Berlin, Heidelberg: Springer Berlin Heidelberg, 2003, pp. 144–158.
[57] M. Belkin, I. Matveeva, and P. Niyogi, “Tikhonov regularization and
semi-supervised learning on large graphs,” Proc. IEEE Intl. Conf.
Acoust.Speech Signal Process., May 2004.
[58] ——, “Regularization and semi-supervised learning on large graphs,”
in Proc. 17th Annual Conf. Learning Theory, J. Shawe-Taylor and
Y. Singer, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2004,
pp. 624–638.
[59] A. Argyriou, M. Herbster, and M. Pontil, “Combining graph laplacians
for semi-supervised learning,” NIPS’05, pp. 67–74, 2005.
[60] M. Belkin, P. Niyogi, and V. Sindhwani, “Manifold regularization: A
geometric framework for learning from labeled and unlabeled examples,”
J. Mach. Learn. Res., vol. 7, pp. 2399–2434, Dec. 2006.
[61] K. Zhang, L. Lan, J. T. Kwok, S. Vucetic, and B. Parvin, “Scaling up
graph-based semisupervised learning via prototype vector machines,”
IEEE Trans. Neural Networks Learning Systems, vol. 26, no. 3, pp.
444–457, March 2015.
[62] K. I. Kim, J. Tompkin, H. Pfister, and C. Theobalt, “Context-guided
diffusion for label propagation on graphs,” in IEEE Intl. Conf. Computer
Vision (ICCV), Dec 2015, pp. 2776–2784.
[63] D. Romero, M. Ma, and G. B. Giannakis, “Kernel-based reconstruction
of graph signals,” IEEE Trans. Signal Process., vol. 65, no. 3, pp. 764–
778, Feb 2017.

11

[64] ——, “Estimating signals over graphs via multi-kernel learning,” in
IEEE Statist. Signal Process. Workshop (SSP), June 2016, pp. 1–5.
[65] V. N. Ioannidis, D. Romero, and G. B. Giannakis, “Kernel-based
reconstruction of space-time functions via extended graphs,” in Asilomar
Conf. Signals Syst. Comput., Nov 2016, pp. 1829–1833.
[66] D. Romero, V. N. Ioannidis, and G. B. Giannakis, “Kernel-based
reconstruction of space-time functions on dynamic graphs,” CoRR, vol.
abs/1612.03615, 2016. [Online]. Available: http://arxiv.org/abs/1612.
03615
[67] Y. Shen, B. Baingana, and G. B. Giannakis, “Kernel-based structural
equation models for topology identification of directed networks,” IEEE
Trans. Signal Process., vol. 65, no. 10, pp. 2503–2516, May 2017.
[68] X. Dong, D. Thanou, P. Frossard, and P. Vandergheynst, “Learning
graphs from signal observations under smoothness prior,” CoRR, vol.
abs/1406.7842, 2014. [Online]. Available: http://arxiv.org/abs/1406.7842
[69] C. F. V. Loan, “The ubiquitous Kronecker product,” J. Computat. Appl.
Mathematics, vol. 123, no. 1–2, pp. 85–100, 2000.
[70] P. K. Shivaswamy and T. Jebara, Laplacian Spectrum Learning. Berlin,
Heidelberg: Springer Berlin Heidelberg, 2010, pp. 261–276.
[71] X. Dong, D. Thanou, P. Frossard, and P. Vandergheynst, “Laplacian
matrix learning for smooth graph signal representation,” in IEEE Intl.
Conf. Acoust. Speech Signal Processing, April 2015, pp. 3736–3740.
[72] H. E. Egilmez, E. Pavez, and A. Ortega, “Graph learning from data under
laplacian and structural constraints,” IEEE J. Selected Topics Signal
Process., vol. 11, no. 6, pp. 825–841, Sept 2017.
[73] V. Kalofolias, “How to learn a graph from smooth signals,” in
Proc. Intl. Conf. Artificial Intell. Statistics, ser. Proc. Machine
Learning Research, A. Gretton and C. C. Robert, Eds., vol. 51.
PMLR, 09–11 May 2016, pp. 920–929. [Online]. Available: http:
//proceedings.mlr.press/v51/kalofolias16.html
[74] F. R. K. Chung, Spectral Graph Theory. AMS, 1996.
[75] M. E. J. Newman, Networks: An Introduction. Oxford University Press,
2010.
[76] Swedish meteorological and hydrological institute (smhi). [Online].
Available: http://opendata-download-metobs.smhi.se/
[77] A. L. Goldberger, L. A. N. Amaral, L. Glass, J. M. Hausdorff, P. C.
Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C.-K. Peng, and H. E.
Stanley, “Physiobank, physiotoolkit, and physionet,” vol. 101, no. 23,
pp. e215–e220, 2000.
[78] S. Sanei and J. A. Chambers, EEG Signal Processing. Wiley, 2007.
[79] ”European Tracer Experiment (ETEX),” 1995, [online]. [Online].
Available: https://rem.jrc.ec.europa.eu/RemWeb/etex/.
[80] H. Behjat, N. Leonardi, L. S. a, and D. V. D. Ville, “Anatomicallyadapted graph wavelets for improved group-level fmri activation mapping,” NeuroImage, vol. 123, pp. 185 – 199, 2015.
[81] J. Diedrichsen, J. H. Balsters, J. Flavell, E. Cussans, and N. Ramnani, “A
probabilistic mr atlas of the human cerebellum,” NeuroImage, vol. 46,
no. 1, pp. 39 – 46, 2009.

