Unified Statistical Theory of Spectral Graph Analysis
Subhadeep Mukhopadhyay
Department of Statistical Science, Temple University
Philadelphia, Pennsylvania, 19122, U.S.A.
Dedicated to the beloved memory of Emanuel (Manny) Parzen.

arXiv:1602.03861v4 [] 20 Sep 2016

Abstract
The goal of this paper is to show that there exists a simple, yet universal1 statistical
logic of spectral graph analysis by recasting it into a nonparametric function estimation
problem. The prescribed viewpoint appears to be good enough to accommodate most
of the existing spectral graph techniques as a consequence of just one single formalism
and algorithm.

Keywords and phrases: Nonparametric spectral graph analysis; Graph correlation density field (GraField); Empirical and smoothed spectral graph analysis; High-dimensional
discrete data smoothing.

Contents
1 Introduction
1.1 Spectral Graph Analysis: A Practitioner’s Guide
1.2 Previous Theoretical Treatments . . . . . . . . .
1.3 A Set of Questions . . . . . . . . . . . . . . . . .
1.4 Unified Construction Principle . . . . . . . . . .

.
.
.
.

2
2
3
3
4

2 Fundamentals of Statistical Spectral Graph Analysis
2.1 Graph Correlation Density Field . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Karhunen-Loéve Representation of Graph . . . . . . . . . . . . . . . . . . . . . . .
2.3 Nonparametric Spectral Approximation Theory . . . . . . . . . . . . . . . . . . . .

5
5
7
9

3 Empirical Spectral Graph Analysis
3.1 Laplacian Spectral Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Diffusion Map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Modularity Spectral Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11
11
13
14

4 Smoothed Spectral Graph Analysis
4.1 High-dimensional Undersampled Regime
4.2 Spectral Smoothing . . . . . . . . . . . .
4.3 Type-I Regularized Graph Laplacian . .
4.4 Type-II Regularized Graph Laplacian .
4.5 Google’s PageRank Method . . . . . . .
4.6 Other Generalizations . . . . . . . . . .

15
15
15
16
17
17
18

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

5 Application to Graph Regression

19

6 Concluding Remarks

22

1

By universal, I mean the theory should not be specific to any particular spectral graph analysis technique
but be meaningfully applicable to a class of models.

1

1

Introduction

Spectral graph data analysis is undoubtedly the most favored technique for graph data
analysis, both in theory and practice. It has been extraordinarily successful in many applications including data mining, web search, quantum computing, computer vision, image
segmentation, and others. Over the last several decades, extensive efforts have been made
by researchers and practitioners to develop a suit of spectral graph analysis techniques (e.g.,
Laplacian, Modularity, Diffusion map, regularized Laplacian, Google PageRank model etc.)
with increasing sophistication and specialization. However, no single algorithm can be regarded as a panacea for dealing with the evolving complexities of modern graphs. Therefore,
the most important and pressing question for the field today appears to be whether we can
develop a unifying language to establish “bridges” between a wide variety of spectral graph
analysis techniques, and thus providing logically connected means for reaching different
ends. Undoubtedly, any such formulation would be of great theoretical significance and
practical value, that will ultimately provide the applied data scientists clear guidance and a
systematic strategy for selecting the proper spectral tools to arrive at a confident conclusion.
To that end, this work attempts to unify the theories of spectral graph analysis by purely
statistical means. That is, we seek a constructive (rather than confirmatory) theoretical
framework integrating classical and modern spectral graph analysis algorithms, which to
date have been viewed as distinct methods.

1.1

Spectral Graph Analysis: A Practitioner’s Guide

The way spectral graph analysis is currently taught and practiced can be summarized as
follows (also known as spectral heuristics):
1. Let G = (V, E) denotes a (possibly weighted) undirected graph with a finite set of
vertices |V | = n, and a set of edges E. Represent the graph using weighted adjacency
matrix A where A(x, y; G) = wxy if the nodes x and y are connected by an edge and
0 otherwise; weights are non-negative and symmetric.
2. Define “suitable” spectral graph matrix (also known as graph “shift” operator). Most
popular and successful ones are listed below:
• L = D−1/2 AD−1/2 ;

Chung (1997)

• B = A − N −1 ddT ;

Newman (2006)

• T = D−1 A;

Coifman and Lafon (2006)
−1/2

• Type-I Reg. Lτ = Dτ

−1/2

• Type-II Reg. Lτ = Dτ

−1/2

A Dτ

;

Chaudhuri et al. (2012)

−1/2

Aτ Dτ

; Amini et al. (2013)

• Google’s PageRank Tα = αD−1 A + (1 − α)F ; (Brin and Page, 1999)
D = diag(d1 , . . . , dn ) ∈ Rn×n , di denotes the degree of a node, τ > 0 regularization
P
parameter, F ∈ Rn×n with all entries 1/n, and N = 2|E| = x,y A(x, y).
2

3. Perform spectral decomposition of the matrix selected at step 2, whose eigenvectors
form an orthogonal basis of Rn . Spectral graph theory seeks to understand the interesting properties and structure of a graph by using leading nontrivial eigenvectors and
eigenvalues, first recognized by Fiedler (1973), which provide the basis for performing
harmonic analysis on graphs.
4. Compute graph Fourier transform (GFT) by expanding functions (or signals) defined
over the vertices of a graph as a linear combination of the selected eigenbasis (from step
3) and carry out learning tasks such as regression, clustering, classification, smoothing,
kriging, etc.

1.2

Previous Theoretical Treatments

There has been a flood of publications on statistical spectral graph theory, most of which
have focused on asymptotic justifications of different techniques (mentioned in the previous
section). See, in particular, those based on Laplacian (Von Luxburg et al., 2008, Belkin
and Niyogi, 2001), modularity (Davis and Sethuraman, 2016), diffusion map (Coifman and
Lafon, 2006, Coifman et al., 2005), and regularized Laplacian (Chaudhuri et al., 2012, Le
et al., 2015) methods. More refined theoretical results in the context of a specific parametric
graph model (such as stochastic block model or its variants) are discussed in Joseph and Yu
(2016), Rohe et al. (2011), Sarkar et al. (2015) and Qin and Rohe (2013).
What is this paper not about? It is important to draw a distinction between the theory
that purports to provide the constructions of different spectral graph techniques and, the
theory that does not. This article focuses on the former. Instead of statistical confirmation,
here our main interest lies in understanding the statistical origins of various spectral graph
techniques. This is mostly ill-understood and neglected territory. The current theoretical
treatments do not give us any clues on this issue.

1.3

A Set of Questions

Despite many advances, the theory and practice of spectral graph analysis are still very
much compartmentalized and unsystematic. The graph signal processing engineers are often
left bewildered by the vastness of the existing literature and huge diversity of methods
(developed by the machine learning community, applied harmonic analysts, physicists, and
statisticians). This has led to a need to develop a broader perspective on this topic, lest
we be guilty of not seeing the forest for the trees. The question how different spectral
graph techniques can naturally originate from some underlying basic principle, plays the
key role in this regard by clarifying the mystery of why and when to use them. The reality
is we still lack general tools and techniques to attack this question in a statistical way.
Given the very different character of the existing spectral models, it seems that new kinds
of abstractions need to be formulated to address this fundamental challenge. What would
these abstractions look like? How can we discover a statistical path that naturally leads
to these different spectral methods? One such promising theory is discussed in this paper,
3

taking inspiration from celebrated history on nonparametric spectral analysis of time series.
The prescribed modern viewpoint shows an exciting confluence of nonparametric function
estimation, quantile domain methods, and Fourier harmonic analysis to unveil a more simple
conceptual structure of this subject, avoiding undue complexities.

1.4

Unified Construction Principle

We start by giving a glance to the unified spectral representation scheme, whose (nonparametric) statistical interpretation and justification is deferred until the next section, as it
requires a more technical background. It turns out that all known spectral graph techniques
are just different manifestations of this single general algorithm in some way.
Probability Notations and Definitions. We will stick to the following probabilistic notations.
P
Define network probability mass function by P (x, y; G) = A(x, y; G)/ x,y A(x, y; G); Vertex
P
P
probability mass function by p(x; G) = y P (x, y; G) and p(y; G) = x P (x, y; G) with the
associated quantile functions Q(u; X, G) and Q(v; Y, G). Finally, define the important Graph
Interaction Function by GIF(x, y; G) = P (x, y; G)/p(x; G)p(y; G), x, y ∈ V (G).
Algorithm 1. Nonparametric Spectral Graph Analysis: A Unified Algorithm
Step 1. For given discrete graph G of size n, construct GraField kernel function C :
[0, 1]2 → R+ ∪ {0} defined a.e by



P Q(u; X), Q(v; Y ); G

