Introduction to Machine Learning and Clustering

Credits
AstroStat Academy logo This notebook contains original work by the authors unless stated otherwise. Any external material is properly credited to its sources.
References to papers, datasets, and software are acknowledged. Original content is licensed under the GNU General Public License v3.0 (GNU GPLv3).

Introduction to Machine Learning and Clustering #

In this session we provide a short introduction to Machine Learning, and one of its applications (Clustering).
The goals are:

  • to get a grasp of what machine learning is

  • to show the basic concept of data exploration and interpretation

  • to use clustering algorithms in a particular problem.

Setup#

from collections import defaultdict
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import colors
import seaborn as sns

from IPython.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))

import requests

url = "https://raw.githubusercontent.com/AstroStat-Academy/assets-public/main/colab/clone_and_cd_colab.py"
colab = requests.get(url).text
exec(colab)

url = "https://raw.githubusercontent.com/AstroStat-Academy/assets-public/main/styles/plot_style.py"
style = requests.get(url).text
exec(style)

import warnings
warnings.filterwarnings(action='ignore')

# some handy functions
def add_plot_info():
    """Just add some information for the plot"""
    plt.title('Clusters')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.grid(True)

Running locally, skipping Colab setup.
Imported matplotlib.
Imported seaborn.
Plotting style set.

A bit of introduction to ML#

dive_to_python_ml.jpg


Figure 1. Dive into Python and machine learning! (Source: pinterest random image)

What is machine learning ?#

Arthur Samuel (1959):

“Field of study that gives computers the ability to learn without being explicitly programmed.”

How did we get here?#

ml_history.png


