Laplacian Eigenmaps: Complete Guide — Principles, Examples, and Python Implementation

Summary

Laplacian Eigenmaps are one of the most elegant nonlinear dimensionality reduction methods. Introduced by Belkin and Niyogi in 2001, they rely on a fundamental principle: preserving the local structure of the data by projecting them into a lower-dimensional space while keeping neighboring points close to each other.

Unlike PCA, which seeks to maximize global variance, or Isomap, which preserves global geodesic distances, Laplacian Eigenmaps focus exclusively on local proximity preservation. This approach relies on spectral graph theory and Riemannian geometry, making it a powerful tool for exploring high-dimensional manifolds.

In this article, we will explore the mathematical principle of Laplacian Eigenmaps in depth, the underlying geometric intuition, and provide a complete Python implementation with scikit-learn.

Mathematical principle

The operation of Laplacian Eigenmaps can be broken down into four fundamental steps.

Step 1: Neighborhood graph construction

We have a set of n data points X = {x₁, x₂, …, xₙ} in a ℝᴰ space. The first step is to build a weighted graph G = (V, E) where each vertex vᵢ corresponds to a point xᵢ.

To establish edges, we typically use the k-nearest neighbors (k-NN) method: each point xᵢ is connected to its k nearest neighbors according to a Euclidean metric. We denote N(i) as the set of indices of the k nearest neighbors of xᵢ.

Step 2: Edge weighting

Once connections are defined, we assign each edge a weight wᵢⱼ that measures the similarity between points xᵢ and xⱼ. The most common weighting scheme uses a Gaussian kernel:

wᵢⱼ = exp(-‖xᵢ - xⱼ‖² / t)   if j ∈ N(i) or i ∈ N(j)
wᵢⱼ = 0                       otherwise

where t is a scale parameter (temperature) that controls the decay of the weight as a function of distance. The larger t is, the wider the kernel and the more distant points still have a noticeable influence. The smaller t is, the more only very close neighbors have significant weight.

A simple alternative is to use binary weights: wᵢⱼ = 1 if the points are connected, 0 otherwise (“nearest_neighbors” affinity in scikit-learn).

Step 3: Laplacian matrix construction

We then define the degree matrix D, which is a diagonal matrix where each diagonal element is:

Dᵢᵢ = Σⱼ wᵢⱼ

The Laplacian matrix (or graph Laplacian) is then defined as:

L = D - W

where W is the weight matrix (weighted adjacency matrix). This matrix L has remarkable mathematical properties: it is symmetric, positive semi-definite, and its eigenvalues are real and positive or zero.

There is also the normalized version of the Laplacian, denoted L_sym = I – D⁻¹⁄² W D⁻¹⁄², which has better numerical properties and is the one used by default in most practical implementations.

Step 4: Minimization and eigenvalue problem

The fundamental goal of Laplacian Eigenmaps is to find a low-dimensional representation Y = {y₁, y₂, …, yₙ} in ℝᵈ (where d ≪ D) that preserves neighborhood relationships as much as possible.

The function to be minimized is as follows:

min Σᵢⱼ ‖yᵢ - yⱼ‖² × wᵢⱼ

This expression can be elegantly rewritten in matrix form:

min yᵀ L y

under the normalization constraints yᵀ D y = 1 and yᵀ D 1 = 0 (to avoid the trivial solution).

Solving this optimization problem leads to a generalized eigenvalue problem:

L y = λ D y

The d smallest nonzero eigenvalues and their associated eigenvectors directly give the representation in dimension d. The eigenvector associated with the smallest eigenvalue (λ = 0) is ignored because it is constant, while the next d eigenvectors form the coordinates of the points in the projected space.

Geometric intuition: the spring model

To truly understand the idea behind Laplacian Eigenmaps, imagine the following situation.

You have a large number of cities on a very detailed, three-dimensional map, with mountains and valleys. Your goal is to relocate all these cities onto a flat surface (two dimensions) in such a way that cities that are geographically close remain close on the flat map.

Now, visualize each neighborhood graph link as a mechanical spring:

• Cities that are very close are connected by short springs that exert strong pull.
• Cities that are less close but still neighbors are connected by longer, more flexible springs.
• Distant cities have no spring connecting them.