,
(1.1)
C (u, v; G) = GIF Q(u; X), Q(v; Y ); G =
p Q(u; X) p Q(v; Y )
where u = F (x; G), v = F (y; G) for x, y ∈ {1, 2, . . . , n}.
Step 2. Let {ξj }j≥0 be a complete orthonormal system on [0, 1]. Construct discrete transform basis ξj (x; F (X; G)) := ξj (F (x; G)), evaluated at the vertex-rank-transforms. They
satisfy the following orthonormality conditions (degree-weighted):
X
ξj (x; F (X)) ξk (x; F (X)) p(x; G) = δjk ,
x

thus constitutes an orthonormal basis for L2 (F ; G).
Step 3. Transform coding of graphs. Construct generalized spectral graph matrix M(G, ξ) ∈
Rn×n with respect to ηj (u) := ξj (Q(u; X); F ):
D

Z1

M[j, k; G, ξ] = ηj ,

(C − 1)ηk

E
L2 [0,1]

0

=

X

ξj (`; F )ξk (m; F )P (`, m; G).

(1.2)

`,m

Mξ can be viewed as a transform coefficient matrix of the orthogonal series expansion of
C (u, v; G) with respect to the product bases {ηj ηk }1≤j,k≤n .
4

P
Step 4. Perform the singular value decomposition (SVD) of Mξ = U ΛU T = k uk µk uTk ,
where where uij are the elements of the singular vector of moment matrix U = (u1 , . . . , un ),
and Λ = diag(µ1 , . . . , µn ), µ1 ≥ · · · µn ≥ 0.
Step 5. Obtain approximate Karhunen-Loéve (KL) representation basis (which acts as a
graph Fourier basis) of the graph G by
φek =

n
X

ujk ξj , for k = 1, . . . , n − 1,

j=1

which can be directly used for subsequent signal processing on graphs.

Organization. The next section is devoted to the fundamental principles that underlie this
algorithm. The unifying power will be shown in Sections 3 and 4. By doing so, we will
introduce empirical and smooth spectral graph analysis techniques. We also address the
open problem of obtaining a formal interpretation of “spectral regularization”. A deep
connection between high-dimensional discrete data smoothing and spectral regularization
is discovered. This new perspective provides, for the first time, the theoretical motivation
and fundamental justification for using regularized spectral methods, which were previously
considered to be empirical guesswork-based ad hoc solutions. Important application towards
spatial graph regression is discussed in Section 5. In Section 6, we end with concluding
remarks. Additional applications are available as supplementary materials.

2
2.1

Fundamentals of Statistical Spectral Graph Analysis
Graph Correlation Density Field

Wiener’s generalized harmonic analysis formulation on spectral representation theory of
time series starts by defining the autocorrelation function (ACF) of a signal. In particular,
Wiener–Khintchine theorem asserts ACF and the spectral density are Fourier duals of each
other. Analogous to Wiener’s correlation technique, we describe a new promising starting
point, from which we can develop the whole spectral graph theory systematically by bringing
nonparametric function estimation and harmonic analysis perspectives. Both statistical and
probabilistic motivations will be given.
Definition 1. For given discrete graph G of size n, the piecewise-constant bivariate kernel
function C : [0, 1]2 → R+ ∪ {0} is defined almost everywhere through



p Q(u; X), Q(v; Y ); Gn

 , 0 < u, v < 1. (2.1)
C (u, v; Gn ) = GIF Q(u; X), Q(v; Y ); Gn =
p Q(u; X) p Q(v; Y )
Theorem 1. GraField defined in Eq. (2.1) is a positive piecewise-constant kernel satisfying
ZZ
ZZ
X
C (u, v; Gn ) du dv =
C (u, v; Gn ) du dv = 1,
[0,1]2

(i,j)∈{1,...,n}2 I
ij

5

where
(
Iij (u, v) =

1,
0,

if (u, v) ∈ (F (i; X), F (i + 1; X)] × (F (j; Y ), F (j + 1; Y )]
elsewhere.

Note 1. The bivariate step-like shape of the GraField kernel is governed by the (piecewiseconstant left continuous) quantile functions Q(u; X, Gn ) and Q(v; Y, Gn ) of the discrete vertex
probability measures. As a result, in the continuum limit (let (Gn )n≥0 be a sequence of graphs
whose number of vertices tends to infinity n → ∞), the shape of the piecewise-constant discrete Cn approaches a “continuous field” over unit interval – a self-adjoint compact operator
on square integrable functions (defined on the graph) with respect to vertex probability
measure p(x; G) (see Section 2.2 for more details).
Motivation 1. We start with a diffusion based probabilistic interpretation of the GraField
kernel. The crucial point to note is that the “slices” of the C (2.1) can be expressed
as p(y|x; G)/p(y; G) in the vertex domain. This alternative conditional probability-based
viewpoint suggests a connection with the random walk on the graph. Interpret p(y|x; G) as
transition probability from vertex x to vertex y in one time step. Also note that p(y; G) is
the stationary probability distribution on the graph, as we have limt→∞ p(t, y|x; G) = p(y; G)
regardless of the initial starting point x (moreover, for connected graphs the stationary
distribution is unique). Here p(t, y|x; G) denotes the probability distribution of a random
walk landing at location y at time t, starting at the vertex x. See Lovász (1993) for an
excellent survey on the theory of random walks on graphs.
The graph affinity function measures how the transition probability p(y|x; G) is different from
the “baseline” stationary distribution (long-run stable behavior) p(y; G). That comparison
ratio is the fundamental interaction function for graphs which we denote by GIF(x, y; G).
This probabilistic interpretation along with Theorem 2 and 5 will allow us to integrate the
diffusion map (Coifman and Lafon, 2006) technique into our general statistical framework
in Section 3.2.
Motivation 2. GraField compactly represents the affinity or strength of ties (or interactions) between every pair of vertices in the graph. To make this clear, let us consider the
following adjacency matrix of a social network representing 4 employees of an organization


A=

0
2
0
0

2
0
3
3

0
3
0
3

0
3
,
3
0

where the weights reflect numbers of communication (say email messages or coappearances in social events etc.). Our interest lies in understanding the strength of association between the employees i.e., Strength(x, y) for all pairs of vertices. Looking at the
matrix A (or equivalently based on the histogram network estimator p(x, y; G) = A/N
P
with N =
x,y A(x, y) = 22) one might be tempted to conclude that the link between
6

employee 1 and 2 is the weakest, as they have communicated only twice, whereas employees 2, 3, and 4 constitute strong ties, as they have interacted more frequently. Now, the
surprising fact is that (i) Strength(1, 2) is twice that of Strength(2, 3) and Strength(2, 4);
also (ii) Strength(1, 2) is 1.5 times of Strength(3, 4)! To understand this paradox, compute the vertex-domain empirical GraField kernel matrix (Definition 1) with (x, y)th entry
N ·A(x, y; G)/d(x)d(y)


0

22/8

0

0

 22/8
0
22/16 22/16 
.
Cn = 
 0 22/16
0
22/12 
0

22/16 22/12

0

This toy example is in fact a small portion (with members 1, 9, 31 and 33) of the famous
Zachary’s karate club data, where the first two members were from Mr. Hi’s group and the
remaining two were from John’s group. The purpose of this illustrative example is not to
completely dismiss the adjacency or empirical graphon (Lovász and Szegedy, 2006) based
analysis but to caution the practitioners so as not to confuse the terminology “strength of
association” with “weights” of the adjacency matrix – the two are very different objects.
Existing literature uses them interchangeably without paying much attention. The crux of
the matter is: association does not depend on the raw edge-density, it is a “comparison
edge-density” that is captured by the GraField; see Section 3.2 for its intriguing connection
with diffusion based graph-distance.
Motivation 3. GraField can also be viewed as properly “renormalized Graphon,” which
is reminiscent of Wassily Hoeffding’s “standardized distributions” idea (Hoeffding, 1940).
Thus, it can be interpreted as a discrete analogue of copula (the Latin word copula means
“a link, tie, bond”) density for random graphs that captures the underlying correlation field.
We study the structure of graphs in the spectral domain via this fundamental graph kernel
C that characterizes the implicit connectedness or tie-strength between pairs of vertices.
Fourier-type spectral expansion results of the density matrix C are discussed in the ensuing
section, which is at the heart of our approach. We will demonstrate that this correlation
density operator-based formalism provides a useful perspective for spectral analysis of graphs
that allows unification.

2.2

Karhunen-Loéve Representation of Graph

We define the Karhunen-Loéve (KL) representation of a graph G based on the spectral
expansion of its GraField function C (u, v; G). Schmidt decomposition (Schmidt, 1907) of
C yields the following spectral representation theorem of the graph.
Theorem 2. The square integrable graph correlation density kernel C : [0, 1]2 → R+ ∪ {0}
of two-variables admits the following canonical representation
C (u, v; Gn ) = 1 +

n−1
X
k=1

7

λk φk (u)φk (v),

(2.2)

where the non-negative λ1 ≥ λ2 ≥ · · · λn−1 ≥ 0 are singular values and {φk }k≥1 are the
orthonormal singular functions hφj , φk iL 2 [0,1] = δjk , for j, k = 1, . . . , n − 1, which can be
evaluated as the solution of the following integral equation relation
Z
[C (u, v; G) − 1] φk (v) dv = λk φk (u),
k = 1, 2, . . . , n − 1.
(2.3)
[0,1]

Remark 1. By virtue of the properties of Karhunen-Loéve (KL) expansion (Loève, 1955),
the eigenfunction basis φk satisfying (2.3) provides the optimal low-rank representation of a
graph in the mean square error sense. In other words, {φk } bases capture the graph topology in the smallest embedding dimension and thus carries practical significance for graph
compression. Hence, we can call those functions the optimal coordinate functions or Fourier
representation bases. Accordingly, the fundamental statistical modeling problem hinges on
finding approximate solutions to the optimal graph coordinate system {φ1 , . . . , φn−1 } satisfying the integral equation (2.3).
Definition 2. Any function or signal y ∈ Rn defined on the vertices of the graph y : V 7→ R
P
such that kyk2 = x∈V (G) |y(x)|2 p(x; G) < ∞, can be represented as a linear combination of
the Schmidt bases of the GraField density matrix C . Define the generalized graph Fourier
transform of y
n
X
yb(λk ) := hy, φk i =
y(x)φk [F (x; G)].
x=1

This spectral or frequency domain representation of a signal, belonging to the square integrable Hilbert space L 2 (F, G) equipped with the inner product
hy, ziL 2 (F ;G) =

X

y(x)z(x)p(x; G),

x∈V (G)

allows us to construct efficient graph learning algorithms. As {φk }’s are KL spectral bases,
the vector of projections onto this basis function decay rapidly, and hence may be truncated
aggressively to capture the structure in a small number of bits.
Definition 3. The entropy (or energy) of a discrete graph G, is defined using the Parseval
relation of the canonical representation
ZZ
X
Entropy(G) =
(C − 1)2 du dv =
|λk |2 .
k

[0,1]2

This quantity, which captures the departure of uniformity of the C , can be interpreted as
a measure of ‘structure’ or the ‘compressibility’ of the graph. Entropy measure can also be
used to (i) define graph homogeneity and (ii) design fast algorithms for graph isomorphism.
For homogeneous graphs, the shape of the correlation density field is uniform over the unit
square. The power at each harmonic components, as a function of frequency, is called
the power spectrum of the graph. Flat spectrum is equivalent to the analogous notion of
white-noise process.
8

2.3

Nonparametric Spectral Approximation Theory

We view the spectral graph learning algorithm as a method of approximating (λk , φk )k≥1 that
satisfies the integral equation (2.3), corresponding to the graph kernel C (u, v; G). In practice,
often the most important features of a graph can be well characterized and approximated by
a few top (dominating) singular-pairs. The statistical estimation problem can be summarized
as follows:
n

o
b1 , φb1 , . . . , λ
bn−1 , φbn−1
An×n 7→ C 7→
λ
that satisfies Eq. (2.3).
2.3.1

From Continuous to Discrete Basis and Back

The Fourier-type nonparametric spectral approximation method starts by choosing an expansion basis. Let ξj denote an orthonormal basis of L 2 [0, 1]. Construct the discrete transform basis of Rn by evaluating ξj at the vertex-rank-transforms ξj (F (x; G)) := ξj (x; F, G).
Verify that ξj (x; F, G) are orthogonal with respect to the measure p(x; G) satisfying
X
ξj (x; F, G) ξk (x; F, G) p(x; G) = 0, for j 6= k.
x

Thus form an orthogonal basis of the Hilbert space L 2 (F ; G). Define η-functions (quantile
domain unit bases), generated from mother ξj by


ηj (u; G) = ξj Q(u; X); F, G , 0 < u < 1,
as piecewise-constant (left-continuous) functions over the irregular grid {0, p(1), p(1)+p(2), . . . ,
Pn
j=1 p(j) = 1} satisfying hηj , ηk iL 2 [0,1] = 0, if j 6= k. They will provide a useful tool to recast conventional matrix calculus-based approaches as a functional statistical problem. One
unique aspect of our construction is that in the continuum limit (as the size of the graph
n → ∞) the discrete η-basis of Rn approaches the mother ξ-function, a basis of L 2 [0, 1].
We call it asymptotic “reproducing” property.
2.3.2

Projection Methods for Eigenvector Approximation

We are interested in the nonparametric estimation of eigenpairs {λk , φk }k≥0 . Approximate
the unknown eigenvectors by the projection, Pn φk , on the span{ηj , j = 1, . . . , n} defined by
φk (u) ≈ Pn φk =

n
X

θjk ηj (u),

0<u<1

(2.4)

j=1

where θjk are the unknown coefficients to be estimated.
Definition 4 (Orthogonal Discrete Graph Transform). We introduce a generalized concept
of matrices associated with graphs called the G-matrix. Define discrete graph transform
with respect to an orthonormal system η as
M[j, k; η, G] =

D

Z1
(C − 1)ηk

ηj ,
0

9

E
L 2 [0,1]

for j, k = 1, . . . , n.

(2.5)

Equivalently, we can define the discrete graph transform to be the coefficient matrix of the
orthogonal series expansion of the GraField kernel C (u, v; G) with respect to the product
bases {ηj ηk }1≤j,k≤n . As a practical significance, this generalization provides a systematic
recipe for converting the graph problem into a “suitable” matrix problem:
Gn (V, E)

−→

An×n

C (u, v; Gn )

−→

{η1 ,...,ηn }

−−−−−→
Eq. (2.5)

M(η, G) ∈ Rn×n .

Theorem 3. The G-matrix (2.5) can also be interpreted as a “covariance” operator for a
discrete graph by recognizing the following equivalent representation for x, y ∈ {1, 2, . . . . , n}
X





M[j, k; ξ, G] = EP ξj F (X; G) ξk F (Y ; G)
=
P (x, y; G) ξj F (x; G) ξj F (y; G) ,
x,y

This can be proved using the basic quantile mechanics fact that Q(F (X)) = X holds with
probability 1 (see Parzen (1979)). Next, we present a general approximation scheme that
provides an effective method of discrete graph analysis in the frequency domain.
Theorem 4 (Nonparametric spectral approximation). The Fourier coefficients {θjk } of the
projection estimators (2.4) of the GraField eigenfunctions (eigenvalues and eigenvectors),
satisfying the integral equation (2.3), can be obtained by solving the following generalized
matrix eigenvalue problem
MΘ = SΘ∆,
(2.6)
R1
where Mjk = ηj , 0 (Cn − 1)ηk L 2 [0,1] , and Sjk = ηj , ηj L 2 [0,1] .
To prove define the residual of the governing equation (2.3) by expanding φk as series
expansion (2.4),
R(u) ≡

X
j

θjk

h Z1

i

C (u, v; G) − 1 ηj (v) dv − λk ηj (u) = 0.

(2.7)

0

Now for complete and orthonormal {ηj } requiring the error R(u) to be zero is equivalent to
the statement that R(u) is orthogonal to each of the basis functions
R(u), ηk (u)

L 2 [0,1]

= 0,

k = 1, . . . , n.

(2.8)

This leads to the following set of equations:
X
j

θjk

h ZZ

1
i
X Z


C (u, v; Gn )−1 ηj (v)ηk (u) dv du − λk
θjk
ηj (u)ηk (u) du = 0. (2.9)
j

[0,1]2

0

Theorem 4 plays a key role in our statistical reformulation. In particular, we will show how
the fundamental equation (2.6) provides the desired unity among different spectral graph
techniques by systematically constructing a “suitable” class of coordinate functions.

10

Note 2. The fundamental idea behind the Rietz-Galerkin (Galerkin, 1915) style approximation scheme for solving variational problems in Hilbert space played a pivotal inspiring
role to formalize the statistical basis of the proposed computational approach.
Note 3. Our nonparametric spectral approximation theory based on eigenanalysis of G-Matrix,
remains unchanged for any choice η-function (following the recipe for construction given in
Sec 2.3.1), which paves the way for the generalized harmonic analysis of graphs.
The next two sections investigate this general scheme under various choices of η-functions,
and nonparametric estimation methods of Cn . By doing so, many “mysterious similarities”
among different spectral graph algorithms are discovered which were not known before.

3

Empirical Spectral Graph Analysis

Three popular traditional spectral graph analysis models will be synthesized in this section.

3.1

Laplacian Spectral Analysis

Laplacian is probably the most heavily used spectral graph technique in practice. Here we
will demonstrate for the first time how the Laplacian of a graph naturally originates by
purely statistical reasoning, totally free from the classical combinatorial based logic.
Degree-Adaptive Block-pulse Basis Functions. One of the fundamental, yet universally valid (for any graph) choice for {ηj }1≤j≤n is the indicator top hat functions (also known
as block-pulse basis functions, or in short BPFs). However, instead of defining the BPFs
on a uniform grid (which is the usual practice) here (following Sec 2.3.1) we define them on
P
the non-uniform mesh 0 = u0 < u1 · · · < un = 1 over [0,1], where uj = x≤j p(x; X) with
local support
(
p−1/2 (j) for uj−1 < u ≤ uj ;
(3.1)
ηj (u) =
0
elsewhere.
They are disjoint, orthogonal, and a complete set of functions satisfying
Z1
0

Z1
Z1
p
2
ηj (u) du = 1, and
ηj (u)ηk (u) du = δjk .
ηj (u) du = p(j),
0

0

Note 4. The shape (amplitudes and block lengths) of our specially designed BPFs depend
on the specific graph structure via p(x; G) as shown in Fig 1. In order to obtain the spectral
domain representation of the graph, it is required to estimate the spectra of GraField kernel
φk , by representing them as block pulse series. The next result describes the required
computational scheme for estimating the unknown expansion coefficients {θjk }.
Theorem 5. Let φ1 , . . . , φn the canonical Schmidt bases of L 2 graph kernel C (u, v; G),
satisfying the integral equation (2.3). Then the empirical solution of (2.3) for block-pulse

11

2

2

1
3
1

4

3

0.0

0.0

0.6

1.0

1.2

2.0

4

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0
0.0
0.0

0.0

1.0

1.0

2.0

2.0

0.0

1.0

1.0

2.0

2.0

0.0

1.0

1.0

2.0

2.0

0.0

Figure 1: Two graphs and their corresponding degree-adaptive block-pulse functions. The
amplitudes and the block length of the indicator basis functions (3.1) depends on the degree
distribution of that graph.
orthogonal series approximated (3.1) Fourier coefficients {θjk } can equivalently be written
down in closed form as the following matrix eigen-value problem
L∗ [θ] = λθ,

(3.2)
1/2

where L∗ = L − uuT , L is the Laplacian matrix, u = Dp 1n , and Dp = diag(p1 , . . . , pn ).
Note that the discrete GraField kernel takes value GIF(j, k; Gn ) = P (j, k; Gn )/p(j; Gn )p(k; Gn )
and the tensor–product bases ηj (u)ηk (v) take value p−1/2 (j; Gn )p−1/2 (k; Gn ) over the rectangle Ijk (u, v) for 0 < u, v < 1. This observation reduces the master equation (2.9) to the
following system of linear algebraic equations expressed in the vertex domain:
#
"
X
p
p
p(j, k)
θjk p
− p(j) p(k) − λk δjk = 0.
(3.3)
p(j)p(k)
j
12

Empirical plugin nonparametric estimator. The estimating equation (3.3) contains
unknown network and vertex probability mass functions which need to be estimated from
the data. The most basic nonparametric estimates are Pe(j, k; Gn ) = A[j, k; Gn ]/N and
pe(j; Gn ) = d(j; Gn )/N . By plugging these empirical estimators into Eq. (3.3), said equation
can be rewritten as the following compact matrix form:
h
√ √ Ti
(3.4)
L − N −1 d d θ̃ = λ̃θ̃,
where L = D−1/2 AD−1/2 is the graph Laplacian matrix.
Significance 1. Theorem 5 allows us to interpret the graph Laplacian as the empirical
f Gn ) under degree-adaptive indicator basis choice for the η shape function.
G-matrix M(η,
Our technique provides a completely nonparametric statistical derivation of an algorithek versus k can be considered as a “raw
mic spectral graph analysis tool. The plot of λ
periodogram” analogue for graph data analysis.

3.2

Diffusion Map

We provide a statistical derivation of Coifman’s diffusion map (Coifman and Lafon, 2006)
algorithm, which hinges upon the following key result.
Theorem 6. The empirical GraField admits the following vertex-domain spectral diffusion
decomposition at any finite time t
X
pe(t, y|x; Gn )
λ̃tk φ̃k (x; F )φ̃k (y; F ),
= 1+
pe(y; Gn )
k

(3.5)

−1/2

where φek = Dp̃ uk , (uk is the kth eigenvector of the Laplacian matrix L), (φk ◦ F )(·) is
abbreviated as φk (·; F ), pe(y|x; G) = T (x, y), and T = D−1 A is the transition matrix of a
random walk on G with stationary distribution p̃(y; Gn ) = d(y; Gn )/N .
−1/2
Replacing the estimated coefficients from Theorem 5 into (2.4) yields φek = Dp uk , where
uk is the kth eigenvector of the Laplacian matrix L, immediately leads to the following
vertex-domain spectral decomposition result of the empirical GraField:

X
pe(y|x; G)
= 1+
λ̃k φ̃k (x)φ̃k (y),
pe(y; G)
k

(3.6)

For an alternative proof of the expansion (3.6) see Appendix section of Coifman and Lafon
(2006) by noting φek are the right eigenvectors of random walk Laplacian T .
Significance 2. In light of Theorem 6, define diffusion map coordinates at time t as the
mapping from x to the vector


x 7−→ λt1 φ1 (x; F ), . . . , λtk φk (x; F ) , x ∈ {1, 2, . . . , n},
13

which can be viewed as an approximate ‘optimal’ Karhunen-Loéve representation basis and
thus can be used for non-linear embedding of graphs. Define diffusion distance, a measure
of similarity between two nodes of a graph, as the Euclidean distance in the diffusion map
space
X 
2
0
(3.7)
λ2t
Dt2 (x, x0 ) =
j φj (x; F ) − φj (x ; F ) .
j≥1

This procedure is known as diffusion map, which has been extremely successful tool for
manifold learning2 . Our approach provides an additional insight and justification for the
diffusion coordinates by interpreting it as the strength of connectivity profile for each vertex,
thus establishing a close connection with empirical GraField.

3.3

Modularity Spectral Analysis

Theorem 7. To approximate the KL graph basis φk =
u ≤ uj ) to be the characteristic function satisfying
Z1

Z1
ηj (u) du =

0

P

j

θjk ηj , choose ηj (u) = I(uj−1 <

ηj2 (u) du = p(j; G).

0

Then the corresponding empirically estimated spectral graph equation (2.9) can equivalently
be reduced to the following generalized eigenvalue equation in terms of the matrix B =
A − N −1 ddT
Bα = λDα.
(3.8)
Significance 3. The matrix B, known as modularity matrix, was introduced by Newman
(2006) from an entirely different motivation. Our analysis reveals that the Laplacian and
Modularity based spectral graph analyses are equivalent in the sense that they inherently use
the same underlying basis expansion (one is a rescaled version of the other) to approximate
the optimal graph bases.
Remark 2. Solutions (eigenfunctions) of the GraField estimating equation based on the
G-matrix under the proposed specialized nonparametric approximation scheme provides a
systematic and unified framework for spectral graph analysis. As an application of this
general formulation, we have shown how one can synthesize the well-known Laplacian, diffusion map, and modularity spectral algorithms and view them as “empirical” spectral graph
analysis methods. It is one of those rare occasions where one can witness the convergence
of statistical, algorithmic and geometry-motivated computational models.
2

Manifold learning: Data-driven learning of the “appropriate” coordinates to identify the intrinsic nonlinear structure of high-dimensional data. We claim the concept of GraField allows decoupling of the
geometrical aspect from the probability distribution on the manifold.

14

4

Smoothed Spectral Graph Analysis

Smoothness is a universal requirement for constructing credible nonparametric estimators.
Spectral graph analysis is also no exception. An improved smooth version of raw-empirical
spectral graph techniques will be discussed, revealing a simple and straightforward statistical
explanation of the origin of regularized Laplacian techniques.

4.1

High-dimensional Undersampled Regime

Recall from Theorem 4 that the generalized matrix eigenvalue equation (2.6) depends on the
unknown network and vertex probability mass functions. This leads us to the question of
estimating the unknown distribution P = (p1 , p2 , . . . , pn ) (support size = size of the graph
P
= n) based on N sample, where N = ni=1 di = 2|E|. Previously (Theorems 5-7) we have
used the unsmoothed maximum likelihood estimate (MLE) pe(x; G) to construct our empirical spectral approximation algorithms, which is the unique minimum variance unbiased
estimator. Under standard asymptotic setup, where the dimension of the parameter space n
is fixed and the sample size N tends to infinity, the law of large numbers ensures optimality
of pe. As a consequence, empirical spectral analysis techniques are expected to work quite
well for dense graphs.
Estimation of probabilities from sparse data. However, the raw pe is known to be strictly suboptimal (Witten and Bell, 1991) and unreliable in the high-dimensional sparse-regime where
N/n = O(1) (i.e., when parameter dimension and the sample size are comparably large).
This situation can easily arise for modern day large sparse graphs where the ratio N/n is
small and there are many nodes with low degree, as is the case of degree sparsity. The naive
MLE estimator can become unacceptably noisy (high variability) due to the huge size and
sparse nature of the distribution. In order to reduce the fluctuations of “spiky” empirical
estimates, some form of “smoothing” is necessary. The question remains: How to tackle
this high-dimensional discrete probability estimation problem, as this directly impacts the
quality of our nonparametric spectral approximation.
My main purpose in the next section is to describe one such promising technique for smoothing raw-empirical probability estimates, which is flexible and in principle can be applied to
any sparse data.

4.2

Spectral Smoothing

We seek a practical solution for circumventing this problem that lends itself to fast computation. The solution, that is both the simplest and remarkably serviceable, is the Laplace/Additive
smoothing (Laplace, 1951) and its variants, which excel in sparse regimes (Fienberg and
Holland, 1973, Witten and Bell, 1991). The MLE and Laplace estimates of the discrete
distribution p(j; Gn ) are respectively given by
Raw-empirical MLE estimates:

pe(j; Gn ) =
15

dj
;
N

Smooth Laplace estimates:

pbτ (j; Gn ) =

dj + τ
N + nτ

(j = 1, . . . , n).

Note that the smoothed distribution pbτ can be expressed as a convex combination of the
empirical distribution pe and the discrete uninform distribution 1/n
 
N
nτ
1
pbτ (j; Gn ) =
pe(j; Gn ) +
,
(4.1)
N + nτ
N + nτ n
which provides a Stein-type shrinkage estimator of the unknown probability mass function
p. The shrinkage significantly reduces the variance, at the expense of slightly increasing the
bias.
Choice of τ . The next issue is how to select the “flattening constant” τ . The following
choices of τ are most popular in the literature3 :

1
Laplace estimator;



 1/2
Krichevsky–Trofimov estimator;
τ=
 1/n
Perks estimator;


 √
N /n Minimax estimator (under L2 loss).
Note 5. Under increasing-dimension asymptotics, this class of estimator is often difficult
to improve without imposing additional smoothness constraints on the vertex probabilities;
see Bishop et al. (2007, Chapter 12). The latter may not be a valid assumption as nodes of
a graph offer no natural order in general.
With this understanding, smooth generalizations of empirical spectral graph techniques will
be discussed, which have a close connection with recently proposed spectral regularized
techniques.

4.3

Type-I Regularized Graph Laplacian

Construct τ -regularized smoothed empirical ηj;τ basis function by replacing the amplitude
−1/2
p−1/2 (j) by p̂τ (j) following (4.1). Incorporating this regularized trial basis, we have the
following modified G-matrix based linear algebraic estimating equation (2.6):
"
#
X
p
p
pe(j, k)
θjk p
(4.2)
− p̂τ (j) p̂τ (k) − λk δjk = 0.
p̂
(j)p̂
(k)
τ
τ
j
Theorem 8. The τ -regularized block-pulse series based spectral approximation scheme is
equivalent to representing or embedding discrete graphs in the continuous eigenspace of
Type-I Regularized Laplacian = Dτ−1/2 A Dτ−1/2 ,

(4.3)

where Dτ is a diagonal matrix with i-th entry di + τ .
Note 6. It is interesting to note that this exact regularized Laplacian formula was proposed
by Chaudhuri et al. (2012) and Qin and Rohe (2013), albeit from a very different motivation.
3

For more details on selection of τ see Fienberg and Holland (1973) and references therein.

16

4.4

Type-II Regularized Graph Laplacian

Theorem 9. Estimate the joint probability p(j, k; G) by extending the formula given in (4.1)
for the two-dimensional case as follows:
 
nτ
1
N
pe(j, k; G) +
,
(4.4)
pbτ (j, k; G) =
N + nτ
N + nτ n2
which is equivalent to replacing the original adjacency matrix by Aτ = A + (τ /n)11T . This
modification via smoothing in the estimating equation (4.2) leads to the following spectral
graph matrix
Type-II Regularized Laplacian = Dτ−1/2 Aτ Dτ−1/2 .

(4.5)

Note 7. Exactly the same form of regularization of Laplacian graph matrix (4.5) was proposed by Amini et al. (2013) as a fine-tuned empirical solution.
Significance 4. Recently there has been a lot of discussion on how to choose the spectral
regularization parameter τ . Using Davis–Kahan’s theorem, Joseph and Yu (2016) have
proposed a data-driven technique called DKest. Their method involves repeated eigendecomposition of regularized Laplacian matrix (which is a dense matrix) over a grid of
values of τ . There are two major concerns for DKest: heavy computational burden and
theoretical validity (which only holds for stochastic block models and its extensions). In
contrast, our analysis provides a strikingly simple recommendation without any additional
computational overhead (by connecting it to large sparse distribution smoothing), and thus
is much easier to construct, implement, and use in practice. Similar to Joseph and Yu
(2016), the critical parameter of our theory is the ratio N/n–the average degree of the
nodes, contrary to the previous results (Chaudhuri et al., 2012, Qin and Rohe, 2013) that
relied on the minimum degree.

4.5

Google’s PageRank Method

Smoothing of network and vertex probability distributions appearing in generalized matrix
equation (2.6) resulted in Type-I and Type-II regularized Laplacian methods. The third
possibility is to directly smooth the conditional or transitional probability matrix to develop a regularized version of random walk Laplacian (which we call Type-III regularized
Laplacian) method for large sparse graphs.
Smoothing conditional probability function. Consider a random walk on G with transition
probability T (i, j; G) = Pr(Xt+1 = j|Xt = i) ≥ 0. Note that the smoothing (4.4) can
equivalently be represented as
 
A(i, j; G) + τ /n
1
Tτ (i, j; G) =
= (1 − ατ )T (i, j; G) + ατ
,
(4.6)
di + τ
n

17

where the degree-adaptive regularization parameter ατ = τ (di + τ )−1 . One can construct
an equivalent and more simplified (non-adaptive) estimator by directly smoothing each row
of T (as it is a row-stochastic or row-Markov matrix) via Stein-like shrinkage


1/n 1/n . . . 1/n
 .
..
.. 
..
Tα = (1 − α)T + α  ..
(4.7)
.
.
.  , 0 < α < 1.
1/n 1/n . . . 1/n
Spectral analysis can now proceed on this regularized transition probability matrix (either
(4.6) or (4.7)) by substituting empirical transition matrix Te = D−1 A.
Significance 5. The non-adaptive regularized random-walk Laplacian (transition matrix)
Tα was introduced by Larry Page and Sergey Brin in 1996 (Brin and Page, 1999) and is
called the Google’s PageRank matrix. What seems natural enough from our nonparametric
smoothing perspective, is in fact known to be highly surprising adjustment–the “Teleportation trick”.
‘Teleporting is the essential distinguishing feature of the PageRank random walk that had
not appeared in the literature before’ – Vigna (2005).
Note 8. The Google’s PageRank matrix (which is different from the empirical random
walk Laplacian in an important way) is probably the most famous, earliest, and spectacular
example of spectral regularization that was originally introduced to counter the problem of
dangling nodes (nodes with no outgoing edges) by “connecting” the graph4 .
Web is a huge heterogeneous sparse graph where dangling nodes (dead-ends) and disconnected components are quite common. As a result random walk can get stuck, and may
cycle around an isolated set of pages. Smoothing allows the random walk to teleport to a
web page uniformly at random (or adaptively based on degree-weights) whenever it hits a
dead end. The steady state vector describes the long term visit rate–the PageRank score,
computed via eigen-decomposition of Tα .

4.6

Other Generalizations

The beauty of our statistical argument is that it immediately opens up several possibilities
to construct new types of spectral regularization schemes, which are otherwise hard to guess
using previous understanding. Two such promising techniques are discussed here. The first
one deals with Stein smoothing with data-driven shrinkage parameter.
Theorem 10. Under the sparse asymptotic setup where the domain size n → ∞ at the same
rate as the sample size N , the risk (expected squared-error loss) of smooth add-b
τ probability
estimate pbτb (4.1) with the following data-driven choice
P
N 2 − ni=1 d2i
(4.8)
τb = Pn 2
n i=1 di − N 2
4

From a statistical estimation viewpoint, this can be thought of as a way to escape from the “zerofrequency problem” for discrete probability estimation.

18

is uniformly smaller than the risk of the unsmoothed MLE estimator pe = N −1 (d1 , . . . , dn ).
In addition, the risk of pbτb is uniformly smaller than the risk of the estimator formed by
choosing τ = 1/2.
For proof one can use the technique described in Fienberg and Holland (1973). As a consequence of this result, we can further improve the spectral regularization algorithm by
selecting τ in a data-driven way, with no extra computation.
While the Laplace or additive smoothing performs well in general, there are situations where
they perform poorly (Gale and Church, 1994). The Good-Turing estimator (Good, 1953) is
often the next best choice which is given by
pbGT (i; G) =

$di +1 di + 1
,
·
$di
N

(4.9)

where $k denotes the number of nodes with degree k. An excellent discussion on this topic
can be found in Orlitsky et al. (2003). One can plug in the estimate (4.9) into the equation
(4.2) to generate new spectral graph regularization technique.
Significance 6. We wish to emphasize that our novelty lies in addressing the open problem
of obtaining a rigorous interpretation and extension of spectral regularization. To the best of
our knowledge this is the first work that provides a more formal and intuitive understanding
of the origin of spectral regularization. We have shown how the regularization naturally
arises as a consequence of high-dimensional discrete data smoothing. In addition, this
point of view allows us to select appropriate regularization parameter τ with no additional
computation.

5

Application to Graph Regression

We study the problem of graph regression as an interesting application of the proposed
nonparametric spectral analysis algorithm. Unlike traditional regression settings, here, one
is given n observations of the response and predictor variables over the graph. The goal is
to estimate the regression function by properly taking into account the underlying graphstructured information along with the set of covariates.
Meuse Data Modeling. We apply our frequency domain graph analysis to address the
spatial prediction problem. Fig 2 describes the Meuse data set, a well known geostatistical
dataset. There is a considerable spatial pattern one can see from Fig 2. We seek to estimate
a smooth regression function of the dependent variable Y (zinc concentration in the soil)
via generalized spectral regression that can exploit this spatial dependency. The graph was
formed according to the geographic distance between points based on the spatial locations of
the observations. We convert spatial data into signal supported on the graph by connecting
two vertices if the distance between two stations is smaller than a given coverage radius. The
maximum of the first nearest neighbor distances is used as a coverage radius to ensure at least
19

Zinc concentrations (ppm)

Flooding frequency classes

Distance to river

Soil type classes
1000

●
● ●

1800

● ● ● ●
●
● ● ●●
●
●
●
● ●
●
●
●
●
●
●●
●
● ● ●
● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●
● ● ●
●
● ●
● ●
●
●
●
●
● ●
●
● ● ●
●
● ●
●
●
●
●
●● ●
●
●
●
●
●
●
● ● ●
●
●
●
●
●
● ●
● ● ●
●
●
●
●
●
●
●
●● ●
●
●● ●
● ● ●
● ●
●● ●
● ●
●
●●
●
●●
●
●
●
●
●
● ●
●
●
●
● ●
●●
●
●
●
● ●
●
●
●
●

●
●

●
●

●
● ●

●

1600

1400

1200

1000

800

600

400

● ● ● ●
●
● ● ●●
●
●
●
● ●
●
●
●
●
●
●●
●
● ● ●
● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●
● ● ●
●
● ●
● ●
●
●
●
●
● ●
●
● ● ●
●
● ●
●
●
●
●
●● ●
●
●
●
●
●
●
● ● ●
●
●
●
●
●
● ●
● ● ●
●
●
●
●
●
●
●
●● ●
●
●● ●
● ● ●
● ●
●● ●
● ●
●
●●
●
●●
●
●
●
●
●
● ●
●
●
●
● ●
●●
●
●
●
● ●
●
●
●

●
●

200

●

●
●

●

●
● ●

●

●
●
●

1
2
3

● ● ● ●
●
● ● ●●
●
●
●
● ●
●
●
●
●
●
●●
●
● ● ●
● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●
● ● ●
●
● ●
● ●
●
●
●
●
● ●
●
● ● ●
●
● ●
●
●
●
●
●● ●
●
●
●
●
●
●
● ● ●
●
●
●
●
●
● ●
● ● ●
●
●
●
●
●
●
●
●● ●
●
●● ●
● ● ●
● ●
●● ●
● ●
●
●●
●
●●
●
●
●
●
●
● ●
●
●
●
● ●
●●
●
●
●
● ●
●
●
●

●

●

●

●

●

●

●

●
● ●

●

800

600

400

200

●
●

●

● ● ● ●
●
● ● ●●
●
●
●
● ●
●
●
●
●
●
●●
●
● ● ●
● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●
● ● ●
●
● ●
● ●
●
●
●
●
● ●
●
● ● ●
●
● ●
●
●
●
●
●● ●
●
●
●
●
●
●
● ● ●
●
●
●
●
●
● ●
● ● ●
●
●
●
●
●
●
●
●● ●
●
●● ●
● ● ●
● ●
●● ●
● ●
●
●●
●
●●
●
●
●
●
●
● ●
●
●
●
● ●
●●
●
●
●
● ●
●
●
●
●

●
●

●

●

●
●
●

●
●

●

Figure 2: The Meuse dataset consists of n = 155 observations taken on a support of 15 × 15
m from the top 0−20 cm of alluvial soils in a 5×2 km part of the right bank of the floodplain
of the the river Meuse, near Stein in Limburg Province (NL). The dependent variable Y
is the zinc concentration in the soil (in mg kg−1 ), shown in the leftmost figure. The other
three variables, flooding frequency class (1 = once in two years; 2 = once in ten years; 3 =
one in 50 years), distance to river Meuse (in metres), and soil type (1= light sandy clay; 2
= heavy sandy clay; 3 =silty light clay), are explanatory variables.
one neighbor for each node. Fig 3(C) shows the sizes of neighbours for each node ranging
from 1 to 22. Three nodes (# 82, 148, and 155) have only 1 neighbour; additionally one
can see a very weakly connected small cluster of three nodes, which is completely detached
from the bulk of the other nodes. The reason behind this heterogeneous degree distribution
(as shown in Fig 3) is the irregular spatial pattern of the Meuse data.
We model the relationship between Y and spatial graph topology by incorporating nonparametrically learned spectral representation basis. Thus, we expand Y in the eigenbasis of C
and the covariates for the purpose of smoothing, which effortlessly integrates the tools from
harmonic analysis on graphs and conventional regression analysis. The model (5.1) described
in the following algorithm simultaneously promotes spatial smoothness and sparsity.
Algorithm 2. Nonparametric Spectral Graph Regression
Step 1. Input: We observe {yi ; xi1 , . . . , xip } at vertex i of the graph G = (V, A) with size
|V | = n. The regularization parameters τ .
−1/2

Step 2. Construct τ -regularized block-pulse shape functions ηj;τ = p̂j;τ I(uj−1 < u ≤ uj )
for j = 1, . . . , n.
Step 3. Compute the smoothed G-Matrix with respect to trial bases {ηk;τ }1≤k≤n
Z1
D
E
Mτ [j, k; η, G] = ηj;τ , (Cn − 1)ηk;τ
for j, k = 1, . . . , n.
L2 [0,1]

0

20

1
2
3

●
● ●● ●
●
●●
●
● ● ●●
●
●
●
● ●●
●
● ●
●●
●●
●●
●●
●
●
●
● ●●
●
●
●
●
●
●
●
●
●
●
●
●●●●●
●
● ● ●
●
●●●
●
●●
●
●
●
●
●
●●
●
●
●
●
●●
●
●● ●●
●
●●
●
●
●●● ●
●
●
●
●
●
●
●
●●● ●
● ●●●
●●●
●●
●●
●●●●● ● ● ● ●●
●
●
●
● ●● ●● ●●●
●
● ● ● ●●●
●
●
●● ●● ●
●
●● ● ● ●
●
●
102

101

69
103

108
137

94
154

111
68

136
112
113
104
67131
85
114
135
86 115
107
65
116
64
63

97

142 143

96

100
98 99
74 76138 140
71
77 91 152
70 72 87 88 75 139
81
73
80
89
153
79
78
90

83

66

95

144

141 145

151

149

150

146

15

110

147

93

92

148

105

106

51

61

Frequency

62
117

50

12260 121

49
128 134
11948 127
132
126
47 44
125

12358
55

54

59

53

46

124
43 42
52 45
40
41

20

13
16
14 8
17
9
129
18
15
84
130
12
133
11 10
29
24
28 27
26
30
25

19

21
22
23

33 32

155

1

2

39
38
37
36

35 34

3

7

4

5

0

57

56

5

120

10

109

6

5

31

118

10

15

20

degree of the spatial Meuse data graph

82

Figure 3: From spatial to graph representation of the Meuse river data set. For each node
we observe {yi ; x1i , x2i , x3i }ni=1 . We address the problem of approximating the regression
function by simultaneously incorporating the effect of explanatory variables X and the
underlying spatial graph dependency.
Step 4. Construct smooth graph Fourier basis
φbk;τ (u) =

n
X

ujk ηj;τ , for k = 1, . . . , n − 1,

j=1

P
T
where M = U ΛU T =
k uk µk uk , and uij are the elements of the singular vector U =


(u1 , . . . , un ). Construct spectral basis matrix Φ = φb1;τ , . . . , φbk;τ ∈ Rn×k for the graph G.


Step 5. Construct combined graph regressor matrix XG = Φ; X , where X = [X1 , . . . Xp ] ∈
Rn×p is the matrix of predictor variables.
Step 6. Solve for βG = (βΦ , βX )T
βbG = argmin ky − XG βG k22 + λkβG k1 .

(5.1)

βG ∈Rk+p

Algorithm 2 extends traditional regression to data sets represented over a graph. The proposed frequency domain smoothing algorithm efficiently captures the spatial graph topology
via the spectral coefficients βbΦ ∈ Rk and can be interpreted as covariate-adjusted discrete
graph Fourier transform of the response variable Y . The `1 sparsity penalty automatically
selects the coefficients with largest magnitudes thus provides compression.
The following table shows that incorporating the spatial correlation in the baseline unstructured regression model using spectral orthogonal basis functions (which are estimated from
the spatial graph with k = 25) boosts the model fitting from 62.78% to 80.47%, which is an
improvement of approximately 18%.
21