Figure 2. Timeline from Artificial Intelligence, to Machine Learning, and Deep learning.
(Credit: Nvidia blog - (an interesting read!)

ml_keymoments.png


Figure 3. Key moments in AI.
(Credit: Sai Venkat Javvadi (2025)

and … why now?#

  • More powerful, abundant, and cheap computation (CPUs/GPUs).

  • Growing data sets.

  • Advancements in underlying algorithms and implementation.

This is true for both everyday and Astronomy applications !

ML branches#

ml_bracnhes.png


Figure 4. Branches and example applications of Machine Learning
(Credit: CogHub)

  • Supervised: Labelled data where the algorithms learn to predict the output from the input data.

  • Unsupervised: Non-labelled data where the algorithms learn to identify structures from the input data.

  • Semi-supervised: Some labelled data - most is not - and a mixture of supervised and unsupervised techniques can be implemented.

Unsupervised approaches in particular:#

  • Clustering: discover groupings or/and structures in the data, i.e. concentrations of datapoints or overdensities (e.g. the locations of galaxies in a BPT diagram)

  • Association: discover the rule(s) describing between variables or features in a dataset (e.g. customer recommendations)

  • Dimensionality reduction: tool to reduce the number of input variables or features in a dataset (the more you have the more challenging it becomes to build a predictive model - referred to as curse of dimensionality) - also usefulf for visualization purposes.

Pros#

  • They can make new discoveries, as often enough we don’t know what they’re looking for in data.

  • They do not require training, which saves (huge) time on producing labels (manual classification tasks such as spectroscopic classification).

  • It reduces the chance of human error and bias, which could occur during manual labeling processes.

  • Unlabeled data is much easier and faster to get.

Cons#

  • Output needs careful and proper interpretation:

    • the groups may not match informational classes

    • extra effort has to be made to validate the groups.

  • Less accurate predictive results, as the labels are not part of the process and the method has to learn it by itself.

  • More time is needed to fit these algorithms:

    • they need time to analyze and calculate all possibilities

    • the deal with huge datases that may increase computational complexity.

The most critical take-home points:#

  • No matter which algorithm you pick, the goal of ML is to make predictions and classifications.

  • There is no optimal algorithm, it all depends on your specific problem!

Clustering#

or…

“What can data tell us ?”

Selecting a tool - the K-means algorithm#

The K-means algorithm tries to partition a sample of N observations (with each observation being a \(d\)-dimensional vector) into \(k\) individual clusters \(C_k\).

We first need to find a metric, i.e. to define the loss function, which (in this case) is the within-cluster sum-of-squares of the observations:

\[ Loss = \sum_{k=1}^{K} \sum_{i \in C_k} \|x_i-\mu_k\|^2 \]

where \(\mu_k=\frac{1}{N_k}\sum_{i \in C_k} x_i\) is the mean/centroid of the \(N_k\) points included in each of the \(C_k\) clusters.

The solution comes from minimizing the above function, i.e.:

\[ \min_{\mu_k} \left( \sum_{k=1}^{K} \sum_{i \in C_k} \|x_i-\mu_k\|^2 \right) \]

Steps:

  1. Initiate algorithm by selecting \(k\) means
    e.g. select randomly \(k\) observations as initial means - see also Wiki:K-means initialization

  2. Assign each observation to the nearest cluster

  3. Calculate the new mean value for each cluster \(C_k\) according to the new observations assiged

  4. Repeat steps 2 and 3 up to the point that there are no updates in the assigments to the clusters.

kmeans.gif


Figure 5. Evolution of K-means centroids through iterations.
(Credit: Clustering with Scikit with GIFs, by David Sheehan)

A globally optimal minimum is not guaranteed (might converge to a local minimum). This is highly dependent on the initialization of the centroids. This is why, in practice, K-means is run multiple times with different starting values selecting the result with the lowest sum-of-squares error. To improve on that we can initially select centroids that are generally distant from each other (sklearn implementation by using init='k-means++' parameter). For more see Grouping data points with k-means clustering, by Jeremy Jordan.

Complexity

\(O(knT)\), where k, n and T are the number of clusters, samples and iterations, respectively.

Pros

  • Simple and intuitive

Cons

  • The number of clusters (K) must be provided (or cross-validated)

  • There is an inherit assumption of isotropic clusters (i.e. not well fitted for elongated clusters, or manifolds with irregular shapes)

  • Inertia is not a normalized metric: lower values are better, but as the dimensions increase so does the inertia

An alternative for faster implementation - Mini Batch K-means

For faster computations the sklearn offers the Mini Batch K-means method which simply breaks the initial set of observations/data points to smaller randomly selected subsamples.

For each subsample in the mini batch the assigned centroid is updated by taking into account the average of that subsample and all previous subsamples assigned to that centroid. This is repeated until the predefined number of iterations is reached. Its results are generally only slightly worse then the standard algorithm.

Example use of K-means#

In the following example we artificially create two groups of data and we use the K-means algorithm.

Creating some data

# Set random seed for reproducibility
np.random.seed(2026)

# select number of samples/points per cluster
num_samples = 50

# Generate cluster 1 with normal distribution
cluster1_mean = [3, 3]
cluster1_cov = [[2, 0], [0, 2]]
cluster1_samples = np.random.multivariate_normal(cluster1_mean, 
                                                 cluster1_cov, 
                                                 num_samples)

# Generate cluster 2 with normal distribution
cluster2_mean = [8, 8]
cluster2_cov = [[2, 0], [ 0, 2]]
cluster2_samples = np.random.multivariate_normal(cluster2_mean, cluster2_cov, num_samples)

plt.plot(cluster1_samples[:,0],cluster1_samples[:,1], 
         'o', color='tab:blue', label='Cluster 1')

plt.plot(cluster2_samples[:,0],cluster2_samples[:,1], 
         'o', color='tab:orange', label='Cluster 2')

add_plot_info()
../_images/fbd939d006ab865f5062a3ef3ca063b28a02390aac774fff39cba1cf8c3bc5f3.png

The above plot is the result of the mock data creation. In reality what we would see is the following plot (i.e. only the data without any idea of the original group).

X = np.concatenate((cluster1_samples,cluster2_samples))

plt.plot(X[:,0], X[:,1] , '.k', ms=10) 

add_plot_info()
../_images/9e0459d42735df6ac765f66e5df38c6cf8caf9ca4e43578508e807963ba5a40e.png

Implementation

Throughout the machine-learning part we will be using the scikit-learn module developed on top of NumPy, SciPy, and matplotlib, offering efficient tools on machine-learning applications, data mining, and data analysis.

We will use sklearn.cluster.KMeans, for which the most important (currently) option is the number of groups (we guess). A valid option is of course 2.

from sklearn.cluster import KMeans

# number of clusters
Clusters_kmeans = 2

# preparing the model with this number of clusters
kmeans_model = KMeans(n_clusters=Clusters_kmeans)
# fitting the model to the data
kmeans_model.fit( X )

KMeans(n_clusters=2)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Check how we apply the algorithm:

  1. Selecting the parameters (e.g. number of clusters)

  2. Run the algorithm using the .fit() method (model learns from data)

# returning the cluster identified
print("Cluster centers:")
print(kmeans_model.cluster_centers_)

# comparing with the real ones
print()
print("Real centers:")
print(cluster1_mean, cluster2_mean)

Cluster centers:
[[3.08491726 3.10191329]
 [8.21070614 8.33597448]]

Real centers:
[3, 3] [8, 8]

The algorithm returns the class assigned to each data point:

print(kmeans_model.labels_)

[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]

# coordinates of the clusters for plotting
cc_x = kmeans_model.cluster_centers_[:,0]
cc_y = kmeans_model.cluster_centers_[:,1]

# plotting the centers
plt.plot(cc_x, cc_y, 'k+', ms=100)

# plotting the datapoints color-coded according 
# to the cluster they have been assigned to
new_map = matplotlib.cm.gray.from_list('clustering', ('blue', 'red'), N=Clusters_kmeans)
scat = plt.scatter( X[:,0], X[:,1], c=kmeans_model.labels_, edgecolors='face', cmap=new_map)
cb = plt.colorbar(scat, ticks=range(0,Clusters_kmeans+1,1))   # number of clusters
cb.set_ticklabels(range(1,Clusters_kmeans+2,1))
cb.set_label('Cluster Label')

add_plot_info()

plt.show()
../_images/9daf54b858b3c60bcc872166cc5f676fad78a2a0fe28bb0c7b3bbd75135ba30f.png

Exercise [10 min]

Objective: Play around with the number of clusters as an input to K-means and explore the effect.

Task: Use the following function (cluster_generator) to create some clusters. It does it automatically for you so that you do not know how many clusters it creates (well… you can find out but you shouldn’t!). It returns an array of positions for all samples from all clusters.

  • Plot the samples and guess the number of clusters.

  • Put this into K-means and run the algorithm.

  • Discuss with your partner what you notice. You can run it a few times.

def cluster_generator(num_samples=50, verbose=False, seed=None):
    """ A function to automatically generate random 
    clusters. It returns one single dataset (points 
    with x and y) for all clusters (as a test set 
    for exercise) and the cluster centers. 
    
    num_samples: 50    , numper of points per cluster
    verbose:     False , not printing any information on 
                       the clusters
    seed:        None  , generate random values, else  
                       assign a value for reproducability. 
    """
    # for reproducability
    if seed is not None:
            np.random.seed(seed)
            
    # automatically create a number of clusters in this range
    num_clusters = np.random.randint(3, 7)

    # select number of samples/points per cluster
    num_cluster_samples = num_samples

    # Generate clusters with normal distribution
    center_clusters = []    
    data_clusters = []
    for c in range(num_clusters):        
        cluster_mean = [np.random.randint(-10, 15), np.random.randint(-10, 15)]
#        print(cluster_mean)
        cluster_cov = [[np.random.randint(1,3), 0], [0, np.random.randint(1,3)]]
#        print(cluster_cov)
        cluster_samples = np.random.multivariate_normal(cluster_mean, 
                                                         cluster_cov, 
                                                         num_cluster_samples)
        data_clusters.append(cluster_samples)
        center_clusters.append(cluster_mean)    

        if verbose==True:
            print(cluster_mean)
            plt.scatter(cluster_samples[:,0],cluster_samples[:,1], label=f'Cluster {c}')
            
    
    data_clusters = np.concatenate(data_clusters, axis=0)
#    print(center_clusters)
    return data_clusters, center_clusters


# Generating and visualizing data
X, Xc = cluster_generator()

plt.plot(X[:,0], X[:,1] , '.k', ms=10) 
add_plot_info()
../_images/7fcbfd5ead66a60f0828c9864e4184e1e4b98e650c315064043a85644d4caf8d.png 
# Fill in the blanks

from sklearn.cluster import KMeans


# number of clusters -GUESS !
Clusters_kmeans = ...

# preparing the model with this number of clusters
kmeans_model = KMeans(n_clusters=...)
# fitting the model to the data
kmeans_model.fit( ... )

# returning the cluster identified
print(f"Cluster centers from K-means ({len(...)} found):")
print(...)

# comparing with the real ones
print()
print(f"Real centers ({len(Xc)} clusters):")
for c in Xc:
    print(c)
    plt.plot(c[0],c[1], 'x', c='tab:green', ms=80)

# coordinates of the clusters for plotting
cc_x = ...[:,0]
cc_y = ...[:,1]

# plotting the centers
plt.plot(cc_x, cc_y, 'k+', ms=100)

# plotting the datapoints color-coded according 
# to the cluster they have been assigned to
new_map = matplotlib.cm.gray.from_list('clustering', ('blue', 'red'), N=Clusters_kmeans)
scat = plt.scatter( ..., ..., c=..., edgecolors='face', cmap=new_map)
cb = plt.colorbar(scat, ticks=range(0,Clusters_kmeans+1,1))   # number of clusters
cb.set_ticklabels(range(1,Clusters_kmeans+2,1))
cb.set_label('Cluster Label')

plt.text(0.75, 1.05, "x: real", color='g',
         transform=plt.gca().transAxes,
         verticalalignment='top')

plt.text(0.84, 1.05, "+: predicted", color='k',
         transform=plt.gca().transAxes,
         verticalalignment='top')


add_plot_info()

plt.show()

Questions [2 min]

How does the algorithm performs with respect to your guess ?
How dependant is on the number of clusters?

[Spoiler] (click here to expand)

When your guessed number of clusters matches the true number, K-means’ centroids (black +) typically line up closely with the real centers (green x) — cluster_generator produces fairly well-separated, roughly circular Gaussian blobs, exactly the geometry K-means handles well.

The result is highly sensitive to the number of clusters you supply: guess too low and two real clusters get merged into one; guess too high and a single real cluster gets artificially split into pieces. Since your only guide is eyeballing the scatter plot, a close-but-wrong guess (easy to make when two clusters sit near each other) directly produces one of these two failure modes.

How to find the best k value? - the “elbow” approach#

K-means needs a certain number of clusters to run. Therefore, our selection affects our interpretation (and any potential discoveries!). So, how can we improve (automatically) the selection of clusters that the K-means algorithm will search for?

One answer can be given through the “elbow method”. This method is based on the following process:

  1. Running the algorithm for a range of different clusters k (from 1 up to a maximum number).

  2. Calculating the ‘inertia’ for each k choice.

  3. Plotting the above quantity as a function of k:

  4. Identifying the ‘elbow’ point in the plot, i.e., the point on the graph where the rate of decrease significantly slows down. This indicates that adding more clusters beyond this point does not substantially reduce inertia, and thus, we do not gain much.

  5. The k value corresponding to this point is considered the optimal choice.

☛ This method does not always provide a clear answer, especially if the dataset does not exhibit well-defined clusters or if they have different sizes and densities.

Note: we are using the latest X dataset (from the previous set). You may see different results from others!

from sklearn import metrics

inert = []
xclusters = range(2,14)
for i in xclusters:
    print(f'-- working with {i} clusters...')
    kmeans_model = KMeans(n_clusters=i, random_state=0, n_init='auto')
    kmeans_model.fit(X)
    inert.append(kmeans_model.inertia_)

-- working with 2 clusters...
-- working with 3 clusters...
-- working with 4 clusters...
-- working with 5 clusters...
-- working with 6 clusters...
-- working with 7 clusters...
-- working with 8 clusters...
-- working with 9 clusters...
-- working with 10 clusters...
-- working with 11 clusters...
-- working with 12 clusters...
-- working with 13 clusters...

plt.plot( xclusters, inert)
plt.xlabel('Number of clusters')
plt.ylabel('Inertia (loss) ')
plt.show()
../_images/15d32d2981dc520371a1a162168e81958f4dfe206373cb32918af8a942d77d30.png

Let’s plug in the number of clusters corresponding to the “elbow” found:

from sklearn.cluster import KMeans


# number of clusters - from "elbow" plot
Clusters_kmeans = 3

# preparing the model with this number of clusters
kmeans_model = KMeans(n_clusters=Clusters_kmeans)
# fitting the model to the data
kmeans_model.fit( X )

# returning the cluster identified
print(f"Cluster centers from K-means ({len(kmeans_model.cluster_centers_)} found):")
print(kmeans_model.cluster_centers_)

# comparing with the real ones
print()
print(f"Real centers ({len(Xc)} clusters):")
for c in Xc:
    print(c)
    plt.plot(c[0],c[1], 'x', c='tab:green', ms=80)

# coordinates of the clusters for plotting
cc_x = kmeans_model.cluster_centers_[:,0]
cc_y = kmeans_model.cluster_centers_[:,1]

# plotting the centers
plt.plot(cc_x, cc_y, 'k+', ms=100)

# plotting the datapoints color-coded according 
# to the cluster they have been assigned to
new_map = matplotlib.cm.gray.from_list('clustering', ('blue', 'red'), N=Clusters_kmeans)
scat2 = plt.scatter( X[:,0], X[:,1], c=kmeans_model.labels_, edgecolors='face', cmap=new_map)
cb = plt.colorbar(scat2, ticks=range(0,Clusters_kmeans+1,1))   # number of clusters
cb.set_ticklabels(range(1,Clusters_kmeans+2,1))
cb.set_label('Cluster Label')

plt.text(0.75, 1.05, "x: real", color='g',
         transform=plt.gca().transAxes,
         verticalalignment='top')

plt.text(0.84, 1.05, "+: predicted", color='k',
         transform=plt.gca().transAxes,
         verticalalignment='top')

add_plot_info()

plt.show()

Cluster centers from K-means (3 found):
[[-10.02707225   5.18145856]
 [ 10.80310614  -3.87195188]
 [  3.04618638   3.19399703]]

Real centers (3 clusters):
[-10, 5]
[11, -4]
[3, 3]
../_images/92a255798456ac22ece10afe2ca39d938fed6610d2e05910ddffa08deaf88743.png

Other metrics to automate the cluster number#

Now, given that we do not know anything apriori for the numbers of clusters, there are some metrics that can used to determine the performance in each case.

Silhouette Score#

Using the distances between the points in the same cluster and theirs with all other points in the next nearer cluster:

\[ s = \frac{b - a}{\max(a,b)} \]

where \(a\) corresponds to the mean distance between a sample and all other samples within the same cluster, while \(b\) is the mean distance of a sample (in the examined cluster) with all other points in the next nearest cluster.

⤍ A higher value indicates better performance (-1…1).

Calinski-Harabasz Index#

Similar to the silhouette score, but using the ratio of the sum of dispersion of the samples within each cluster and the dispersion of the clusters in total. For a set of data \(E\) of size \(n_E\) which has been clustered into \(k\) clusters, it is defined as:

\[ s = \frac{\mathrm{tr}(B_k)}{\mathrm{tr}(W_k)} \times \frac{n_E-k}{k-1} \]

where \(\mathrm{tr}(B_k)\) is the trace of the between group dispersion matrix, and \(\mathrm{tr}(W_k)\) the trace of the within-cluster dispersion matrix, defined as:

\[ W_k = \sum_{q=1}^{K} \sum_{x \in C_q} (x - c_q)(x - c_q)^T \]
\[ B_k = \sum_{q=1}^{k} n_q(c_q-c_E)(c_q-c_E)^T \]

with \(C_q\) the set of points in cluster \(q\), \(c_q\) the center of cluster \(q\), \(c_E\) the center of \(E\), and \(n_q\) the number of points in cluster \(q\).

⤍ A higher value indicates better performance (denser and well separated clusters).

Davies-Bouldin Index#

It is an index that calculates the ‘similarity’ between clusters, and actually measures how the distance between clusters compares with the sizes of the clusters themselves.

\[ DB = \frac{1}{k} \sum_{i=1}^{k} \max_{i \neq j} R_{ij} \]

where \(R_{ij} = \frac{s_i+s_j}{d_{ij}}\) corresponds to the similarity measurement, \(s_i\) (\(s_k\)) which is the average distance between each point of cluster \(i\) (k) and the centroid of that cluster (diameter), and \(d_{ij}\) is the distance between cluster centroids \(i\) and \(j\). The form chosen for \(R_{ij}\) is a nonnegative and symmetric.

⤍ A lower value (closer to 0) indicate better performance.

☛ You can (should!) check sklearn’s clustering metrics page for more details regarding the drawbacks and advantages of each metric.

Let’s see them in action!

from sklearn.metrics import silhouette_score, calinski_harabasz_score, davies_bouldin_score

# Define range of cluster numbers to test
k_values = range(2, 10)

# Initialize lists to store scores
inertia = []
silhouette = []
calinski = []
davies = []

# Fit KMeans and compute metrics for each k
for k in k_values:
    kmeans = KMeans(n_clusters=k, n_init=10, random_state=42)
    labels = kmeans.fit_predict(X)

    inertia.append(kmeans.inertia_)
    silhouette.append(silhouette_score(X, labels))
    calinski.append(calinski_harabasz_score(X, labels))
    davies.append(davies_bouldin_score(X, labels))

# Plot all metrics for visual comparison
fig, axs = plt.subplots(2, 2, figsize=(12, 10))

axs[0, 0].plot(k_values, inertia, 'bo-')
axs[0, 0].set_title('Elbow Method (Inertia)')
axs[0, 0].set_xlabel('k')
axs[0, 0].set_ylabel('Inertia')

axs[0, 1].plot(k_values, silhouette, 'go-')
axs[0, 1].set_title('Silhouette Score')
axs[0, 1].set_xlabel('k')
axs[0, 1].set_ylabel('Score')

axs[1, 0].plot(k_values, calinski, 'ro-')
axs[1, 0].set_title('Calinski-Harabasz Index')
axs[1, 0].set_xlabel('k')
axs[1, 0].set_ylabel('Score')

axs[1, 1].plot(k_values, davies, 'mo-')
axs[1, 1].set_title('Davies-Bouldin Index')
axs[1, 1].set_xlabel('k')
axs[1, 1].set_ylabel('Score')

plt.tight_layout()
plt.show()
../_images/c4435eb4d5c4c7568e54150e01bba83d31cb05e3ca8ce9164ed44b9678cd8757.png

Our pick would be 3 clusters in this case. We note the following:

  1. In this case all the metrics are consistent - but they could provide slightly different results, as they are sensitive to different cluster properties.

  2. Combining different approaches leads to robust results.

  3. But careful interpretation is required to derive non-misleading conclusions.

Selecting a tool - the DBSCAN algorithm#

The DBSCAN (Density Based Spatial Clustering of Applications with Noise) algorithm views clusters as areas of high density separated by areas of low density. Thus, clusters found by DBSCAN can be of any shape - as opposed to other algorithms (K-means for example) which assume that clusters are convex shaped.

It uses 2 parameters:

  • \(eps\) : neighborhood size

  • \(minPts\) : minimum number of points for a neighborhood to be considered dense

With a sinlge scan we can label the points as: core, border, noise. How?

  • A point \(p\) is defined core if at least \(minPts\) points are within the area defined by \(eps\) (including itself).

  • A border point is a non-core point that has at least one core point in its neighborhood

  • A noise point is neither a core nor a border point. There represent outliers in the data set

dbscan.png


Figure 6. Definition of core, border, and noise points, according to DBSCAN,
(Credit: Clustering applied to showers in the OPERA, by Alex Rogozhnikov)

Having defined these points, DBSCAN operates as follows:

Steps:

  1. The algorithm selects a core point that has not been explored.

  2. It finds the points that are within an eps distance from this point and examines them.

  3. If these points are also core points, it examines their eps-neighborhood, expanding the cluster.

  4. This expansion continues recursively: each core point found this way has its own eps-neighborhood checked in turn, pulling in its neighbors (whether core or border) into the same cluster. The cluster keeps growing until no more reachable core points are left to expand from.

  5. When these points are exhausted, it selects the next unvisited core point (in another cluster) and continues the previous process.

  6. The process ends when all points have been assigned to a cluster or identified as noise.

dbscan-smiley.gif


Figure 7. An example of how it constructs clusters (right).
(Credit: Clustering applied to showers in the OPERA, by Alex Rogozhnikov)

Complexity

\(O(n^2)\), where n is the number of points.

Pros

  • Any number of clusters

  • Clusters of varying size and shape

  • Finds and ignores outliers

Cons

  • Relatively slow

  • Extremely sensitive to parameters choice

  • In rare cases, border points move to an other cluster when DBSCAN is re-run

  • Serious troubles with clusters with varying density
    (OPTICS and HDBSCAN are variations which address this problem)

Example use of DBSCAN#

In the following example we artificially create two groups of data and we use the DBSAN algorithm.

Loading some data

# import this to create some fancy data
from sklearn.datasets import make_moons

# Generate synthetic data (moons dataset)
Xd, y = make_moons(n_samples=200, noise=0.1, random_state=2026)

num_classes = len(set(y))
new_map = matplotlib.cm.gray.from_list('clustering', ('blue', 'red'), N=num_classes)
scat2 = plt.scatter( Xd[:,0], Xd[:,1], c=y, 
                    edgecolors='face', cmap=new_map)
cb = plt.colorbar(scat2, ticks=range(1,num_classes+2,1))   # number of clusters
cb.set_ticklabels(range(1,num_classes+2,1))
cb.set_label('Cluster Label')

add_plot_info()
../_images/6dc632a5013417167bf6d375506b190dde0ad0caf54974e833dc7199fb619e91.png

What we should see if we combine into one dataset.

plt.plot(Xd[:,0], Xd[:,1], 'ok') 

add_plot_info()
../_images/ec54f3acb8dd1dfa349d542039bbc7edc3758746cb269e1267f0c991d17a6f02.png

Implementation

We will use sklearn.cluster.DBSCAN. The most important parameters of course are eps and min_samples.

from sklearn.cluster import DBSCAN

# define the model first
dbscan = DBSCAN(eps=0.08, min_samples=3)
# fit the data
dbscan.fit(Xd)

DBSCAN(eps=0.08, min_samples=3)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Notice how all algorithms are applied using the same methods.

Let’s see the number of classes returned by the algorithm:

print(set(dbscan.labels_))

{np.int64(0), np.int64(1), np.int64(2), np.int64(3), np.int64(4), np.int64(5), np.int64(6), np.int64(7), np.int64(8), np.int64(9), np.int64(10), np.int64(11), np.int64(12), np.int64(13), np.int64(14), np.int64(15), np.int64(16), np.int64(17), np.int64(18), np.int64(19), np.int64(20), np.int64(21), np.int64(22), np.int64(23), np.int64(24), np.int64(25), np.int64(-1)}

Let’s visualize the predictions now.

num_classes = len(set(dbscan.labels_))
new_map = matplotlib.cm.gray.from_list('clustering', ('blue', 'red'), N=num_classes)
scat2 = plt.scatter( Xd[:,0], Xd[:,1], c=dbscan.labels_, 
                    edgecolors='face', cmap=new_map)
cb = plt.colorbar(scat2, ticks=range(-1,num_classes,1))   # number of clusters excluding noise
cb.set_ticklabels(range(-1,num_classes,1))
cb.set_label('Cluster Label')

add_plot_info()
../_images/7fc35c34be83e5fae37d6bef05b8d93c0418aa99c60766e8ec882df3b63d50f7.png

Well…this is kind of embarrassing ! DBSCAN should work, what’s the problem?

The 3 panels below show what happens geometrically when eps is far from the right value.

Each run uses the same dataset (Xd) with min_samples = 5. Grey crosses mark points classified as noise (label = −1).

# eps sensitivity demo: same data, 3 different eps values
from sklearn.cluster import DBSCAN

fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=True)

eps_values = [0.05, 0.22, 0.60]
subtitles  = [
    'eps = 0.05  (too small)\nmost points become noise',
    'eps = 0.22  (about right)\n2 clusters recovered',
    'eps = 0.60  (too large)\neverything is 1 cluster'
]

for ax, eps_val, subtitle in zip(axes, eps_values, subtitles):
    db     = DBSCAN(eps=eps_val, min_samples=5).fit(Xd)
    labels = db.labels_

    # noise points in grey
    noise = labels == -1
    ax.scatter(Xd[noise, 0], Xd[noise, 1],
               c='lightgrey', s=15, marker='x', zorder=1)

    # one colour per cluster
    for k in sorted(set(labels)):
        if k == -1:
            continue
        mask = labels == k
        ax.scatter(Xd[mask, 0], Xd[mask, 1], s=15, alpha=0.7, zorder=2)

    n_cl = len(set(labels) - {-1})
    n_no = int(noise.sum())
    ax.set_title(f'{subtitle}\n({n_cl} cluster(s), {n_no} noise pts)', fontsize=9)
    ax.set_xticks([])
    ax.set_yticks([])

plt.suptitle('Effect of eps on DBSCAN  (min_samples = 5)', fontsize=11)
plt.tight_layout()
plt.show()
../_images/3a05a00e04da7abeb7eaaf1dfc2f70f8991e57ee38fa41b4b7f2d67013a36d7d.png

Question [5 min]

Which are the optimal parameters?

In order to identify the two clusters play around with the eps and mins_sample.

HINT: suggested ranges 0.05 – 0.30, 2 – 7, respectively

NOTE: the group of points with -1 indicate noise.

[Spoiler] (click here to expand)

A working pair of parameters for this dataset:

  • eps = 0.22

  • min_samples = 5

Comparing the two algorithms#

Let’s now check if K-means is able to recover the clusters. It is kind of obvious what should be the number of clusters…

# number of clusters
Clusters_kmeans = 2

# preparing the model with this number of clusters
kmeans_model = KMeans(n_clusters=Clusters_kmeans)
# fitting the model to the data
kmeans_model.fit( Xd )

# coordinates of the clusters for plotting
cc_x = kmeans_model.cluster_centers_[:,0]
cc_y = kmeans_model.cluster_centers_[:,1]

# plotting the centers
plt.plot(cc_x, cc_y, 'k+', ms=100)

# plotting the datapoints color-coded according 
# to the cluster they have been assigned to
new_map = matplotlib.cm.gray.from_list('clustering', ('blue', 'red'), N=Clusters_kmeans)
scat2 = plt.scatter( Xd[:,0], Xd[:,1], c=kmeans_model.labels_, edgecolors='face', cmap=new_map)
cb = plt.colorbar(scat2, ticks=range(0,Clusters_kmeans+1,1))   # number of clusters
cb.set_ticklabels(range(1,Clusters_kmeans+2,1))
cb.set_label('Cluster Label')

add_plot_info()

plt.show()
../_images/9166b025df402df92f76b60284cbd3bc6b23b016e8048f45b84a7fe6a7fc9163.png

Question [2 min]

What would you comment regarding the performance of K-means instead of DBSCAN?

[Spoiler] (click here to expand)

K-means does find 2 clusters, but they are geometrically wrong.

K-means assumes clusters are roughly spherical (isotropic), so it divides the crescent-shaped moons along a straight boundary rather than following their curved geometry.

DBSCAN has no such assumption: it groups points by local density regardless of shape, which is why it succeeds where K-means fails on this dataset.

Application#

The following dataset originates from measurements of two spectral lines for a sample of stars. The strength of a spectral line can be measured by its equivalent width (EQW).

equivalent_width-wiki.jpg


Figure 8. Defining EQW: the width of a line with intensity equal to the local continuum and total flux equal to that of the line.
(Credit: Wikipedia: Equivalent Width, by Szdori )

Therefore, the larger the EQW the stronger the line is. The presence of spectral lines depends on the temperature of the stellar sources. Because of this, we see a developmet of different spectral lines as we move from the hottest to the cooler stars. This corresponds to moving from earlier spectral types (O-type stars; 50-25kK) to later ones (M-type; 3.5-2.5kK) see the Morgan-Keenan spectral classification scheme].

