In general it does not make much sense to cluster features. In an ideal world for your features to be the best they can be they should actually be independent, thus there should be no relationship between them. Typically when we talk about clustering it is clustering the instances. To attribute some associative labels to a subset of the instances based on the similarity of their feature values.
Many clustering algorithms exist, I would say that the most popular is K-means however spectral clustering and Gaussian mixtures are also frequently used. As always, each algorithm is best suited for a specific type of dataset, it is up to you to choose which is best suited, or you can just try all of them and see which is best.
Here you can find a list of clustering algorithms with their respective usecases. Always use the libraries when you want to implement standard algorithms, they are highly optimized. But for education sake it is good to look at what is happening.
I will describe a homebrew version of the K-means algorithm such that you can understand what is happening under the hood and perhaps you will see why we cluster instances and not features.
The K-means algorithm
import numpy as np
import matplotlib.pyplot as plt
First we will make some artificial data. These will consist of $n$ 2D Gaussian clusters with a given mean and variance. Here $n=5$ and we will have 300 instances per Gaussian distribution.
params = [[[ 0,1], [ 0,1]],
[[ 5,1], [ 5,1]],
[[-2,5], [ 2,5]],
[[ 2,1], [ 2,1]],
[[-5,1], [-5,1]]]
n = 300
dims = len(params[0])
data = []
y = []
for ix, i in enumerate(params):
inst = np.random.randn(n, dims)
for dim in range(dims):
inst[:,dim] = params[ix][dim][0]+params[ix][dim][1]*inst[:,dim]
label = ix + np.zeros(n)
if len(data) == 0: data = inst
else: data = np.append( data, inst, axis= 0)
if len(y) == 0: y = label
else: y = np.append(y, label)
num_clusters = len(params)
print(y.shape)
print(data.shape)
(1500,)
(1500, 2)
So now we have 1500 instances (records) in 2D space. This can be extended to any number dimensions. 2 is easiest to plot.
plt.scatter(data[:,0], data[:,1])
plt.show()
I will paste the entire algorithm at the bottom of the answer for a quick copy and paste but here I will go through its different parts so you can see how it works. The algorithm goes as follows: first we initialize some centroids within the range of our data. These are the red dots in the image below
def train(self, data, verbose=1):
shape = data.shape
ranges = np.zeros((shape[1], 2))
centroids = np.zeros((shape[1], 2))
for dim in range(shape[1]):
ranges[dim, 0] = np.min(data[:,dim])
ranges[dim, 1] = np.max(data[:,dim])
if verbose == 1:
print('Ranges: ')
print(ranges)
centroids = np.zeros((self.k, shape[1]))
for i in range(self.k):
for dim in range(shape[1]):
centroids[i, dim] = np.random.uniform(ranges[dim, 0], ranges[dim, 1], 1)
if verbose == 1:
print('Centroids: ')
print(centroids)
plt.scatter(data[:,0], data[:,1])
plt.scatter(centroids[:,0], centroids[:,1], c = 'r')
plt.show()
Then we will calculate the distance from each instance (record) to each centroid.
distances = np.zeros((shape[0],self.k))
for ix, i in enumerate(data):
for ic, c in enumerate(centroids):
distances[ix, ic] = np.sqrt(np.sum((i-c)**2))
Then we will attribute each instance as belonging to the centroid it is closest to.
labels = np.argmin(distances, axis = 1)
Now we will update the position of the new centroids by finding the mean position of all instances which are closest to the given centroid in each dimension.
new_centroids = np.zeros((self.k, shape[1]))
for centroid in range(self.k):
temp = data[labels == centroid]
if len(temp) == 0:
return 0
for dim in range(shape[1]):
new_centroids[centroid, dim] = np.mean(temp[:,dim])
We then repeat this process until the centroids no longer move significantly. Usually if the difference in position for all the centroids is less than machine epsilon we consider the algorithm to have converged.
if np.linalg.norm(new_centroids - centroids) < np.finfo(float).eps:
print("DONE!")
break
For this dataset it took 16 iterations for convergence, this is the final result.
What we see here is that each instance is grouped into a cluster with similar attributes as itself. For example if the dimensions of our data represented height (x-axis) and weight (y-axis), then we can group people into 5 different BMI indices.
- Purple: short and skinny
- Green: short and fat
- Blue: mid height and mid weight
- Yellow: Mid height and fat
- Turquoise: Tall and mid weight.
Each member in our population belongs to a specific BMI index which is clustered based on two measurable attributes he possesses (features), his height and his weight.
The full K-means algorithm
This is the algorithm
class Kmeans(object):
def __init__(self, k=1):
self.k = k
def train(self, data, verbose=1):
shape = data.shape
ranges = np.zeros((shape[1], 2))
centroids = np.zeros((shape[1], 2))
for dim in range(shape[1]):
ranges[dim, 0] = np.min(data[:,dim])
ranges[dim, 1] = np.max(data[:,dim])
if verbose == 1:
print('Ranges: ')
print(ranges)
centroids = np.zeros((self.k, shape[1]))
for i in range(self.k):
for dim in range(shape[1]):
centroids[i, dim] = np.random.uniform(ranges[dim, 0], ranges[dim, 1], 1)
if verbose == 1:
print('Centroids: ')
print(centroids)
plt.scatter(data[:,0], data[:,1])
plt.scatter(centroids[:,0], centroids[:,1], c = 'r')
plt.show()
count = 0
while count < 100:
count += 1
if verbose == 1:
print('-----------------------------------------------')
print('Iteration: ', count)
distances = np.zeros((shape[0],self.k))
for ix, i in enumerate(data):
for ic, c in enumerate(centroids):
distances[ix, ic] = np.sqrt(np.sum((i-c)**2))
labels = np.argmin(distances, axis = 1)
new_centroids = np.zeros((self.k, shape[1]))
for centroid in range(self.k):
temp = data[labels == centroid]
if len(temp) == 0:
return 0
for dim in range(shape[1]):
new_centroids[centroid, dim] = np.mean(temp[:,dim])
if verbose == 1:
plt.scatter(data[:,0], data[:,1], c = labels)
plt.scatter(new_centroids[:,0], new_centroids[:,1], c = 'r')
plt.show()
if np.linalg.norm(new_centroids - centroids) < np.finfo(float).eps:
print("DONE!")
break
centroids = new_centroids
self.centroids = centroids
self.labels = labels
if verbose == 1:
print(labels)
print(centroids)
return 1
def getAverageDistance(self, data):
dists = np.zeros((len(self.centroids),))
for ix, centroid in enumerate(self.centroids):
temp = data[self.labels == ix]
dist = 0
for i in temp:
dist += np.linalg.norm(i - centroid)
dists[ix] = dist/len(temp)
return dists
def getLabels(self):
return self.labels
To use the algorithm use the following, where data is artificial data we made above but can also be any numpy matrix where the records are the rows and the features are the columns.
kmeans = Kmeans(5)
kmeans.train(data)