X

e0 + X
Φ

e Iτ =1 + X
Φ

e Iτ =.5 + X
Φ

e Iτ =.28 + X
Φ

e II
Φ
τ =1 + X

e II
Φ
τ =.5 + X

e II
Φ
τ =.28 + X

62.78

74.80

78.16

80.47

80.43

80.37

80.45

80.43

e 0, Φ
e I and Φ
e II are the graph-adaptive piecewise constant orthogonal basis functions
Here Φ
τ
τ
derived respectively from the ordinary, Type-I regularized, and Type-II regularized Laplacian matrix. Using other forms of regularization the performance does not change much,
and thus we do not show it here. Our spectral method can also be interpreted as a kernel
smoother, where the spatial dependency is captured by the discrete GraField C . We finally
conclude that extension from traditional regression to graph-structured spectral regression
significantly improves the model accuracy.

6

Concluding Remarks

What is the “common core” behind all existing spectral graph techniques? Despite half a
century of research, it remains one of the most formidable open issues, if not the core problem
in modern spectral graph theory. The need for a unified theory become particularly urgent
and pressing in the era of “algorithm deluge” in order to systematize the practice of graph
data analysis.
As an attempt to shed some light on this previously unasked question, I have introduced
a framework for statistical reasoning that appears to offer a complete and autonomous description (totally free from the classical combinatorial or algorithmic linear algebra-based
languages) with the intention of leading the reader to a broad understanding of how different
techniques interrelate. The prescribed approach brings a fresh perspective by appropriately
transforming the spectral problem into a nonparametric approximation and smoothing problem to arrive at new algorithms, as well as both practical and theoretical results. It was a
great surprise for me to be able to deduce existing techniques from some underlying basic
principles in a self-consistent way. At the end, there is always the hope that a comprehensive understanding gained through unification will inspire new tools, to attack some of the
important open problems in this area.
Acknowledgments: The author benefited greatly from many fruitful discussions and exchanges of ideas he had at the John W. Tukey 100th Birthday Celebration meeting, Princeton University September 18, 2015. The author is particularly grateful to Professor Ronald
Coifman for suggesting the potential connection between our formulation and Diffusion map
(random walk on graph). We also thank Isaac Pesenson for helpful questions and discussions.

