Implementing k-means clustering in Ruby

Inspired by the partitioning problem I set about to implement a well known algorithm, k-means clustering, from memory, for fun! ... and, for science... Interestingly, this is somehow the opposite of the partitioning problem. In the partitioning problem we tried to maximize variation of categorical variables within groups, whereas here we're trying to find groups of elements that are the most "similar" to each other in a n-dimensional continuous space.

The main idea of k-means is the following - if you have a bunch of elements and a given number of clusters:

1. Create the initial clusters randomly within the space spanned by the elements (typically (always?) you would pick randomly from your elements).
2. Lob all elements into the cluster with the nearest center (using some Euclidean distance metric typically).
3. Recenter each cluster on the average of its elements, 
4. If necessary move the elements to their now nearest clusters. 
5. Repeat the "re-centering" and moving of elements until we have a "stable" (enough) result.

Simple enough, but I don't think I have ever implemented it first hand. (At least not in Ruby.)

First I implemented some simple data structures, Points:

	class Point
	def initialize(coords)
	@coords = coords
	end
	def distance_to(point)
	Math.sqrt( @coords.keys.inject(0) { |sum, key| sum + (@coords[key] - point.send(key))**2 } )
	end
	def to_s
	"#{@coords.values.join(", ")}"
	end
	def dimensions
	@coords.keys
	end
	def method_missing(m, *args)
	if m.to_s.end_with? '='
	@coords[m.to_s.chop.to_sym]= args.first
	else
	@coords[m]
	end
	end
	end

view raw point.rb hosted with ❤ by GitHub

I use the snazzy method missing as a way to reference named dimensions of points easily, e.g. point.x and point.y (or point.z, point.t, point.w1, point.hausdorff, or point.scary).

Clusters are simply:

	class Cluster
	attr_reader :center
	def initialize(center)
	@center = center
	@points = []
	@moved = true
	end
	def add_point(point)
	@points << point
	end
	def update_center(delta = 0.001)
	@moved = false
	averages = {}
	@center.dimensions.each do |dimension|
	averages[dimension] =
	@points.inject(0.0) {|sum, point| sum + point.send(dimension)} /
	@points.length unless @points.length == 0
	end
	unless Point.new(averages).distance_to(@center) < delta
	@center = Point.new(averages)
	@moved = true
	end
	end
	def clear_points
	@points = []
	end
	def collect(dimension)
	@points.collect {|p| p.send(dimension) }
	end
	def distance_to(point)
	@center.distance_to point
	end
	def number_of_points
	@points.length
	end
	def to_s
	"#{@center.to_s}: #{number_of_points} points, cost: #{ @points.inject(0) { |sum, point| sum + point.distance_to(@center) }}"
	end
	def moved?
	@moved
	end
	end

view raw cluster.rb hosted with ❤ by GitHub

The definition of "moved" was the first design choice I could not remember entirely - is it when a certain number of points "jump" from a cluster to another or when the cluster center moves more than a certain distance? (Most of the time I would assume this boils down to the same.) I went for the latter. The clusters themselves keep track of if they have moved or not.

I also implemented a special case of points, 2D points:

	class Point2D < Point
	def initialize(x, y = nil)
	y.nil? ? super(x) : super({x: x, y: y})
	end
	end

view raw point2d.rb hosted with ❤ by GitHub

Last helper method in here is simply to plot these clusters (in 2D), and I use gnuplot like so:

	def cluster_plot_2D(clusters, seq = 1, title = "Random noise", get_x = 'x', get_y = 'y')
	if defined? Gnuplot
	Gnuplot.open do |gp|
	Gnuplot::Plot.new(gp) do |plot|
	plot.terminal "png"
	plot.output File.expand_path("../outfiles/clusters-#{seq}.png", __FILE__)
	plot.title title
	# Plot each cluster's points
	clusters.each do |cluster|
	# Collect all x and y coords for this cluster
	x = cluster.collect(get_x.to_sym)
	y = cluster.collect(get_y.to_sym)
	# Plot w/o a title (clutters things up)
	plot.data << Gnuplot::DataSet.new([x,y]) do |ds|
	ds.notitle
	end
	end
	# Plot each cluster's centers
	x = clusters.collect {|p| p.center.send(get_x.to_sym) }
	y = clusters.collect {|p| p.center.send(get_y.to_sym) }
	plot.data << Gnuplot::DataSet.new([x,y]) do |ds|
	ds.notitle
	end
	end
	end
	end
	end

view raw cluster_plot_2d.rb hosted with ❤ by GitHub

The real meat is in the k-means method itself of course:

	def k_means (points, k = 5, delta = 0.001, plot_on = false)
	k = points.length if points.length < k
	clusters = []
	k.times { clusters << Cluster.new(points.sample) }
	iterations = 0
	while (clusters.any?(&:moved?))
	clusters.each(&:clear_points)
	points.each do |point|
	shortest = Float::INFINITY # requires Ruby 1.9.2 or later
	cluster_found = nil
	clusters.each do |cluster|
	distance = cluster.distance_to(point)
	if distance < shortest
	cluster_found = cluster
	shortest = distance
	end
	end
	cluster_found.add_point point unless cluster_found.nil?
	end
	clusters.delete_if { |cluster| cluster.number_of_points == 0 }
	clusters.each { |cluster| cluster.update_center delta}
	iterations += 1
	end
	clusters
	end

view raw k_means.rb hosted with ❤ by GitHub

This takes as parameters the points themselves (as an Enumerable of Points), the number of clusters (k), the delta (defines what is a "move"), and if we want to plot each iteration. (Especially interesting while debugging the algorithm... :-)) No real surprises here, I guess. I chose to delete clusters with no points. I'm not entirely sure if that is what's in the original algorithm. (I guess by now I should duckduckgo it. :-))

Then I created a very silly little script that generates 10000 random points in 2D and (tries to) put them into 7 clusters, gnuplot'ing each iteration, to test the algorithm.

	# parameters
	k = 7
	number_of_points = 10000
	delta = 0.001
	plot_on = false
	require './k-means.rb'
	points = []
	number_of_points.times do
	points << Point2D.new(10*rand()-5.0, 10*rand()-5.0)
	end
	clusters = k_means points, k, delta, plot_on
	clusters.each_with_index do |cluster, index|
	puts "#{ cluster.center.to_s }\t#{ index }"
	end
	cluster_plot_2D clusters, "#{k}-final" if plot_on

view raw k-means-2d.rb hosted with ❤ by GitHub

Voila, the results are in:

montage -geometry "320"x"240" -tile 7x7 outfiles/clusters-7* clusters-all.png

Or maybe slightly more interesting, as an animated GIF (with sound?):

 convert outfiles/clusters-7* clusters-all.gif

Epiolog
I guess a variant of this could serve as a solution to the partitioning problem if we just map the categorical variables to the continuous space (using the weights?) and sample one element from each resulting (population/number of groups + 1) cluster (where the cluster size would be locked - Lloyd's algorithm?) to form groups.

Moral
It is useful to (try to) implement an algorithm you think you know from scratch from time time to see if you really do know all the details in the implementation of it... Especially in new languages.

Update: Clone this github if you for any reason would want to play around with this yourself: https://github.com/mortenjohs/k-means...