Self-Organizing Map (SOM): Complete Guide to Kohonen Self-Organizing Maps

Summary

Self-Organizing Maps (SOM) are one of the most elegant unsupervised learning algorithms ever designed. Introduced by Finnish professor Teuvo Kohonen in 1982, these artificial neural networks possess a remarkable property: they project high-dimensional data onto a two-dimensional map while preserving the topology of the input space. In other words, similar samples in the original space are found close to each other on the map. This characteristic makes SOM a powerful tool for multidimensional data visualization, clustering, and exploring hidden structures in complex datasets. Unlike traditional clustering methods like K-Means, a SOM does not simply produce isolated groups: it builds an organized spatial representation where geographic proximity on the map reflects similarity between data points.

Mathematical Principle of the Algorithm

The operation of a Self-Organizing Map relies on a competitive learning mechanism structured around three fundamental steps. Each step is essential to understanding how the map manages to organize itself.

Step 0: Initialization

We have a grid of neurons (usually rectangular or hexagonal). Each neuron i has a weight vector wᵢ of the same dimension as the input data. These weights are initialized randomly, meaning that at the start, the map is completely disorganized and reflects no structure of the data.

Step 1: Finding the Best Matching Unit (BMU)

For each training sample x, we search for the neuron whose weight is closest to x according to a distance metric, typically the Euclidean distance:

$$c = \arg\min_i ||\mathbf{x} – \mathbf{w}_i||$$

The neuron c thus identified is called the Best Matching Unit (BMU), i.e., the best matching unit. It is the neuron whose weight vector most closely resembles the current sample. This competition between neurons is at the heart of SOM learning.

Step 2: Weight Update

Once the BMU is identified, we update not only its weights but also those of its neighbors on the grid. The update rule is as follows:

$$\mathbf{w}_i(t+1) = \mathbf{w}_i(t) + \alpha(t) \cdot h(c, i, t) \cdot (\mathbf{x} – \mathbf{w}_i(t))$$

where:

• α(t) is the learning rate at iteration t, which decreases over time.
• h(c, i, t) is the Gaussian neighborhood function that measures the influence of the BMU c on the neighboring neuron i.

Step 3: Neighborhood Function

The neighborhood function determines how a neuron i is affected by the presence of the BMU c:

$$h(c, i, t) = \exp\left(-\frac{||\mathbf{r}_c – \mathbf{r}_i||^2}{2 \cdot \sigma(t)^2}\right)$$

where r_c and r_i are the respective positions of neurons c and i on the two-dimensional grid. The Gaussian function ensures that neurons close to the BMU are strongly updated, while distant neurons are updated very little.

Parameter Decay

To ensure algorithm convergence, two parameters decay exponentially over the iterations:

$$\alpha(t) = \alpha_0 \cdot \exp(-t / \text{max_iterations})$$

$$\sigma(t) = \sigma_0 \cdot \exp(-t / \text{max_iterations})$$

• α₀ is the initial learning rate.
• σ₀ is the initial neighborhood radius.

This decay is crucial: in the early iterations, a large neighborhood radius allows the map to organize globally. Progressively, the radius shrinks and learning becomes more local, refining the details of the map. It is this process that allows the SOM to preserve topology: neurons that are close on the grid end up with similar weights, thus reflecting the structure of the input data.

Intuition: The Archipelago of Ecosystems

Imagine an archipelago of islands arranged on a regular grid in the middle of the ocean. Each island represents a neuron of the self-organizing map. At the beginning, each island has a completely random habitat: some are desert, others are covered in tropical forests, some are snowy, others are marshlands. No order is visible from the sky.

Now imagine that animals representing your data arrive progressively, one by one. Each animal seeks the island whose habitat best suits it. A polar bear will naturally choose the snowy island, a toucan will prefer the tropical island. This elected island is the Best Matching Unit (BMU) for that animal.

But here is the magic of the algorithm: when an island receives a visitor, it adapts its habitat to even better match the animal, and its neighboring islands do the same, but to a lesser extent. The island itself changes a lot, its direct neighbors change a little, and the more distant islands barely change at all. This is the neighborhood function in action.