References
Amini, A. A., Chen, A., Bickel, P. J., Levina, E. (2013). Pseudo-likelihood methods
for community detection in large sparse networks. The Annals of Statistics, 41 2097–2122.
22

Belkin, M. and Niyogi, P. (2001). Laplacian eigenmaps and spectral techniques for
embedding and clustering. NIPS, 14 585–591.
Bishop, Y. M., Fienberg, S. E. and Holland, P. W. (2007). Discrete multivariate
analysis: theory and practice. Springer Science & Business Media.
Brin, S. and Page, L. (1999). Reprint of: The anatomy of a large-scale hypertextual web
search engine. Computer networks, 56 3825–3833.
Chaudhuri, K., Graham, F. C. and Tsiatas, A. (2012). Spectral clustering of graphs
with general degrees in the extended planted partition model. Journal of Machine Learning Research Workshop and Conference Proceedings, 35 1–23.
Chung, F. R. (1997). Spectral graph theory, vol. 92. American Mathematical Soc.
Coifman, R. R. and Lafon, S. (2006). Diffusion maps. Applied and computational
harmonic analysis, 21 5–30.
Coifman, R. R., Lafon, S., Lee, A. B., Maggioni, M., Nadler, B., Warner, F.
and Zucker, S. W. (2005). Geometric diffusions as a tool for harmonic analysis and
structure definition of data: Diffusion maps. Proceedings of the National Academy of
Sciences of the United States of America, 102 7426–7431.
Davis, E. and Sethuraman, S. (2016). Consistency of modularity clustering on random
geometric graphs. arXiv preprint arXiv:1604.03993.
Fiedler, M. (1973). Algebraic connectivity of graphs. Czechoslovak Mathematical Journal,
23 298–305.
Fienberg, S. E. and Holland, P. W. (1973). Simultaneous estimation of multinomial
cell probabilities. Journal of the American Statistical Association, 68 683–691.
Gale, W. and Church, K. (1994). What is wrong with adding one. Corpus-based research
into language, 1 189–198.
Galerkin, B. (1915). Series development for some cases of equilibrium of plates and beams
(in Russian). Wjestnik Ingenerow Petrograd, 19 897–908.
Good, I. J. (1953). The population frequencies of species and the estimation of population
parameters. Biometrika, 40 237–264.
Hoeffding, W. (1940). Massstabinvariante korrelationstheorie. Schriften des Mathematischen Seminars und des Instituts für Angewandte Mathematik der Universität Berlin, 5
179–233.
Joseph, A. and Yu, B. (2016). Impact of regularization on spectral clustering. Ann.
Statist., 44 1765–1791.
23