Exercise [15 min]

Objective: Identify the clusters of spectral types for hot massive stars based on the ratio of spectral lines.

Taks: To do this follow the automated approach for K-means. In particular:

  • Run the automated estimate of the identified clusters we run previously for K-means.

  • Check all metrics and decide which one to use.

  • Run the K-means algorithm with the selected number of clusters and visualize the result.

  • Think over your results (what could they mean?) and cross-check with the real data.

def flospecConv(arg):
    """
    Function to convert from spectral types to 
    float numbers (e.g. B0,O9.5 to 20.0,19.5)
    and backwards.
    """ 
    try:
        float(arg)
        if str(arg)[0]=='1':
            sp = 'O'
        elif str(arg)[0]=='2':
            sp = 'B'
        elif str(arg)[0]=='3':
            sp = 'A'
        else:
            sys.exit(' ! ERROR: more than O/B stars! Adjust conversion function.')
        new_arg = sp+str(arg)[1:]
    except ValueError:
        if arg[0]=='O' or arg[0]=='o':
            fl = '1'
        elif arg[0]=='B' or arg[0]=='b':
            fl = '2'
        elif arg[0]=='A' or arg[0]=='a':
            fl = '3'
        else:
            sys.exit(' ! ERROR: Check input! If more than O/B stars adjust conversion function.')
        new_arg = float(arg.replace(arg[0],fl))

    return new_arg