Over time, with thousands of animals visiting the archipelago, something remarkable happens: islands close to each other end up with similar habitats. The archipelago organizes itself spontaneously! Snowy islands cluster in one corner, tropical forests in another, deserts in a third. The map becomes an organized representation of the ecosystems present in the data. Similar animals always end up in the same region of the archipelago. Topology is preserved: if you walk across the map, you traverse ecosystems that change gradually, without abrupt jumps.

This is exactly what a Self-Organizing Map does with your data: it builds a map where spatial proximity reflects semantic similarity.

Python Implementation

From-Scratch Implementation of a 2D SOM

Here is a complete implementation of a Self-Organizing Map from scratch, with no dependencies other than NumPy and Matplotlib. This version uses a rectangular grid and a Gaussian neighborhood function.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.preprocessing import MinMaxScaler

class SelfOrganizingMap:
    """
    From-scratch implementation of a Self-Organizing Map (SOM)
    by Kohonen with a 2D rectangular grid.
    """

    def __init__(self, map_width, map_height, input_dim,
                 sigma0=None, learning_rate0=0.5, n_iterations=10000,
                 seed=42):
        self.map_width = map_width
        self.map_height = map_height
        self.input_dim = input_dim
        self.n_iterations = n_iterations
        self.learning_rate0 = learning_rate0
        # Initial radius: max of half the map dimensions
        self.sigma0 = sigma0 if sigma0 else max(map_width, map_height) / 2.0

        rng = np.random.RandomState(seed)
        # Random weight initialization
        self.weights = rng.rand(map_width, map_height, input_dim)

        # Position grid (2D coordinates)
        self.grid_x = np.arange(map_width)
        self.grid_y = np.arange(map_height)
        self.pos_x, self.pos_y = np.meshgrid(self.grid_x, self.grid_y)
        self.positions = np.column_stack([self.pos_x.ravel(), self.pos_y.ravel()])

    def _compute_distances(self, x):
        """Computes the Euclidean distance between x and all weights."""
        diff = self.weights - x[np.newaxis, np.newaxis, :]
        return np.sqrt(np.sum(diff ** 2, axis=2))

    def find_bmu(self, x):
        """Finds the Best Matching Unit for a sample x."""
        distances = self._compute_distances(x)
        idx = np.argmin(distances)
        bmu_x, bmu_y = np.unravel_index(idx, (self.map_width, self.map_height))
        return bmu_x, bmu_y

    def _gaussian_neighborhood(self, bmu_x, bmu_y, sigma):
        """Computes the Gaussian neighborhood function for all neurons."""
        dist_sq = ((self.pos_x - bmu_x) ** 2 + (self.pos_y - bmu_y) ** 2)
        neighborhood = np.exp(-dist_sq / (2 * sigma ** 2))
        return neighborhood

    def fit(self, X, verbose=True):
        """Trains the SOM on data X."""
        for t in range(self.n_iterations):
            # Random sample selection
            idx = np.random.randint(0, len(X))
            x = X[idx]

            # Parameter decay
            alpha = self.learning_rate0 * np.exp(-t / self.n_iterations)
            sigma = self.sigma0 * np.exp(-t / self.n_iterations)

            # Find the BMU
            bmu_x, bmu_y = self.find_bmu(x)

            # Neighborhood function
            neighborhood = self._gaussian_neighborhood(bmu_x, bmu_y, sigma)
            neighborhood = neighborhood[:, :, np.newaxis]  # Broadcasting

            # Weight update
            self.weights += alpha * neighborhood * (x - self.weights)

            if verbose and t % 1000 == 0:
                print(f"Iteration {t}/{self.n_iterations} | "
                      f"α={alpha:.4f} | σ={sigma:.4f}")

        return self

    def transform(self, X):
        """Projects data X onto the map (returns BMU coordinates)."""
        bmu_coords = []
        for x in X:
            bx, by = self.find_bmu(x)
            bmu_coords.append((bx, by))
        return np.array(bmu_coords)

    def u_matrix(self):
        """Computes the U-Matrix (map of average distances between neighbors)."""
        u = np.zeros((self.map_width, self.map_height))
        for i in range(self.map_width):
            for j in range(self.map_height):
                neighbors = []
                for di, dj in [(-1,0),(1,0),(0,-1),(0,1)]:
                    ni, nj = i+di, j+dj
                    if 0 <= ni < self.map_width and 0 <= nj < self.map_height:
                        dist = np.linalg.norm(self.weights[i,j] - self.weights[ni,nj])
                        neighbors.append(dist)
                if neighbors:
                    u[i,j] = np.mean(neighbors)
                else:
                    u[i,j] = 0
        return u