Laplace, P. S. (1951). A philosophical essay on probabilities, translated from the 6th
French edition by Frederick Wilson Truscott and Frederick Lincoln Emory. Dover Publications (New York, 1951).
Le, C. M., Levina, E. and Vershynin, R. (2015). Sparse random graphs: regularization
and concentration of the laplacian. arXiv:1502.03049.
Loève, M. (1955). Probability Theory; Foundations, Random Sequences. New York: D.
Van Nostrand Company.
Lovász, L. (1993). Random walks on graphs. Combinatorics, Paul erdos is eighty, 2 1–46.
Lovász, L. and Szegedy, B. (2006). Limits of dense graph sequences. Journal of Combinatorial Theory, Series B, 96 933–957.
Newman, M. E. (2006). Finding community structure in networks using the eigenvectors
of matrices. Physical review E, 74 036104.
Orlitsky, A., Santhanam, N. P. and Zhang, J. (2003). Always good turing: Asymptotically optimal probability estimation. Science, 302 427–431.
Parzen, E. (1979). Nonparametric statistical data modeling (with discussion). Journal of
the American Statistical Association, 74 105–131.
Perks, W. (1947). Some observations on inverse probability including a new indifference
rule. Journal of the Institute of Actuaries, 73 285–334.
Qin, T. and Rohe, K. (2013). Regularized spectral clustering under the degree-corrected
stochastic blockmodel. In Advances in Neural Information Processing Systems. 3120–3128.
Rohe, K., Chatterjee, S. and Yu, B. (2011). Spectral clustering and the highdimensional stochastic blockmodel. The Annals of Statistics 1878–1915.
Sarkar, P., Bickel, P. J. (2015). Role of normalization in spectral clustering for stochastic blockmodels. The Annals of Statistics, 43 962–990.
Schmidt, E. (1907). Zur Theorie der linearen und nicht linearen Integralgleichungen Zweite
Abhandlung. Mathematische Annalen, 64 433–476.
Vigna, S. (2005). Trurank: taking pagerank to the limit. 14th international conference on
World Wide Web 976–977.
Von Luxburg, U., Belkin, M. and Bousquet, O. (2008). Consistency of spectral
clustering. The Annals of Statistics 555–586.
Witten, I. H. and Bell, T. (1991). The zero-frequency problem: estimating the probabilities of novel events in adaptive text compression. Information Theory, IEEE Transactions
on, 37 1085–1094.
24

