Sprint 1: K-means spectral clustering

For the past few weeks I’ve been performing the usual tasks of getting my development environment set up, interfacing with the Apache Mahout repository, learning more about Maven’s idiosyncrasies (gotta love missing environment variables), and building Mahout from scratch. I’ve also been working on the sprint 1 deliverable: a working implementation of k-means spectral clustering, a vanilla version of the final EigenCuts algorithm I’m poised to deliver by the end of the summer.

K-means spectral clustering is much simpler than the full EigenCuts algorithm in that its computations effectively cease after performing eigen-decompositions on the normalized Laplacian matrix. Again, I’ll delve into the specifics of EigenCuts in a later post, but what is outlined in the previous post is, essentially, all this algorithm requires. I’ll run through the steps here, along with some choice Python code that I’ve been working on as a prototype before implementing it in the Mahout framework.

But first, I need to define some terms that I’ll be using.

• R: Raw data matrix. This is the K × n-dimensional data the user is interested in clustering.
• S: Similarity matrix, a K × K transformation of R by showing how “related” each point is, pairwise. This “relation” function can be anything, from pixel intensity to radial Euclidean distance.
• A: Adjacency/Affinity matrix, also K × K, this time a transformation of S by applying a k-nearest neighbor filter to build a representation of the graph (or, for a fully-connected graph, A = S). This matrix is critical to our calculations.
• D: Diagonal degree matrix, also K × K, formed by summing the degree of each vertex and placing it on the diagonal.
• L: Normalized and symmetric Laplacian matrix, formed in an operation with A and D. This is the matrix on which we will perform eigen-decompositions.
• U: Matrix of eigenvectors of L. We’ll perform k-means clustering on their components.

Step 1: Process raw data

This is a somewhat tricky process by itself. The data can be comprised of some arbitrary K points, all in some arbitrary n-dimensional space. Each data point in this abstract matrix R needs to be processed into a pairwise similarity score via a “relation” function; in most cases, this function takes the form

exp{ -(xixj)2 / c }

where each x is the raw data point, and c is some scaling factor that can be constant or relative to the data itself (such as the median). These processes construct the similarity matrix S.

Here is a Python implementation of the construction of S, where the scaling factor c = 2σ2, where σ = pg. p is a constant of 1.5 that is used in the EigenCuts paper, and g is the raw median of the Euclidean distances of all neighboring data points.

```for i in range(0, len(points)):
row = []
for j in range(0, len(points)):
# scaled pairwise comparison
if i == j:
row.append(0)
else:
d = 0
diffs = []
for k in range(0, len(points[j])):
diff = ((points[j][k] - points[i][k]) ** 2)
diffs.append(diff)
d += diff
row.append(math.exp((-d) / \
(((numpy.median(diffs) * 1.5) ** 2) * 2)))
S.append(row)
```

This step is a simple transformation of the similarity matrix we calculated in the previous step, and the specifics depend on the needs of the user. Specifically, if the user requires a fully-connected graph, then in this case the adjacency matrix A will be exactly the same as the similarity matrix S. Within the EigenCuts paper, the choice was made to use an 8-neighborhood, or take the 8 closest points out of all the neighboring points.

Given our current relation function, taking the 8 nearest points will equate to taking the 8 points with the largest values in the current row of S. Here is a Python implementation:

```if knn < len(points):
for i in range(0, len(S)):
A.append([])
row = S[i][:]
row.sort(reverse=True)
row = row[:knn]
for j in range(0, len(S[i])):
if S[i][j] in row:
A[i].append(S[i][j])
else:
A[i].append(0)
else:
A = S
```

Step 3: Construct diagonal degree matrix

This step is probably the most straightforward: simply sum the degrees of each vertex and place the value along the diagonal. Again, a Python implementation:

```for j in range(0, len(S)):
D.append([])
dj = 0
index = 0
for i in range(0, len(S[j])):
D[j].append(0)
if i == j:
index = i
dj += A[j][i]
D[j][index] = dj
```

Step 4: Construct normalized Laplacian

Here’s where it all starts coming together. We now want to build the normalized Laplacian matrix out of the matrices D and A that we previously built. This is accomplished via the formula from the previous post:

This is still pretty simple matrix algebra; the Python package numpy can handle this very well.

Step 5: Perform eigen-decomposition

Another straightforward computation: using the matrix L, we want to decompose it into its eigenvectors and eigenvalues. We’ll use the eigenvectors (or the M largest ones, again depending on the user’s need) to build a matrix U of column eigenvectors of L.

Step 6: Perform k-means clustering

Using the matrix U of eigenvectors, we perform a k-means clustering on the rows, essentially creating a new data set of M dimensions, with the nth component of each eigenvector acting as a single data point and, by proxy, representing the nth data point in R, S, or A. Essentially, we have created a “conglomerate” data set of the eigenvector components. We perform k-means clustering on this new K × M data set, and whichever cluster the nth “conglomerate” point is assigned, the corresponding nth data point in the original set also receives the same cluster assignment.

Step 7: ???

Take what has been written in Python and finish translating it on the Mahout framework!

Oh hai!

3 Responses to Sprint 1: K-means spectral clustering

1. sirus says:

For the order of eigenvectors should it be using the absolute value?

• magsol says:

Yes. The absolute values of the eigenvalues should be ordered in decreasing magnitude, such that: 1 ≥ |λ1| ≥ |λ2| ≥ … ≥ |λn|

2. Rose says:

Hi! I am new to Mahout and I have been trying to use the K-means Spectral Clustering, but I ran into the problem described in the comments from the jira MAHOUT-524: the Lancsoz solver tries to input the output of the VectorMatrixMultipliction as a “calculations/laplacian-166/tmp/data” file, instead of the “calculations/laplacian-166/part-m-00000″.
From what I read in the comments, there were some solutions but I could not really understand if currently there is a running version of the Spectral Clustering.
Thank you.