Application on the Iris Dataset

Let’s now illustrate the use of our SOM on the famous Fisher’s Iris dataset, a machine learning classic containing 150 samples of three iris species.

# Load and normalize data
iris = load_iris()
X = iris.data
y = iris.target
scaler = MinMaxScaler()
X_scaled = scaler.fit_transform(X)

# Create and train the SOM
som = SelfOrganizingMap(
    map_width=10,
    map_height=10,
    input_dim=4,
    sigma0=4.0,
    learning_rate0=0.5,
    n_iterations=50000,
    seed=42
)
som.fit(X_scaled, verbose=True)

# Visualize the U-Matrix
u_mat = som.u_matrix()
plt.figure(figsize=(8, 6))
plt.imshow(u_mat.T, cmap='viridis')
plt.colorbar(label="Average distance between neighboring neurons")
plt.title("U-Matrix of the SOM applied to the Iris dataset")
plt.xlabel("X coordinate on the grid")
plt.ylabel("Y coordinate on the grid")
plt.tight_layout()
plt.savefig("som_iris_umatrix.png", dpi=150)
plt.show()

# Mark the BMU of each class on the grid
bmu_coords = som.transform(X_scaled)
plt.figure(figsize=(8, 6))
colors = ['red', 'green', 'blue']
species = ['Setosa', 'Versicolor', 'Virginica']
for cls in range(3):
    mask = y == cls
    bx, by = bmu_coords[mask].T
    plt.scatter(bx, by, c=colors[cls], label=species[cls], alpha=0.6,
                edgecolors='black', s=100)
plt.title("SOM Projection — Iris Species on the Grid")
plt.xlabel("X coordinate")
plt.ylabel("Y coordinate")
plt.legend()
plt.gca().set_xlim(-1, som.map_width)
plt.gca().set_ylim(-1, som.map_height)
plt.tight_layout()
plt.savefig("som_iris_clusters.png", dpi=150)
plt.show()

Using MiniSom

For practical applications, the MiniSom library offers a robust and optimized implementation of SOMs. It installs easily via pip.

pip install minisom

from minisom import MiniSom

# MiniSom configuration
som_minisom = MiniSom(
    x=10, y=10,
    input_len=4,
    sigma=4.0,
    learning_rate=0.5,
    neighborhood_function='gaussian'
)

# Random initialization
som_minisom.random_weights_init(X_scaled)

# Training with decay
som_minisom.train_random(
    X_scaled,
    num_iteration=50000,
    verbose=True
)

# Visualize the weight map
plt.figure(figsize=(10, 8))
for i in range(4):
    plt.subplot(2, 2, i+1)
    plt.imshow(som_minisom.get_weights()[:, :, i].T, cmap='plasma')
    feature_names = ['Sepal length', 'Sepal width', 'Petal length', 'Petal width']
    plt.title(f"Weight Map — {feature_names[i]}")
    plt.colorbar()
plt.tight_layout()
plt.savefig("som_weight_maps.png", dpi=150)
plt.show()

RGB Color Clustering with SOM

A particularly visual example of SOM utility is color clustering: take an image, extract all its RGB pixels, and train a SOM to group them into a representative palette.

from skimage import io

# Load an image (replace with your path)
image = io.imread("photo.jpg")
pixels = image.reshape(-1, 3).astype(np.float64) / 255.0

# SOM to extract a palette of 16 colors (4x4)
color_som = SelfOrganizingMap(
    map_width=4,
    map_height=4,
    input_dim=3,
    sigma0=2.0,
    learning_rate0=0.5,
    n_iterations=20000,
    seed=123
)
color_som.fit(pixels, verbose=False)

palette = color_som.weights.reshape(4, 4, 3)

# Display the palette
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

axes[0].imshow(image)
axes[0].set_title("Original image")
axes[0].axis("off")

axes[1].imshow(palette)
axes[1].set_title("SOM Palette (4×4)")
axes[1].axis("off")

plt.tight_layout()
plt.savefig("som_color_palette.png", dpi=150)
plt.show()