Supporting Online Material
Unified Statistical Theory of Spectral Graph Analysis
Subhadeep Mukhopadhyay
Department of Statistical Science, Temple University
Philadelphia, Pennsylvania, 19122, U.S.A.
deep@temple.edu
To further demonstrate the potential application of smooth spectral graph algorithms, in
this supplementary section we discuss the community detection problem that seeks to divides nodes into into k groups (clusters), with larger proportion of edges inside the group
(homogeneous) and comparatively sparser connections between groups to understand the
large-scale structure of network. Discovering community structure is of great importance in
many fields such as LSI design, parallel computing, computer vision, social networks, and
image segmentation.

Graph Clustering
In mathematical terms, the goal is to recover the graph signals (class labels) y : V 7→
{1, 2, . . . , k} based on the connectivity pattern or the relationship between the nodes. Representation of graph in the spectral or frequency domain via the nonlinear mapping Φ :
G(V, E) 7→ Rm using discrete KL basis of density matrix C as the co-ordinate is the most
important learning step in the community detection. This automatically generates spectral
features {φ1i , . . . , φmi }1≤i≤n for each vertex that can be used for building the distance or
similarity matrix to apply k-means or hierarchical clustering methods. In our examples, we
will apply k-means algorithm in the spectral domain, which seeks to minimizing the withincluster sum of squares. In practice, often the most stable spectral clustering algorithms
determine k by spectral gap: k = argmaxj |λj − λj+1 | + 1.
Algorithm 3. Nonparametric Spectral Graph Partitioning
1. Input: The adjacency matrix A ∈ Rn×n . Number of clusters k. The regularization
parameter τ .
2. Estimate the top k − 1 spectral connectivity profile for each node {φe1i;τ , . . . , φe(k−1)i;τ }
using Algorithm 2. Store it in Φ ∈ Rn×k−1 .
3. Apply k-means clustering by treating each row of Φ as a point in Rk−1 .
4. Output: The cluster assignments of n vertices of the graph C1 , . . . , Ck .

