1

A New Spectral Clustering Algorithm

arXiv:1710.02756v1 [cs.LG] 7 Oct 2017

W.R. Casper1 and Balu Nadiga2
Abstract—We present a new clustering algorithm that is based on searching for natural gaps in the components of the
lowest energy eigenvectors of the Laplacian of a graph. In comparing the performance of the proposed method with a set
of other popular methods (KMEANS, spectral-KMEANS, and an agglomerative method) in the context of the
Lancichinetti-Fortunato-Radicchi (LFR) Benchmark for undirected weighted overlapping networks, we find that the new
method outperforms the other spectral methods considered in certain parameter regimes. Finally, in an application to
climate data involving one of the most important modes of interannual climate variability, the El Niño Southern Oscillation
phenomenon, we demonstrate the ability of the new algorithm to readily identify different flavors of the phenomenon.
Index Terms—Clustering, eigenvectors, machine learning

F

1

I NTRODUCTION

C

LUSTERING is an unsupervised learning technique used to identify natural subgroups of a
set of data wherein the members of a subgroup share
certain properties. On representing the set of data as a
graph, clustering methods group vertices of the graph
into clusters based on the edge structure of the graph
[e.g., see 23]. We note here that the edge structure
itself may be derived purely from interactions between
vertices1 , e.g., as in a social network, or based on
a measure of similarity between the vertices themselves, or some combination of the two. Needless to
mention, clustering methods find use in an extremely
wide range of disciplines (e.g., neuroscience [7], social networks [17], computer vision [8], and digital
forensic analysis [4]). Indeed most modes of climate
variability are typically first identified by applying
clustering methods to climatological data sets [e.g.,
14, 3, 29, 10, 26, and others].
Since approaches to clustering are extremely varied, we do not attempt even a brief overview of such
methods. For that the reader is referred to reviews
such as [23]. Instead, we shift focus directly to the new
algorithm after briefly considering common aspects of
spectral approaches. In this article, we propose a spectral clustering algorithm based on direct analysis of
the magnitudes of components of the eigenvectors (of
the Laplacian matrix; see next section) with smallest
eigenvalues. We call this algorithm spectral maximum
gap method (SP-MGM). We compare the performance
of this algorithm with popular clustering algorithms
1 Department

of Mathematics, Louisiana State University, Baton
Rouge LA. Email:wcasper1@lsu.edu
2 CCS-2, Los Alamos National Lab, Los Alamos, New Mexico.
Email:balu@lanl.gov
1. and in which context clustering is commonly referred to as
community-detection

including k-means (KMEAN), more typical spectral
clustering based on running kmeans on eigenvectors
(SP-KMEANS), and agglomerative clustering (AGG).
To determine the performance of these algorithms,
we rely on a series of LFR benchmark graphs [13].
Notably, it is seen that on a certain series of these
graphs with varying size, average degree, and weight
mixing parameter, our algorithm consistently outperform these other methods. Following our benchmark
analysis, we apply SP-MGM to a collection of monthly
averaged sea surface temperatures in the Niño 3.4
region of the Pacific Ocean. We demonstrate that our
method correctly and readily identifies the various
flavors of El Niño and La Niña events.

2

S PECTRAL C LUSTERING A LGORITHMS

Spectral clustering algorithms rely on using spectral
characteristics of the Laplacian of a (weighted) graph
to partition the vertices of the graph into natural
clusters. Here, the Laplacian can refer to the normalized Laplacian, the non-normalized Laplacian, or
even stranger entities like the p-Laplacian [2]. The
use of spectral characteristics to identify clusters in
a graph may be justified in numerous ways, such as
by considering random walks, minimal cuts, or block
matrix diagonalization. For a recent survey on spectral
clustering methods, the reader is referred to [18]; for
an introduction to spectral clustering itself, see, e.g.,
[31].
To present a brief and intuitive description of spectral clustering methods, we consider the eigenvectors
of the Laplacian of a disconnected graph. Let G be a
graph with adjacence matrix A and degree matrix D
(the diagonal matrix whose nonzero entries are the degrees of the vertices of G). Then the (non-normalized)

Laplacian of G is defined as L = D − A. The matrix
L is positive semidefinite, and the kernel has a basis
spanned by the indicator functions of the connected
components of G. In particular, if G is disconnected
then the connected components of G may be deduced
from the kernel of L.
Clusters in a graph G are intuitively regions in
G where the vertices are connected by edges with
strong weights, relative to the weights of edges going
between clusters. For this reason, the Laplacian L of
e of
the graph G should be very close to the Laplacian L
e formed from G by deleting
the disconnected graph G
intercluster edges. This means that the eigenvalues of
L will correspond approximately to the eigenvalues
e . Since L
e is positive semidefinite, the smallest
of L
e is zero. Therefore the eigenvectors of
eigenvalue of L
the smallest eigenvalues of L should be small perture . Thus by analyzing
bations of vectors in the kernel of L
the structure of the eigenvectors corresponding to the
smallest handful of eigenvalues of L, we should expect
to retrieve information about the clusters of G. It is
exactly this idea that all spectral clustering algorithms
exploit. Note that if G is connected, then the kernel
of L consists of vectors whose entries are all identical,
and therefore does not provide any information. For
this reason, the kernel of L is typically ignored in
spectral clustering algorithms.
To further clarify this idea, consider the example
of three cliques of size 4 joined together by single
bonds in general position, as in Figure 1a. The nonnormalized Laplacian of this graph in block matrix
form is