# Reading the data file and selecting lines
# when selecting lines do it in pairs, and with shorter line first

PATH_data = "data/stellar_types.dat"

sellines = ['HeII/4200', 'HeI/4471'] # in order of wavelength
#sellines = ['HeI/4471', 'MgII/4481'] # in order of wavelength

stars=defaultdict(list)
with open(PATH_data,'r') as inf:
    for line in inf:
        cols = line.split()
        objt = cols[0]
        spln = cols[1]
        if spln in sellines: 
#            print(splin)
            eqw = cols[2]
            stars[objt].append(eqw)             # the shorter line is appended first!
            
#print(stars)

# Creating data structures:

sptype, flosptype, eqwA, eqwB = [], [], [], []
for s in stars.keys():
    sptype.append(s.split('-')[0])
    flosptype.append(flospecConv(s.split('-')[0]))
    eqwA.append(float(stars[s][0]))                 # A is the shorter line
    eqwB.append(float(stars[s][1]))                 # B is the other line (obviously!)
    
# > Organizing data in an analysis-ready fashion:
X = np.column_stack((eqwA,eqwB))

print('Sample shape:')
print("___________________________________")
print('  X  | ' + str(X.shape))
print('     | ' + str(X.shape[0]) + ' samples x ' + str(X.shape[1]) + ' diagnostics' )