S–1

In addition to graph partitioning, the spectral ensemble {λk , φk }1≤k≤m of C contain a wealth
e1 (G; η) for the choice
of information on the graph structure. For example, the quantity 1 − λ
of {ηk } to be normalized top hat basis (3.1), is referred to as the algebraic connectivity,
whose magnitude reflects how well connected the overall graph is. The kmeans clustering after spectral embedding {φej (G; η)}1≤j≤k−1 finds approximate solution to the NP-hard
combinatorial optimization problem based on the normalized cut (Shi and Malik, 2000)
by relaxing the discreteness constraints into one that is continuous (sometimes known as
spectral relaxation).
Data and Results. We investigate four well-studied real-world networks for community
structure detection based on 7 variants of spectral clustering methods.
Example A [Political Blog data, Adamic and Glance (2005)] The data, which contains
1222 nodes and 16, 714 edges, were collected over a period of two months preceding the U.S.
Presidential Election of 2004 to study how often the political blogs refer to one another. The
linking structure of the political blogosphere was constructed by identifying whether a URL
present on the page of one blog references another political blog (extracted from blogrolls).
Each blog was manually labeled as liberal or conservative by Adamic and Glance (2005),
which we take as ground truth. The goal is to discover the community structure based on
these blog citations, which will shed light on the polarization in political blogs.
Table 1 shows the result of applying the spectral graph clustering algorithm on this political
web-blog data. The un-regularized Laplacian performs very poorly, whereas as both typeI/II regularized versions give significantly better results. The misclassification error drops
from 47.95% to 4.7% because of regularization. To better understand why regularization
plays a vital role, consider the degree distribution of the web-blog network as shown in
the bottom panel of Figure 2. It clearly shows the presence of a large number of lowdegree nodes, which necessitates the smoothing of high-dimensional discrete probability
p1 , . . . , p1222 . Thus, we perform the kmeans clustering after projecting the graph in the
Euclidean space spanned by Laplace smooth KL spectral basis {φek;τ }. Regularized spectral
methods correctly identify two dense clusters: liberal and conservative blogs, which rarely
links to a blog of a different political leaning, as shown in the middle panel of Fig 5.
Example B [US College Football, Grivan and Newman (2002)] The American football
network (with 115 vertices, 615 edges) depicts the schedule of football games between NCAA
Division IA colleges during the regular season of Fall 2000. Each node represents a college
team (identified by their college names) in the division, and two teams are linked if they have
played each other that season. The teams were divided into 11 “conferences” containing
around 8 to 12 teams each, which formed actual communities. The teams in the same
conference played more often compared to the other conferences, as shown in the middle
panel of Fig 5. A team played on average 7 intra- and 4 inter-conference games in the season.
Inter-conference play is not uniformly distributed; teams that are geographically close to one
another but belong to different conferences are more likely to play one another than teams
S–2

