I have a data file with thousands of records described by 300 variables. These records belong to 12 different classes. I want to know which classes are the most similar (based on the 300 variables) so that I can group them (I want to make 3 or 4 classes out of the 12 original classes).

I think a good approach could be to calculate the distance between the classes, based on the 300 variables, and then group the classes based on the shortest distance. But I don't know how I could do that practically and where to start. I am using Python or R.

Any ideas?

  • $\begingroup$ Are all your variables numerical ? If so, are they categorical, ordinal, or continuous ? $\endgroup$ – ArnoV Jan 28 at 9:32
  • $\begingroup$ All variables are numerical and continuous (mid-infrared wavelengths values) $\endgroup$ – Falco Jan 28 at 9:40
  • 1
    $\begingroup$ What you want is a similarity between classes. You can compute similarity between classes using the centroids, but that's maybe too approximate. Have you tried k-means with 3 or 4 classes? Another worthy option is the notion of distribution discrepancy: MMD, Wasserstein etc. The Wasserstein metric is very simple to compute if your data is approximately Gaussian ! $\endgroup$ – ArnoV Jan 29 at 13:33
  • $\begingroup$ Also if you wish to compare classes up to a translation (class mean) and rotation (class covariance) you can use Procrustes analysis. $\endgroup$ – ArnoV Jan 29 at 13:35
  • 1
    $\begingroup$ As your question is fairly general I've provided code and concepts for two methods. You can also use k-means, DBSCAN , ... $\endgroup$ – ArnoV Jan 30 at 15:43

My answer supposes you have a basic knowledge of python3, and that for each of the 12 classes you have enough samples to compute representative summary statistics (a mean or covariance on 3 or 10 points does not make much sense).

I'll use randomly generated data without much structure, you'll have to reproduce the experiments on your data. Here I assume your data is a numpy array $$\mathbf{X} = \left[ \mathrm{x}_1\; \mathrm{x}_2\;\dots\;\mathrm{x}_n\right]^T$$ of shape $(n,p)$ where $n$ is the number of samples and $p=300$ the number of variables. A row $i$ is a vector $\mathrm{x}_i$ with 300 dimensions. You also have a label array $\ell = \left[\ell_1\;\ell_2\;\dots\;\ell_n\right]$ of shape $(n,)$ which I will call $\texttt{labels}$ in the python code. I assume $\mathbf{X}$ and $\ell$ are sorted in the same order, namely the $i$-th row $\mathrm{x}_i$ has label $\ell_i$

Eyeballing based on projections

Project your data $\mathbf{X}$ into a two-dimensional object using, for example, t-SNE. This can be done via scikit-learn using the following snippet:

import matplotlib.pyplot as plt
from sklearn.manifold import TSNE

# reduce using t-SNE
tsne = TSNE(n_components=2, random_state=0)
projs = tsne.fit_transform(X)

# plot
colors = [plt.cm.tab20.colors[i] for i in labels]
plt.scatter(projs[:,0], projs[:,1], c=colors)

In my case, the result is not very helpful as my data is quite random, here it is anyway:


If you spot classes, say $C_i$ and $C_j$ which are grouped together in the projected plot, you may assume the datapoints belong to a common class $\tilde{C}$. In my case, merging red-orange-green and brown-grey-purple seems interesting.

Basic: Comparing centroids


A basic approach would be to compute the centroids $\hat{\mu}_k$ of each class $C_k$: $$ \hat{\mu}_k = \dfrac{1}{\lvert C_k\rvert}\sum_{i\in C_k}\mathrm{x}_i$$ where $\lvert C_k\rvert$ is the number of samples in class $k$.

We then group together classes with similar centroids, namely you set a threshold $\delta > 0$ and decide to merge classes $C_i$ and $C_j$ together if $ \lVert{\hat{\mu}_i - \hat{\mu}_j}\rVert < \delta$, with $\lVert\cdot\rVert$ being either the $l_1, l_2$ or $l_\infty$ norms.


The centroids can be computed using numpy masks, for example

centroids = []
ulabels = np.unique(labels)

for l in ulabels:
    centroid = X[labels == l].mean(axis=0)

Then a non optimal way to compute the centroid similarities $\lVert\hat{\mu}_i - \hat{\mu}_j\rVert$ would be