K + K0
Q
QT


,
L=
QT
K + K0
Q


Q
QT
K + K0

(a) A graph with three clusters formed by
cliques of size 4. Vertex color indicates inclusion in a particular clique, and all edges have
weight 1. The nodes are also indexed in the
same order as their appearance in the Laplacian matrix.

(b) A basis for the eigenspace of eigenvectors of
the clique cluster graph with smallest nonzero
eigenvalue λ ≈ 0.5051. The eigenspace is twodimensional, and the basis elements are colorcoded. The structure of these eigenvectors reflects cluster membership in that the values (for
any given eigenvector) at indices belonging to
the same cluster are similar.

K the Laplacian of a 4 × 4 clique, K 0 = QQT + QT Q,
and Q the 4×4 matrix whose only nonzero entry is −1
and occurs in the bottom left corner. The eigenvalues
of this matrix are 0 and 6 with multiplicity 1, 6−λ and
λ with multiplicity 2, and 4 with multiplicity 6, where
here λ ≈ 0.5051. A basis for the for the eigenspace of
the lowest nonzero eigenvalue λ is plotted in Figure
1b. The stucture of these eigenvectors successfully
reveals the presence of the three clusters.
In our algorithm, spectral clustering on a given
data set occurs in three parts. We first create a similarity matrix A describing the similarity or connectedness
between each of our pieces of data. Secondly, we
calculate the eigenvectors of the Laplacian of A, and
sort them in terms of the magnitudes of the associated eigenvalues. Lastly, we use some number of the
eigenvectors with the smallest nonzero eigenvalues
to calculate a natural partition of the vertices of the
graph.

Fig. 1: An illustration of the basis for spectral clustering.

A similarity matrix for a collection of n indexed
data points is an n × n symmetric matrix A whose
i, j ’th entry is a nonnegative value measuring the
similarity between the i and j ’th data values. It is assumed that larger values indicate more similar objects,
and the diagonal elements are typically taken to be
0. Similarity matrices can be created out of data in a
variety of ways. One popular method is to start with
the distance dij between i and j (under some metric),
and then define Aij = exp(−d2ij /2σ) for i 6= j , where
here σ is a parameter.
The (symmetric) normalized Laplacian of A is the
graph Laplacian of the unique weighted graph G(A)
2

whose weight matrix is A. In particular, it is defined
by
L = I − D−1/2 AD−1/2

Algorithm 1: SP-G1 clustering

where here D is a diagonal matrix whose entries are
the sums of each of the column vectors of A. The
matrix L is positive-semidefinite and will always have
at least one eigenvector ~v with eigenvalue 0, given
by ~v = D1/2 [1 . . . 1]T . However, assuming that the
graph G(A) associated to A is connected, this will be
the only eigenvector with eigenvalue 0, up to constant
multiples.
The smallest nonzero eigenvectors of L encode
clustering information for L, and it is these eigenvectors which we wish to parse in the final step.
The reason these eigenvectors contain clustering information can be made intuitive by the following
argument. Clusters consist of regions in the graph with
strong interconnectivity. Specifically, given a vertex
in a cluster we would expect that on average the
strength of its connection with other vertices in the
same cluster should be much greater than the strength
of its connection with vertices in other clusters. From
this point of view, we can decompose the adjacency
matrix A as A = B + C where B is a block diagonal
matrix whose blocks are formed by eliminating intercluster bonds, and where C = B −A is relatively small
compared to B . In this case, the Laplacian L of A may
be written as L = LB + L0 where LB is the Laplacian
of B and L0 is small relative to LB . Consequently
the eigenvectors of LB will be comparable ot the
eigenvectors of L. In particular the 0 eigenvectors of
LB correspond precisely with the eigenvectors of L
with small eigenvalues. Moreover, these eigenvectors
describe the connected components of the graph G(B)
of B , and these are precisely the clusters of G(A).
In this way, using the eigenvectors of the smallest
nonzero eigenvalues makes sense.
Where various clustering algorithms differ is in
this third step, where the eigenvector data is used to
calculate clusters. Before proposing our own clustering scheme, we will recount two common methods
found in the literature, the simplest method SP-G1 and
the most common method SP-KMEANS.

1

2.1

2
3
4
5
6
7
8
9

2.2

Input : weight matrix A
Output: two clusters
Calculate normalized Laplacian L of A
Calculate smallest positive eigenvalue λ1 of L
Calculate an eigenvector ~v of λ1
Choose a threshold r
if vi < r then
put node i in cluster 1
else
put node i in cluster 2
end

SP-KMEANS

A more complicated and very reliable clustering algorithm instead uses eigenvectors ~v1 , . . . , ~vk of the k
smallest positive eigenvalues 0 < λ1 ≤ · · · ≤ λk of
the Laplacian L of A. These eigenvectors are used as
the column vectors of an n × k matrix V . The matrix
V is then normalized so that each of the rows has
norm 1, and the associated row vectors are viewed
as n points on the surface of a sphere in Rk . Then
by running KMEANS on these points, we return k
clusters of the data C1 , . . . , Ck which partition the
original graph [19][16]. The reliability and familiarity
of this algorithm has lead to its popularity as one of
the most common forms of spectral clustering, with
implementations in both R and python [20][11].
Algorithm 2: SP-KMEANS clustering