Table 1: We report % of misclassification error. We compare following seven different
Laplacian variants. K denotes the number of communities.
Type-I Reg. Laplacian
τ = 1/2 τ =

√
N /n

Type-II Reg. Laplacian
τ =1

τ = 1/2 τ =

√

Data

K

Laplacian

τ =1

N /n

PolBlogs

2

47.95%

4.9%

4.8%

5.4%

4.8%

4.7%

5.4%

Football

11

11.3%

7.83%

6.96%

6.96%

6.96%

7.83%

7.83%

MexicoPol

2

17.14%

14.2%

14.2%

14.2%

14.2%

14.2%

14.2%

Adjnoun

2

13.4%

12.5%

12.5%

12.5%

12.5%

12.5%

12.5%

separated by large geographic distances. As the communities are well defined, the American
football network provides an excellent real-world benchmark for testing community detection
algorithms.
Table 1 shows the performance of spectral community detection algorithms to identify the 11
clusters in the American football network data. The regularization boosts the performance
√
by 3-4%. In particular, τ = 1/2 and N /n produces the best result for Type-I regularized
Laplacian, while τ = 1 exhibits the best performance for Type-II regularized Laplacian.
Example C [The Political Network in Mexico, Gil-Mendieta and Schmidt (1996)] The data
(with 35 vertices and 117 edges) represents the complex relationship between politicians in
Mexico (including presidents and their close associates). The edge between two politicians
indicates a significant tie, which can either be political, business, or friendship. A classification of the politicians according to their professional background (1 - military force, 2 civilians: they fought each other for power) is given. We use this information to compare
our 7 spectral community detection algorithms.
Although this is a “small” network, challenges arise from the fact that the two communities cannot be separated easily due to the presence of a substantial number of betweencommunity edges, as depicted in Figs 4 and 5. The degree-sparsity is also evident from Fig 5
(bottom-panel). Table 1 compares seven spectral graph clustering methods. Regularization
yields 3% fewer misclassified nodes. Both the type-I and II regularized Laplacian methods
for all the choices of τ produce the same result.
Example D [Word Adjacencies, Newman (2006)] This is a adjacency network (with 112
vertices and 425 edges) of common adjectives and nouns in the novel David Copperfield
by English 19th century writer Charles Dickens. The graph was constructed by Newman

S–3

●
5

13
24

●
1

●
2

14

30
22

●
3

●
6

●

18

4

●
7

●

21

9

●

23
19

10

●
8

●
12

17

27

34

31

28
25

35

33

20
15
26

32
16

●
11

29

Figure 4:

Mexican political network. Two different colors (golden and blue) denotes the

two communities (military and civilians).
(2006). Nodes represent the 60 most commonly occurring adjectives and nouns, and an edge
connects any two words that appear adjacent to one another at any point in the book. Eight
of the words never appear adjacent to any of the others and are excluded from the network,
leaving a total of 112 vertices. The goal is to identify which words are adjectives and nouns
from the given adjacency network.
Note that typically adjectives occur next to nouns in English. Although it is possible for
adjectives to occur next to other adjectives (e.g., “united nonparametric statistics”) or for
nouns to occur next to other nouns (e.g., “machine learning”), these juxtapositions are less
common. As expected, Fig 5 (middle panel) shows an approximately bipartite connection
pattern among the nouns and adjectives.
A degree-sparse (skewed) distribution is evident from the bottom right of Fig 5. We apply
the seven spectral methods and the result is shown in Table 1. The traditional Laplacian
yields a 13.4% misclassification error for this dataset. We get better performance (although
the margin is not that significant) after spectral regularization via Laplace smoothing.

References
Adamic, L. A. and Glance, N. (2005). The political blogosphere and the 2004 us election:
divided they blog. In Proceedings of the 3rd international workshop on Link discovery.
ACM, 36–43.
Gil-Mendieta, J. and Schmidt, S. (1996). The political network in mexico. Social
Networks, 18 355–381.
S–4

Frequency

20

30

10
8
6

Frequency

10

4

20

30

Frequency

40

50

40

12

60

500
400

50

100

150

200

250

degree of the web−blog network

300

350

0

0

2

10
0

Frequency

300
200
100
0
0

7

8

9

10

11

12

13

degree of the US college football network

5

10

15

degree of the political network in Mexico

0

10

20

30

40

50

degree of the word adjacencies network

Figure 5: The columns denote the 4 datasets corresponding to the political blog, US college
football, politicians in Mexico, and word adjacencies networks. The first two rows display the
un-ordered and ordered adjacency matrix, and the final row depicts the degree distributions.

Grivan, M. and Newman, M. E. J. (2002). Community structure in social and biological
networks. Proc. Natl. Acad. Sci 7821–7826.
Newman, M. E. (2006). Finding community structure in networks using the eigenvectors
of matrices. Physical review E, 74 036104.
Shi, J. and Malik, J. (2000). Normalized cuts and image segmentation. Pattern Analysis
and Machine Intelligence, IEEE Transactions on, 22 888–905.

S–5