When you release this spring system, it stabilizes in a configuration where:

• cities connected by strong springs (very close) do indeed remain close to each other;
• cities without a direct connection can end up anywhere — there is no constraint on their relative position.

This is exactly what Laplacian Eigenmaps do. Minimizing yᵀLy amounts to finding the equilibrium configuration of this spring system, where each weighted edge pulls the connected points together proportionally to their weight.

Where Laplacian Eigenmaps differ from other approaches is that they impose no constraints on distant points. Unlike Isomap, which seeks to preserve global distances, or PCA, which maximizes total variance, this method focuses exclusively on what happens locally. This reveals the intrinsic structure of the underlying manifold, even when it is strongly curved or coiled in the high-dimensional space.

Python implementation

Laplacian Eigenmaps are available in scikit-learn via the SpectralEmbedding class, which implements exactly the algorithm described above with a normalized Laplacian.

Basic example with Swiss Roll

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_swiss_roll
from sklearn.manifold import SpectralEmbedding

# Generate Swiss Roll data
n_samples = 1500
X, color = make_swiss_roll(n_samples, noise=0.1, random_state=42)

# Dimensionality reduction with Laplacian Eigenmaps
embedding = SpectralEmbedding(
    n_components=2,
    affinity="rbf",
    n_neighbors=15,
    gamma=0.5,
    random_state=42
)
X_embedded = embedding.fit_transform(X)

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Original data (3D to 2D projection for visualization)
axes[0].scatter(X[:, 0], X[:, 2], c=color, cmap="viridis", s=10)
axes[0].set_title("Original data (Swiss Roll)")
axes[0].set_xlabel("x₁")
axes[0].set_ylabel("x₃")

# After Laplacian Eigenmaps
axes[1].scatter(X_embedded[:, 0], X_embedded[:, 1], c=color, cmap="viridis", s=10)
axes[1].set_title("After Laplacian Eigenmaps (2D)")
axes[1].set_xlabel("Component 1")
axes[1].set_ylabel("Component 2")

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

Impact of the number of neighbors (n_neighbors)

The choice of k (the number of nearest neighbors) is crucial:

# Comparison of different n_neighbors values
k_values = [5, 15, 50]

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

for i, k in enumerate(k_values):
    embed = SpectralEmbedding(
        n_components=2,
        affinity="rbf",
        n_neighbors=k,
        gamma=1.0,
        random_state=42
    )
    X_emb = embed.fit_transform(X)

    axes[i].scatter(X_emb[:, 0], X_emb[:, 1], c=color, cmap="viridis", s=10)
    axes[i].set_title(f"n_neighbors = {k}")
    axes[i].set_xlabel("Component 1")
    axes[i].set_ylabel("Component 2")

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

With a k that is too small (e.g., k = 5), the graph risks being fragmented into several connected components, and each component will be projected independently, which can produce artificial clusters.

With a k that is too large (e.g., k = 50), the graph becomes too dense and loses its local character: the method then approaches PCA and loses its manifold-unfolding power.

The optimal value is generally between 10 and 30, depending on the density and complexity of the data.

Impact of affinity type

# Comparison of affinity types
affinities = ["rbf", "nearest_neighbors"]

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

for i, aff in enumerate(affinities):
    embed = SpectralEmbedding(
        n_components=2,
        affinity=aff,
        n_neighbors=15,
        gamma=1.0,
        random_state=42
    )
    X_emb = embed.fit_transform(X)

    axes[i].scatter(X_emb[:, 0], X_emb[:, 1], c=color, cmap="viridis", s=10)
    axes[i].set_title(f"Affinity: {aff}")

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

The “rbf” (Radial Basis Function) affinity uses the Gaussian kernel with a gamma parameter that replaces the parameter t in the formula wᵢⱼ = exp(-γ‖xᵢ – xⱼ‖²). The “nearest_neighbors” affinity uses binary weights (wᵢⱼ = 1 for connected neighbors), which is simpler but sometimes less nuanced than the Gaussian kernel.

Hyperparameters

Here are the key hyperparameters of SpectralEmbedding in scikit-learn:

Tuning tips