6

Input : weight matrix A
number of clusters k
Output: k clusters
Calculate normalized Laplacian L of A
Calculate k smallest positive eigenvalues
λ1 < · · · < λk of L
Calculate an eigenvector ~vj of λi
Construct an n × k matrix V = [~v1 . . . ~vk ]
Normalize V so that each row vector has norm 1
Run KMEANS on the row vectors of V

2.3

Spectral Maximum Gap Method (SP-MGM)

1
2

3
4
5

SP-G1 clustering

The oldest and most basic algorithm is SP-G1, which
decomposes the data into two clusters based on the
entries of the eigenvector ~v = [v1 . . . vn ]T with the
lowest nonzero eigenvalue [6][9]. Specifically, a theshold value r needs to be chosen. Then the i’th data point
is put into the first cluster if vi < r, and otherwise
it is put into the second cluster. Of course, there
P are
several natural choices for r, including r = n−1 i vi .
Of course, more sophisticated methods for the choice
of r exist[25].

In the new SP-MGM algorithm we propose, the eigenvectors ~v1 , . . . , ~vk of the k smallest positive eigenvalues 0 < λ1 ≤ · · · ≤ λk of the Laplacian L of
A are first calculated. Next, for each j we create a
new vector ~
uj whose entries are the entries of ~vj in
increasing order. Let vji and uji denote the entries
of ~vj and ~
uj , respectively. Differences uj(i+1) − uji
between subsequent entries of ~
uj represent “gaps" in
~vj .
3

We choose `j so that ~
uj(`j +1) − ~uj`j is the largest
gap. Then we create a partition Pj , Qj of {0, . . . , n−1}
by placing vertex i in Pj if vji <= uj`j , and placing
vertex i in Qj otherwise. Finally, for each subset S of
U = {1, . . . , k} we define a cluster

! 
\
\
CS =
Pi ∩ 
Qi  ,
i∈S

permutations to the Laplacian of the graph G. We
provide full mathematical details in the appendix for
the interested reader. However, we restrict ourselves
to a brief description of the justification here so as to
keep the narrative more generally accessible: First, we
show mathematically that if G consists of several welldefined clusters C1 , . . . , Ck of size n1 , . . . , nk then

ρ1
1
< 1/2 ,
λi,2
4nm
for all 1 ≤ m ≤ k . Here λi,2 and ρ1 are spectral data which measure the quality of the clusters.
Specifically, λi,2 is the smallest nonzero eigenvalue
of cluster Ci and measures the connectivity of the
cluster Ci . For this reason, it is sometimes called the
algebraic connectivity of the cluster Ci . The number
ρ1 is the spectral radius of the Laplacian of the graph
G formed by deleting all of the intra-cluster edges in
G. In this way, it measures the strength of the intercluster connections. For a strongly clustered graph, it
therefore makes sense that ρ1 is small relative to λi,2
for all i. Further, if G consists of several such clusters
C1 , . . . , Ck , then, we can mathematically guarantee that
SP-MGM will correctly identify the clusters. Again, the
reader is referred to the appendix for details.

i∈U \S

max(ni − 1)
i6=m

where here if S or U \S is empty, the associated
intersection is interpreted to be {0, 1, . . . , n − 1}. It is
clear from the construction that {CS : S ⊆ U } forms a
partition of {0, 1, . . . , n − 1}.
The number of clusters in this partition is at least 2,
but can be up to 2n . However, in practice it is observed
that many of the clusters CS are repeated, and the
algorithm tends to find far fewer than 2n clusters. In
fact, the number of clusters created in this way tends
to be close to the actual right number of clusters. This
will be discussed further in the benchmark section
below. Note that we experimented with both using the
normalized and non-normalized Laplacian. We found
that the non-normalized Laplacian performs better for
the benchmark graphs that we considered. We note
that a common and effective strategy in applying spectral methods to large data sets consists of adopting a
hierarchical approach wherein spectral techniques are
used on subsets of the data that are then recombined
[e.g., see 15, 28, 24, 21].

3
3.1

2

3
4
5
6
7
8
9
10
11

Input : weight matrix A
number of modes k
Output: m clusters (2 ≤ m ≤ 2k )
Calculate non-normalized Laplacian L of A
Calculate k smallest positive eigenvalues
λ1 < · · · < λk of L
Calculate an eigenvector ~vj of λi
Set ~
uj = sort(~vj )
Choose `j so that uj(`j +1) − uj`j is largest
if vji < uj`j then
put node i in cluster Pj
else
put node i in cluster Qj
end
For each S ⊆ {1, . . . , k} define

! 
\
\
CS =
Pi ∩ 
Qi  ,
i∈S

LFR Benchmark

To test the quality of our clustering/communitydetection algorithm, we choose to use the LFR benchmark. Benchmark graphs consist of randomly generated graphs with built-in network ties. The properties of the output graph are controlled by a series of parameters. The parameters we consider and
vary are the total number of nodes n, the average
degree davg , the maximum degree dmax , and the
mixing parameter µ. Each node is given a degree
based on a power law distribution. The minimum
and maximum degrees in the network are chosen
so that the mean from the power law distribution
is davg . Nodes are assigned to clusters such that the
(weighted) fraction of intercluster edges is 0 ≤ µ ≤ 1.
Source code for generating LFR Benchmark graphs
is available at https://github.com/eXascaleInfolab/
LFR-Benchmark_UndirWeightOvp.
Note also that each of the algorithms described
above requires an important piece of user input:
namely the number of clusters to look for. In practice,
there are several ways of deciding how many clusters
in a graph to look for, such as the elbow technique, the
jump method, and silhouette analysis [12]. However,
these seem outside the scope of the current paper.
Instead, we supply each algorithm with the correct
number of clusters to look for. The exception to this is
SP-MGM, which does not take in a number of clusters