# Reconstruct the image with the reduced palette
from scipy.spatial.distance import cdist
distances = cdist(pixels, color_som.weights.reshape(16, 3))
assigned = np.argmin(distances, axis=1)
reconstructed = color_som.weights.reshape(16, 3)[assigned].reshape(image.shape)

plt.imshow(reconstructed)
plt.title(f"Image reconstructed with {16} SOM colors")
plt.axis("off")
plt.savefig("som_reconstructed.png", dpi=150)
plt.show()

This application perfectly demonstrates the power of SOMs: the map has automatically extracted a palette of 16 dominant colors that preserve the chromatic distribution of the original image, while organizing these colors in a topological manner on the grid.

Key Hyperparameters

The choice of hyperparameters is critical for the quality of a Self-Organizing Map. Here are the main ones:

Practical Tips

• Map size: an empirical rule suggests approximately 5 × √N neurons, where N is the number of samples. For 1,000 data points, a 15×15 grid (225 neurons) is appropriate.
• σ₀ too small: the map fragments and does not organize globally.
• σ₀ too large: organization is too slow and local details are poorly captured.
• Learning rate too high: instability and weight oscillation.
• Learning rate too low: very slow convergence, suboptimal map.
• Number of iterations: enough iterations are needed for σ(t) and α(t) to decay to very low values, ensuring final stabilization.
• Data normalization: always normalize data before training (MinMaxScaler or StandardScaler). SOMs are sensitive to different scales between variables.

Advantages and Limitations

Advantages

1. Topological preservation: the distinctive characteristic of SOMs. Similar data are close on the map, providing an intuitive visualization of high-dimensional structures.
2. Intuitive visualization: the U-Matrix clearly reveals boundaries between clusters, where distances between neighboring neurons are high.
3. Unsupervised and assumption-free: unlike K-Means, no assumptions about cluster shape are necessary (sphericity, equal size, etc.).
4. Robustness to noise: the collective nature of learning (neighborhood update) makes SOMs less sensitive to outliers than individual methods.
5. Natural cluster detection: the map reveals the underlying structure of the data without requiring the number of clusters to be specified in advance.
6. Interpretability: weight maps can be visualized for each individual feature, providing a fine understanding of what each region of the map represents.

Limitations

1. Non-determinism: random weight initialization and sampling order make results non-reproducible without setting a random seed.
2. Delicate hyperparameter choice: map size, sigma, learning rate, and number of iterations often require experience and multiple trials.
3. Computational cost: each iteration requires computing the distance between the sample and all weights. For very large maps and lots of data, training can be long.
4. No intrinsic quality metric: unlike K-Means which minimizes intra-cluster inertia, SOMs have no single objective score to evaluate the quality of the result.
5. Difficult to evaluate quantitatively: visualization is subjective, and there are few standard metrics to compare two maps.
6. Sensitivity to normalization: non-normalized variables dominate the Euclidean distance and completely bias the result.

4 Concrete Use Cases

1. Customer Data Exploration in Marketing

SOMs are excellent for segmenting a customer base without prior assumptions about the number of segments. By projecting demographic, behavioral, and transactional data onto a 2D map, natural zones are discovered: premium customers, occasional customers, new users, etc. The topology of the map also reveals transitions between segments, identifying borderline customers who could shift from one behavior to another.

2. Document and Text Analysis

By representing each document as a TF-IDF vector or embeddings (Word2Vec, BERT), a SOM can automatically organize a document collection. Articles on similar themes cluster together on the map. A journalist could thus explore archives of thousands of articles by navigating visually on the map, identifying emerging themes and outlier articles that deserve attention.

3. Medical Diagnostic Imaging

SOMs can analyze features extracted from medical images (X-rays, MRIs, CT scans) to detect anomalies. By training a map on healthy tissues, regions corresponding to pathological tissues appear as distinct areas of the U-Matrix, facilitating visual identification of suspicious areas by radiologists.

4. Industrial Quality Control

In manufacturing, SOMs monitor hundreds of sensors simultaneously. Normal operation occupies a specific region of the map. The appearance of samples in new areas signals a drift or an incipient defect. This approach detects multivariate anomalies that univariate thresholds on each individual sensor would completely miss.

See Also

• Is It True That I Can Buy an iPhone with Only 1 Euro?
• Mastering Grid Point Calculation in Lattice Cubes with Python: Complete Guide and Tutorials