• n_components: start with 2 or 3 for visualization. For downstream tasks (clustering, classification), test values between 5 and 50.
• n_neighbors: the default value (5) is often too low. Try 10–30. An empirical rule is to choose k ≈ √n.
• gamma: the automatic value (median) works well in most cases. Decrease gamma to widen the kernel (more significant connections), increase it to narrow it.
• eigen_solver: for small datasets (< 1000 points), leave the automatic choice. For large datasets with sparse matrices, “arpack” is the most efficient.

Advantages and limitations

Advantages

1. Exceptional local structure preservation: neighboring points in the original space remain close in the reduced space, which is ideal for clustering and community detection.
2. Solid theoretical foundation: Laplacian Eigenmaps are intimately related to the Laplace-Beltrami Laplacian on Riemannian manifolds. As the number of points tends to infinity, the eigenvectors converge to the eigenfunctions of the Laplacian on the underlying manifold.
3. Nonlinear structure extraction: unlike PCA, which only captures linear relationships, this method effectively unfolds curved manifolds.
4. Good clustering performance: the matrix obtained after embedding is particularly well-structured for clustering algorithms like k-means (this is actually the foundation of spectral clustering).
5. No “short-circuit” problem: unlike Isomap, which can suffer from shortcuts in geodesic distance computation, Laplacian Eigenmaps are robust to these artifacts because they do not compute shortest paths.

Limitations

1. Sensitive hyperparameter choice: the number of neighbors and the kernel gamma parameter strongly influence the result. A poor choice can produce low-quality embeddings.
2. No explicit projection function: like t-SNE and UMAP, Laplacian Eigenmaps do not provide a function f(x) to project new points. The graph must be reconstructed or Nyström extension techniques must be used.
3. Computational complexity: graph construction is O(n² × D) without an accelerating structure, and diagonalization of the Laplacian is O(n³) in the worst case (though iterative solvers reduce this to O(kn²) for k eigenvectors).
4. Noise sensitivity: outlier points can create misleading connections that distort the embedding. Preprocessing (filtering, denoising) is recommended.
5. Dependence on graph connectivity: if the graph is not connected (which can happen with a k that is too small), each component is treated independently, which can produce a fragmented embedding.

4 concrete use cases

Use case 1: Visualization of genomic data

In bioinformatics, gene expression data involves thousands of genes (dimensions) measured on hundreds of samples. Laplacian Eigenmaps allow these data to be visualized in 2D while preserving natural biological groupings. Cells of the same type or samples from the same experimental conditions naturally form clusters in the reduced space.

# Conceptual example with simulated genomic data
from sklearn.manifold import SpectralEmbedding
import numpy as np

# Simulation: 200 samples, 5000 genes
np.random.seed(42)
X_genes = np.vstack([
    np.random.randn(50, 5000) + 1,   # Cell type A
    np.random.randn(50, 5000) - 1,   # Cell type B
    np.random.randn(50, 5000),        # Cell type C
    np.random.randn(50, 5000) + 0.5  # Cell type D
])

embed = SpectralEmbedding(n_components=2, affinity="rbf",
                          gamma=0.01, n_neighbors=10, random_state=42)
X_gene_emb = embed.fit_transform(X_genes)

Use case 2: Community detection in social networks

Laplacian Eigenmaps are at the heart of spectral clustering, used to detect communities in social networks. By projecting the graph nodes into a low-dimensional space (via Laplacian eigenvectors), we can then apply an algorithm like k-means to identify groups of strongly interconnected nodes.

Use case 3: Image processing and computer vision

In computer vision, Laplacian Eigenmaps are used for dimensionality reduction of images before classification. For example, in face recognition, each image is a very high-dimensional vector (pixels). Spectral embedding reveals the nonlinear structure of the face space (variations in lighting, pose, expression) much better than classical PCA.

Use case 4: Document analysis and text mining

In the field of natural language processing, documents are represented as vectors in a very high-dimensional space (bag of words, TF-IDF, embeddings). Laplacian Eigenmaps make it possible to preserve local semantic relationships: documents dealing with similar topics remain close after projection, which facilitates topic discovery and automatic categorization.

See also

• Optimize Your Python Code: Find the Maximum/Minimum Sum Sub-segment
• How to Calculate the Remainder of Polynomial Division in Python: Practical Guide for Developers