Algorithm 3: SP-MGM clustering

1

B ENCHMARK E XPERIMENTS

i∈U \S

2.3.1 A mathematical justification of the new algorithm
We provide a mathematical justification of the SPMGM algorithm by applying the theory of eigenvector
4

to look for. Instead, it takes in a number of eigenmodes
m. However, in practice, the number of clusters it finds
is m + 1, so we choose m accordingly.
3.2

The Quality of a Clustering Partition

After a clustering algorithm identifies one or more
clusters in a graph, it is natural to want to identify the
quality of the clusters. Are the vertices in a cluster naturally correlated in some fashion, or is the cluster only
an artifact of the specified clustering algorithm? While
there is no single criteria for reliably determining the
quality of a cluster, there exist several metrics which
may be used in combination to judge the reasonability
of ones choice of clusters. For a comparison of several
such metrics, see [1].
The LFR Benchmark produces a weighted graph,
along with a predetermined true value for the neighborhoods. To determine the performance of each clustering algorithm we use two scoring metrics: the normalized mutual info score (NMI) and the adjusted
rand score (ARS). In both metrics, the scores can range
between 0 and 1, with 1 representing perfect cluster
prediction and 0 representing largely independent
cluster prediction from the true clusters.
3.3

Fig. 2: The average of the normalized mutual info
score (NMI) and adjusted rand score (ARS) for various
values of the mixing parameter µ. The benchmark
graphs used here were all 50 nodes each with an
average degree of 5 and a maximum degree of 20.
The scores for each value of µ are averaged over 50
randomly generated benchmark graphs.

Results

We will now compare the performance of our algorithm to the performance of a traditional spectral clustering SP-KMEANS, as well as traditional KMEANS
and an agglomerative clustering algorithm. This collection of algorithms was chosen as a basis of comparison since they are notably diverse in their methodologies. At the same time, they are also popular enough
to have been implemented in the widely-used python
machine learning package scikit-learn [20]. While SPKMEANS has already been described, the reader is
referred to [22] for descriptions of KMEANS and the
agglomerative methods.
3.3.1 Variations in Mixing Parameter
The scores of clustering algorithms for the LFR benchmark on graphs with 50 and 100 nodes with varying
mixing parameter are shown in Figures 2 and 3. From
these figures, it is evident, and notably so, that SPMGM consistently outperforms the other clustering
algorithms considered for small networks over the
entire range of the mixing parameter considered.
Next, we examined the performance of SP-MGM
relative to the other algorithms considered while varying the average vertex degree, and the network size.
Figure 4 shows the results when average vertex degree
was varied. As the average network degree increased,
all algorithms considered showed an increase in score.
However, SP-MGM consistently out-performed the
other algorithms considered. Furthermore, SP-MGM

Fig. 3: The average of the normalized mutual info
score (NMI) and adjusted rand score (ARS) for various
values of the mixing parameter µ. The benchmark
graphs used here were all 100 nodes each with an
average degree of 5 and a maximum degree of 20.
The scores for each value of µ are averaged over 50
randomly generated benchmark graphs.

5

was the only algorithm to perform reasonably when
the average vertex degree was small, and the mixing
parameter was simultaneously low.
The scores of clustering algorithms for the LFR
benchmark on graphs with 50 and 100 nodes with
varying mixing number are shown in Figures 2 and
3. From these figures, it is evident, and notably so,
that SP-MGM consistently outperforms the other clustering algorithms considered for small networks over
the entire range of the mixing parameter considered.
For larger networks (not shown), SP-MGM performs comparably or better than the other algorithms
considered when the algorithms are scoring high
enough. In other words, when any of the algorithms
perform well, SP-MGM performs well as well. We
note here that when edges are distributed uniformly
across the nodes, clusters are not well-defined, and
as such, the computed clustering is likely not robust
(and likely arbitrary) and therefore not meaningful
or useful. Notably, for small values of the mixing
parameter on small graphs, the prediction-skill of SPMGM is far superior.

Fig. 4: The average of the normalized mutual info
score (NMI) and adjusted rand score (ARS) vs. average
vertex degree k . The benchmark graphs used here
all had mixing parameter µ = 0.1 and network size
n = 100. Each score is averaged over 50 randomly
generated benchmark graphs.

3.3.2 Variations in Average Degree and Network Size
Next, we examined the performance of SP-MGM relative to the other algorithms considered while varying
the average vertex degree, and the network size. Figure 4 shows the results when average vertex degree
was varied. As the average network degree increased,
all algorithms considered showed an increase in score.
However, SP-MGM consistently out-performed the
other algorithms considered. Furthermore, SP-MGM
was the only algorithm to perform reasonably when
the average vertex degree was small, and the mixing
parameter was simultaneously low.
Figure 5 shows the results when the network size
was varied. SP-MGM is seen to perform well over
a range of network sizes, with only the agglomerative clustering algorithm out-performing it for larger
graph sizes. We note, however, that the new algorithm tended to outperform its spectral counterpart
SP-KMEAN over the full range of network sizes considered.