distances = np.array([[np.linalg.norm(mu_i - mu_j) for mu_i in centroids] for mu_j in centroids])

This gives you a similarity matrix, which you can plot using e.g. seaborn's heatmap

import seaborn as sns

sns.heatmap(distances, linewidth=2, alpha=0.8, cmap="YlGnBu" )

In my case, this gives


And you can set delta in an ad-hoc fashion (in my case $\delta = 0.5$) to get some results. In my case, it recommends to merge classes $2,4,5,6,7$ together and classes $9,10,11$ together, leaving only 6 classes (as $0,1,3,8$ aren't similar to any other class).

Wasserstein similarity for gaussians


If your data for each class $C_i$ is approximately gaussian $\mathcal{N}(\mu_i, \Sigma_i)$, the Wasserstein-2 similarity between the two classes $C_i, C_j$ is defined as

$$ W_2(C_i,C_j) = \lVert\mu_i - \mu_j \rVert_2^2+ \mathrm{Tr}\left[\Sigma_i + \Sigma_j - 2\left(\Sigma_j^{1/2}\Sigma_i\Sigma_j^{{1/2}}\right)^{1/2}\,\right]$$

Where $\mathrm{Tr}[\cdot]$ is the trace of the matrix. This expression simplifies to the centroid similarity if you assume the variables have the same diagonal covariance $\Sigma_i = \Sigma_j = \sigma^2\mathbf{I}$. Therefore our basic centroid method makes important assumptions.

We can make the model slightly more powerful by assuming the covariances are diagonal $\sigma_i^2\mathbf{I}$, the trace term then becomes $p(\sigma_i^2 + \sigma_j^2 - 2\sigma_i\sigma_j)$ where $p$ is the number of variables (here 300). This is a sort of "variance dot product" which measures the similarity of the two variances $\sigma_i, \sigma_j$.

We can again make the model more powerful by estimating the covariance from the data, using $$ \widehat{\Sigma}_k = \dfrac{1}{\lvert C_k \rvert -1}\sum_{i\in C_k} (\mathrm{x}_i - \hat{\mu}_k)(\mathrm{x}_i - \hat{\mu}_k)^T$$ with as before our centroid $\hat{\mu}_k$ for class $k$. Once the convariances for each class are estimated, just plug them in the formula for $W_2(C_i,C_j)$ to estimate class similarity !


We know how to compute centroids already. If we assume all the covariances are diagonal of the form $\sigma_i^2\mathbf{I}$ then once we know the class, the dimensions have the same variance, therefore we can estimate $\sigma_i^2$ as

ulabels = np.unique(labels)

sigmas = []
centroids = []
for l in ulabels:
    centroid = X[labels == l].mean(axis=0)
    X_l_centered = X[labels == l] - centroid
    sigma = (X_l_centered ** 2).sum() / (X_l_centered.shape[0] - 1)

And then the wassersteins using some helper function

def wass_gauss(mu1,mu2, sq_sig1, sq_sig2):
    return np.linalg.norm(mu1-mu2)**2 + len(mu1)*(sq_sig1 + sq_sig2 - 2*np.sqrt(sq_sig1)*np.sqrt(sq_sig2))

In a list comprehension like before:

wassgauss_dists = np.array([[wass_gauss(mu_i,mu_j, sq_sig_i, sq_sig_j)
                                 for mu_i,sq_sig_i in zip(centroids, sq_sigmas)]
                                    for mu_j,sq_sig_j in zip(centroids, sq_sigmas)])

You can plot a heatmap again, set a threshold $\delta$ again, etc. My results for this wasserstein gaussian simplified case is (by setting $\delta = 50$) that one should merge $0,11$ together and $4,5$ together.

I'll leave you the more general case which isn't that complex.

Wasserstein and MMD

If you assume nothing about your data, you can still keep the idea of similarity between classes $C_i, C_j$ and use more powerful/general tools. A good package is the geomloss package, which makes things very easy! If you need help with this, post a new question and I'll gladly answer.

Clustering algorithms

You can also use well known clustering algorithms such as k-means and DBSCAN. They both have an implementation within scikit-learn.

  • $\begingroup$ Amazing...this is what I was looking for, thank you! $\endgroup$ – Falco Feb 3 at 13:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.