Sample shape:
___________________________________
  X  | (697, 2)
     | 697 samples x 2 diagnostics

# Visualize data

scat = plt.scatter(eqwA, eqwB, c='k')

# automatic labelling from the lines 
plt.xlabel(r""+f"EQW of { sellines[0].split('/')[0]} $\lambda${sellines[0].split('/')[1]}")
plt.ylabel(r""+f"EQW of { sellines[1].split('/')[0]} $\lambda${sellines[1].split('/')[1]}")

# simpler approach
# plt.xlabel(r"EQW of line HeII $\lambda$4200 ")
# plt.ylabel(r"EQW of line HeI $\lambda$4471 ")

plt.show()
../_images/16f7c0ef68af515a51a167cf28eba5e144b6ddb0e8bda138aa1a65c19f129d71.png 
# Enter code here (automated detection 
# of clusters for K-means)

# Enter code here (application of K-means with 
# selected number of clusters and visualization)

Clustering algorithms overview#

But how a good clustering workflow hould look like?

To takle effectively a clustering problem initialy you can follow thise steps:

  1. Load your data

  2. Visualize your data, different stuctures requier different algorithms

  3. Choose the algorithm you want and fit it

  4. Visualize the results and see if it is acceptable

Here we present a set of sklearn clustering algorithms with toy datasets, and then we continue by applying some of them in various astrophyiscal datasets.

This serves as a showcase of the available methods and how they compare. You can easily adapt any of these methods to the following examples or your own problems.

sklearn_clustering.png


Figure 9. A comparison of the clustering algorithms in scikit-learn.
(Credit: Overview of clustering methods, scikit-learn API)

Question [2 min]

Explore the results above - what differences do you notice?

[Spoiler] (click here to expand)

  • There is no single best algorithm.

  • Not all algorithms identify the same number of clusters.

  • Some algorithms detect arbitrary cluster shapes better than others.

  • Some algorithms are faster.

  • Intuitive clustering may not apply to very high-dimensional data.

###EOF