Fig. 5: The average of the normalized mutual info
score (NMI) and adjusted rand score (ARS) for various
network sizes n. The benchmark graphs used here all
had mixing parameter µ = 0.1 and average connectivity k = 10. Each score is averaged over 100 randomly
generated benchmark graphs. The data shows that the
spectral gap method tends to have comparable or better performance than other methods for a wide range
of network sizes. However, performance worsens for
larger network sizes.

3.3.3 The Number of Clusters SP-MGM Detects
From the description of the SP-MGM algorithm, it
is seen that we need to specify only the number m
of eigenmodes to consider. For this reason, a priori
the algorithm should be able to predict anywhere
from 2 and min{2m , n} clusters, where here n is the
total number of vertices in the graph. However, in
practice the algorithm almost always produces m + 1
clusters, at least for all the test graphs examined in this
paper. This result is quite surprising, since from an
intuitive perspective one would suspect the number
of clusters produced to be far more volatile. It seems
6

4.3

like this is additional computational evidence that the
“maximum gaps" we are computing must have a lot to
do with the structure of the graphs themselves, closely
tied to the clusters.

Clustering Results

While ONI is a conventional and simple measure
that is used to classify the state of ENSO, the phenomenon itself is more complicated. That is, although
4 C LUSTERING OF SST DATA IN E L N IÑO 3.4 the ENSO phenomenon can be characterized as a lowdimensional dynamical phenomenon occuring in an
In this section we consider the application of the new
infinite-dimensional dynamical system, it is unlikely
clustering algorithm to a real world data set in the
to be fully or adequately characterized by just ONI. As
context of the climate system. Specifically we consider
such, we conduct a clustering analysis of the spatiallythe El Niño Southern Oscillation (ENSO) since it is one
extended sea surface temperature anomaly (SSTA)
of the most important modes of interannual variability
dataset T(x, y, t) where t runs from years 1955 through
of the climate system. It can be characterised in terms
2016 and x and y correspond to the Niño 3.4 region
of the sea surface temperature anomaly occuring in a
described above. The data vector for each year itself
region of the equatorial Pacific Ocean. The ENSO temis a 3-month average starting from November of the
perature anomaly roughly fluctuates between three
previous year.
phases: a warm El Niño phase, a cold La Niña phase,
In order to conduct cluster analysis, we next transand a neutral phase. The phases of ENSO are correform
the SSTA dataset into a graph whose vertices
lated with changes in precipitation in various regions
are
the
yearly SSTA datapoints embedded in a higharound the globe, with strong correlations in coastal
dimensional
space, the number of dimensions corPacific regions. As such, predictions of ENSO have
responding
to
the number of distinct spatial locahigh socio-economic value, Luckily, ENSO also happens to be one of the few modes of interannual vari- tions in the Niño 3.4 region. Next we consider an
ability that has a useful level of verified predictability. edge structure that reflects the mutual similarity of
pairs of vertices. The resultant similarity matrix has
the form Aij = exp(−dij /σ) for i 6= j , where
4.1 Sea Surface Temperature Data
T
−1 ~
~
We use version 4 of the Extended Reconstructed dij = T (x, y, i) Σ(x, y) T (x, y, j), and where Σ the
matrix for the period
Sea Surface Temperature (ERSST) dataset available spatially-extended covariance −1
considered.
Note
that
Σ(x,
y)
is also sometimes
at https://www.esrl.noaa.gov/psd/data/gridded/
called
the
precision
matrix.
For
simplicity,
we consider
data.noaa.ersst.v4.html. The data has a spatial
a
diagonal
form
for
Σ(x,
y)
that
comprises
the variresolution of 2 degree, and while the monthlyaveraged data that we use are available from 1854- ances of the SSTA at each of the spatial locations in
present, we focus attention on data over the years the Niño 3.4 region. The parameter σ toggles controls
−2
1955-2016, given lower uncertainties over this more localization of the similarity matrix; we use σ = 10 .
In Figure 8, the nodes are colored-coded based on
recent (instrumented) period. For further description
of the data, the reader is referred to https:// cluster-membership and the nodes are placed at the
www.ncdc.noaa.gov/data-access/marineocean-data/ two-dimensional point corresponding to the year (xaxis) and ONI of that year (y-axis). Finally, the thickextended-reconstructed-sea-surface-temperature-ersst-v4.
ness of the edges reflect similarity between the pair
4.2 El Niño Events
of vertices they connect. In this figure, it is seen that
The average sea surface temperature is calculated there is good correlation between the clusters and the
using a 30-year average of the NOAA ERSST data conventional ONI (y-axis). Thus, the new clustering
starting Jan 1, 1950. SST anomaly is calulated relative algorithm SP-MGM is seen to successfully distinguish
between El Niño , La Niña , and neutral events.
to this mean [27].
Next, it is seen that the SP-MGM clustering alWe consider the data from ERSSTv4 in the Niño 3.4
region (the region between 170W and 120W longitude gorithm identifies 7 distinct clusters and that cluster
and between -5N and 5N latitude), as shown in Figure membership is not solely determined by ONI. Indeed,
6. El Niño events are classified by a persistent large that there are different flavors of ENSO besides El
average temperature anomaly in this region. Quanti- Niño and La Niña is now widely recognized and how
tatively, El Niño events are determined by the Oceanic these different flavors influence climate in different
Niño Index (ONI), which is a 3-month running mean regions is an active area of research. The cluster patof the SST anomaly spatially averaged over the Niño terns we find are consistent, e.g., with the patterns
3.4 region. A Niño event is characterized by ONI found in [10] who use a neural-network based analysis
values greater than 0.5 degrees C for at least three in conjunction with other statistical distinguishability
months in a row. The years of El Niño events between tests. We note, however, that there are significant dif1950 and 2000 may be seen from the graph in Figure ferences in the data used by [10] in that they consider
(7) below.
a six month average starting in September and a much
7

Fig. 6: Monthly averaged sea surface temperature on Jan 1, 2017. The region surrounded by the dashed line is
the Niño 3.4 region, which lies in 170W-120W, 5S-5N.
used to automatically detect an appropriate number of
clusters for a specified graph. Finally, we applied this
algorithm to analyze and identify a variety of flavors
of the El Niño Southern Oscillation using spatially
extended sea surface temperature data.

more extensive region of the Pacific.

5

C ONCLUSIONS

In this paper, we presented a new spectral clustering
algorithm which we call the Spectral-Maximum Gap
Method SP-MGM. This algorithm is based on identifying gaps in the structure of eigenvectors. We then
went on to provide a mathematical justification for the
algorithm. While it would be interesting to see where
this algorithm fits in a systematic comparison of spectral algorithms, such as in [30], such a comparison is
outside the scope of the current paper. Further, as with
other spectral schemes, this method can be applied
to large datasets in conjunction with a hierarchical
approach as in [15, 28, 24, 21, and others].
Next, we examined the performance of this algorithm in comparison to a few other popular algorithms
using the LFR benchmark graphs. Results showed that
the SP-MGM algorithm performs either comparbly or
better than its counterparts in a variety of parameter
regimes. We also demonstrated how SP-MGM may be

ACKNOWLEDGEMENT
This work was performed during one authors visit to
CNLS at Los Alamos National Lab.

A PPENDIX
Eigenvector Perturbation
The eigenvalue perturbation problem is the problem
of estimating the eigenvectors and eigenvalues of a
perturbed matrix T = T0 + T1 based on the eigenvalues of the matrix T0 , where here  is a small positive
number which estimates the size of the perturbation.
In such a case one should rightfully assume that the
eigenvalues and eigenvectors of T and T0 are nearly
the same. In this section, we briefly recount the idea
8

Fig. 7: Values of the Oceanic Niño Index (ONI) between 1950 and 2000. Years where an El Niño event occurred
are marked. Since impacts of El Niño are largest in the (northern-hemisphere), winter some of the events span
the calendar year boundary.

~vi of T satisfying

of eigenvector perturbation analysis, deriving in particular the linear perturbation equations for the eigendata of a perturbed matrix. This equation specifically
relates the eigenvectors of T0 and T = T0 + T1 for
 sufficiently small. Our exposition is loosely based
on that found in the standard source [5]. However,
standard perturbation analysis formulas seems to emphasize the reqirement that the eigenvalues of T0 are
distinct, whereas our brief derivation shows that this
is not the case, as long as the formulas are appropriately adjusted. This is important for our application to
graph theory, where we wish T0 to be the Laplacian of
a disconnected graph.

~vi = ~ui +

h~uj , T1 ~ui i
~ui + O(2 ).
λ
−
λ
j
i
6=λi

X
j:λj

(1)

Furthermore the eigenvalue λ0i of ~vi is given by

λ0i = λi + h~uj , T1 ~ui i + O(2 ).

(2)

where here the braces h·, ·i denote the usual inner
product on Rn .
To prove the above two equations, we can use a
standard but important linear algebra tool: the resolvent of a matrix. Given any matrix L, the resolvent
R(L; z) is a matrix-valued rational function on the
complex plane, defined by R(L; z) = (Iz − L)−1 . One
useful property of the resolvent is that acts like a projection operator onto the eigenspaces of L. Specifically,
if λ is an eigenvalue of R(L; z) and C is a closed

Assume that T0 , T1 are symmetric n × n matrices
and let ~
u1 , . . . , ~un be an orthonormal eigenbasis for
T1 , with T1 ~ui = λi ~ui . Then the linear perturbation
equation says that for all i there exists an eigenvector
9

Fig. 8: Results of clustering NDJ SST states in the Niño 3.4 region. Nodes are colored based on the cluster to
which they belong. Years with similar ONI tend to be clustered together.
will be an eigenvector with eigenvalue λ0i . However,
using our expression for R(T, z) we see that

contour in the complex plane, then for any constant
vector ~v ∈ Rn , the integral
I
1
R(L; z)~v dz
2πi C

R(T, z)~ui =

∞
X

m (R(T0 ; z)T1 )m

m=0

1
~ui
z − λi



n
X
1 
h~uj , T ~ui i 
=
~ui +
~uj + O(2 ).
z − λi
(z
−
λ
)
j
j=i

is an eigenvector of L with eigenvalue λ, as long as it
is nonzero.
Now since T = T0 + T1 , a simple calculation
shows

Therefore by Cauchy’s residue theorem
X h~
uj , T1 ~ui i
~vi = ~ui +
~ui + O(2 ).
λ
−
λ
j
i
j:λ 6=λi

R(T ; z) = (I − R(T0 ; z)T1 )−1 R(T0 ; z)
∞
X
=
m (R(T0 ; z)T1 )m R(T0 ; z).

j

This proves Equation 1. To get the associated eigenvalue λ0i , we can use the fact that

m=0

For each i, let λ0i be an eigenvalue of T closest to λi .
Take Ci to be a closed contour in the complex plane
containing λi and λ0i , but no other (distinct) eigenvalue
of T or T0 . Then we know that
I
1
~vi :=
R(T, z)~ui dz
2πi Ci

λ0i h~ui , ~vi i = h~ui , T ~vi i.
From our expression for ~vi , we have h~
ui , ~vi i = 1 +
O(2 ) and h~ui , T ~vi i = λi + h~ui , T1 ~ui i + O(2 ). Thus

λ0i = λi + h~uj , T1 ~ui i + O(2 ).
This proves Equation 2.
10

to order O(2 ), where  = ρ1 /ρ0 . Since fij is supported on Ci and L1 sends fm to a function supported
on the complement of Cm , we find that to order (2 )

Mathematical Justification for SP-MGM
Let G be a graph with n vertices and let C1 , . . . , Ck be
the natural subclusters of G. To show mathematically
how SP-MGM detects the clusters C1 , . . . , Ck from
the graph G, we will consider two subgraphs of G,
which we will call G0 and G1 , and which will be
formed by deleting the intercluster bonds and intracluster bonds of G, respectively. Also to simplify the
notation, throughout this section we will think about
the Laplacian of G as acting not on Rn , but rather
on the space of real-valued functions on G. For this
reason, we will switch from talking about eigenvectors
to talking about eigenfunctions. We will denote the
Laplacians of G, G0 and G1 as L, L0 and L1 . Clearly
L = L0 + L1 .
To better understand why SP-MGM should work
well, we need to consider in more detail the eigenstructure of the Laplacian L0 of G0 . The eigenfunctions of G0 clearly restrict to each cluster Ci to an
eigenfunction on Ci . Each cluster Ci is connected,
and therefore has a kernel spanned by its indicator
function. Thus an orthonormal basis for the kernel of
L0 is given by
( −1/2
ni
, x ∈ Ci
fi (x) =
0
x∈
/ Ci

0
kfm
(x) − fm (x)k ≤ max(ni − 1)
i6=m

Thus we arrive at the conclusion that if  = ρ1 /ρ0 is
small and for all m we have

max(ni − 1)
i6=m

ρ1
1
< 1/2 .
λi,2
4nm

Then SP-MGM will correctly identify the desired clusters.

B IBLIOGRAPHY
[1]

[2]

[3]

for i = 1, . . . , k . More generally, let λi,1 ≤ · · · ≤ λi,ni
be the eigenvalues of Ci . Then the λij are also eigenvalues of L0 and there exists an orthonormal eigenbasis for L0 such that fi,j (x) has eigenvalue λi,j and is
supported on Ci . By definition, fi,1 (x) = fi (x).
Assuming that the clusters are really, honestly clusters and not just part of some artificial partition of G,
we should expect the spectral radius ρ1 of L1 to be
small compared to the spectral radius ρ0 of L0 . Then
the eigendata of L and L0 will be related via the eigenvector perturbation theory discussed above. Each of
the fi (x) will correspond to an eigenfunction fi0 (x)
of the matrix L via this perturbation theory, and the
associated eigenvalues should be expected to be the
smallest k eigenvalues of the matrix L. Moreover, this
eigenvalue perturbation theory gives us an estimate of
fi (x)0 relative to fi (x). If we can show that
1
kfi (x) − fi0 (x)k∞ < 1/2 ,
4ni

[4]

[5]
[6]

[7]

[8]

then the location of the maximum gap in fi (x) and
the location of the maximum gap in fi (x) must be the
same! In this case, SP-MGM is guaranteed to identify a
partition which correctly separates one of the clusters
from the others.
From our eigenvalue perturbation theory, we know
that for m = 1, . . . , k
0
fm
(x) = fm (x) +

ρ1
.
λi,2

[9]
[10]

[11]

ni
k X
X
hfi,j (x), L1 fm1 (x)i
fi,j (x)
λi,j
i=1 j=2

11

Hélio Almeida, Dorgival Guedes, Wagner Meira,
and Mohammed J Zaki. Is there a best quality
metric for graph clusters? In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pages 44–59. Springer, 2011.
Thomas Bühler and Matthias Hein. Spectral clustering based on the graph p-laplacian. In Proceedings of the 26th Annual International Conference on
Machine Learning, pages 81–88. ACM, 2009.
João Corte-Real, Budong Qian, and Hong Xu.
Regional climate change in portugal: precipitation variability associated with large-scale atmospheric circulation. International Journal of Climatology, 18(6):619–635, 1998.
Sergio Decherchi, Simone Tacconi, Judith Redi,
Alessio Leoncini, Fabio Sangiacomo, and Rodolfo
Zunino. Text clustering for digital forensics analysis. Computational Intelligence in Security for Information Systems, pages 29–36, 2009.
Demmel. Applied numerical linear algebra. 1997.
Miroslav Fiedler. A property of eigenvectors of
nonnegative symmetric matrices and its application to graph theory. Czechoslovak Mathematical
Journal, 25(4):619–633, 1975.
Sally Galbraith, James A Daniel, and Bryce Vissel.
A study of clustered data and approaches to its
analysis. Journal of Neuroscience, 30(32):10601–
10608, 2010.
Robert M Haralick and Linda G Shapiro. Image
segmentation techniques. Computer vision, graphics, and image processing, 29(1):100–132, 1985.
Anil K Jain and Richard C Dubes. Algorithms for
clustering data. Prentice-Hall, Inc., 1988.
Nathaniel C Johnson. How many enso flavors can
we distinguish? Journal of Climate, 26(13):4816–
4827, 2013.
Alexandros Karatzoglou, Alex Smola, Kurt
Hornik, and Achim Zeileis. kernlab – an S4 package for kernel methods in R. Journal of Statistical
Software, 11(9):1–20, 2004.

905, 2000.
[26] Karsten Steinhaeuser, Nitesh Chawla, and Auroop Ganguly. Comparing predictive power in
climate data: Clustering matters. Advances in
Spatial and Temporal Databases, pages 39–55, 2011.
[27] Kevin E Trenberth. The definition of el nino.
Bulletin of the American Meteorological Society,
78(12):2771–2777, 1997.
[28] Frederick Tung, Alexander Wong, and David A
Clausi. Enabling scalable spectral clustering
for image segmentation. Pattern Recognition,
43(12):4069–4076, 2010.
[29] Yurdanur Unal, Tayfun Kindap, and Mehmet
Karaca. Redefining the climate zones of turkey
using cluster analysis. International journal of
climatology, 23(9):1045–1055, 2003.
[30] Deepak Verma and Marina Meila. A comparison
of spectral clustering algorithms. University of
Washington Tech Rep UWCSE030501, 1:1–18, 2003.
[31] Ulrike Von Luxburg. A tutorial on spectral clustering. Statistics and computing, 17(4):395–416,
2007.

[12] Trupti M Kodinariya and Prashant R Makwana.
Review on determining number of cluster in kmeans clustering. International Journal, 1(6):90–95,
2013.
[13] Andrea Lancichinetti, Santo Fortunato, and Filippo Radicchi. Benchmark graphs for testing
community detection algorithms. Physical review
E, 78(4):046110, 2008.
[14] Robert Lund and Bo Li. Revisiting climate region definitions via clustering. Journal of Climate,
22(7):1787–1800, 2009.
[15] Ulrike V Luxburg, Olivier Bousquet, and Mikhail
Belkin. Limits of spectral clustering. In Advances
in neural information processing systems, pages 857–
864, 2005.
[16] Marina Meila and Jianbo Shi. Learning segmentation by random walks. In Advances in neural information processing systems, pages 873–879, 2001.
[17] Nina Mishra, Robert Schreiber, Isabelle Stanton,
and Robert E Tarjan. Clustering social networks.
In WAW, volume 4863, pages 56–67. Springer,
2007.
[18] Maria CV Nascimento and Andre CPLF De Carvalho. Spectral methods for graph clustering–a
survey. European Journal of Operational Research,
211(2):221–231, 2011.
[19] Andrew Y Ng, Michael I Jordan, and Yair Weiss.
On spectral clustering: Analysis and an algorithm. In Advances in neural information processing
systems, pages 849–856, 2002.
[20] F. Pedregosa, G. Varoquaux, A. Gramfort,
V. Michel, B. Thirion, O. Grisel, M. Blondel,
P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher,
M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine
Learning Research, 12:2825–2830, 2011.
[21] Diego H Peluffo-Ordónez, John A Lee, and
Michel Verleysen. Short review of dimensionality
reduction methods based on stochastic neighbour
embedding. In Advances in Self-Organizing Maps
and Learning Vector Quantization, pages 65–74.
Springer, 2014.
[22] Lior Rokach and Oded Maimon. The data mining
and knowledge discovery handbook: A complete
guide for researchers and practitioners, 2005.
[23] Satu Elisa Schaeffer. Graph clustering. Computer
science review, 1(1):27–64, 2007.
[24] Theodoros Semertzidis, Dimitrios Rafailidis,
Michael G Strintzis, and Petros Daras. Large-scale
spectral clustering based on pairwise constraints.
Information Processing & Management, 51(5):616–
624, 2015.
[25] Jianbo Shi and Jitendra Malik. Normalized cuts
and image segmentation. IEEE Transactions on
pattern analysis and machine intelligence, 22(8):888–

W.R. Casper received BS degrees in
mathematics and physics, as well as an
MS in mathematics from North Dakota
State University in 2010. He later received
a PhD in mathematics from the University of Washington, Seattle in 2017. He is
currently a postdoctoral researcher in the
Department of Mathematics at Louisiana
State University in Baton Rouge. His research interests are diverse, and include
integrable systems, algebraic geometry,
noncommutative algebra, geophysical fluid dynamics, computational physics and machine learning.

Balu Nadiga is a scientist at the Los
Alamos National Lab. With a specialization in fluid dynamics, he has worked
on problems ranging from climate and
ocean circulation to turbulence modeling and multiphase flows. His other interests include dynamical systems, uncertainty quantification and Bayesian analysis, and using machine learning for turbulence modeling.

12

