Multidimensional Scaling

Imagine you're an archaeologist who has discovered a collection of ancient artifacts at various excavation sites. You don't have GPS coordinates—only a table recording the physical distances between every pair of sites, painstakingly measured by foot. Your goal: reconstruct the spatial layout of these sites on a map. Given only pairwise distances, can you recover the original positions?

This is precisely the problem that Multidimensional Scaling (MDS) solves. It reconstructs a configuration of points in low-dimensional space such that the distances between points match, as closely as possible, a given set of observed dissimilarities. MDS is one of the oldest and most elegant dimensionality reduction techniques, predating PCA's modern formulation and serving as the conceptual ancestor of manifold learning methods like Isomap.

Classical MDS, also known as Principal Coordinate Analysis (PCoA) or Torgerson Scaling, represents the purest form of this idea: when distances are Euclidean, the reconstruction is exact and can be computed via a simple eigenvalue decomposition.

Multidimensional Scaling emerged from the intersection of psychophysics and statistics in the mid-20th century. Psychologists sought to understand how humans perceive similarity—between colors, sounds, concepts, or products. Unlike physical measurements, psychological similarity doesn't yield direct coordinates. Researchers could only ask subjects to rate how similar or different pairs of stimuli were, producing dissimilarity matrices rather than feature vectors.

The fundamental question: Given only pairwise dissimilarities, can we infer an underlying geometric structure?

The breakthrough came from Warren Torgerson in 1952, who showed that when dissimilarities correspond to Euclidean distances, the original configuration can be exactly recovered (up to rotation, reflection, and translation). This result, now called Classical MDS, remains the theoretical foundation for all subsequent MDS variants.

Key insight: Classical MDS doesn't require raw features—only distances. This makes it applicable in domains where:

• Features are unknown or unmeasurable (psychological perception)
• Distance computation is the natural primitive (string edit distances, molecular similarities)
• The original space is extremely high-dimensional (genomic distances)

Let's formalize the Classical MDS problem with mathematical precision.

Given:

• A symmetric dissimilarity matrix $D \in \mathbb{R}^{n \times n}$ where $d_{ij}$ represents the dissimilarity between objects $i$ and $j$
• $d_{ii} = 0$ for all $i$ (zero self-dissimilarity)
• $d_{ij} = d_{ji}$ (symmetry)
• A target dimensionality $k < n$

Goal:

• Find a configuration of $n$ points $X = [\mathbf{x}_1, \ldots, \mathbf{x}_n]^T \in \mathbb{R}^{n \times k}$ such that the Euclidean distances in $X$ match the given dissimilarities as closely as possible:

$$|\mathbf{x}_i - \mathbf{x}_j|2 \approx d{ij}$$

The special case: When the dissimilarities are actually Euclidean distances—meaning there exists some configuration $Z$ in $\mathbb{R}^p$ (possibly high-dimensional) such that $d_{ij} = |\mathbf{z}_i - \mathbf{z}_j|_2$—then the matrix $D$ is called Euclidean. For Euclidean distance matrices, Classical MDS can recover $Z$ exactly (up to the intrinsic dimensionality).

Key mathematical observation:

The relationship between distances and inner products is the key to Classical MDS. If we know the pairwise distances, can we recover the inner products (the Gram matrix)? The answer is yes, through a technique called double centering.

Let $X \in \mathbb{R}^{n \times p}$ be a data matrix (rows are observations). The squared distance between points $i$ and $j$ is:

$$d_{ij}^2 = |\mathbf{x}_i - \mathbf{x}_j|^2 = |\mathbf{x}_i|^2 - 2\mathbf{x}_i^T\mathbf{x}_j + |\mathbf{x}_j|^2$$

This equation links squared distances to inner products. If we had a way to isolate the inner product term $\mathbf{x}_i^T\mathbf{x}_j$, we could recover the Gram matrix $B = XX^T$, and from $B$, we can extract coordinates via eigendecomposition.

The double centering transformation is the mathematical heart of Classical MDS. It converts a squared distance matrix into a matrix of inner products, enabling eigendecomposition.

Setup: Let $D^{(2)} \in \mathbb{R}^{n \times n}$ be the squared distance matrix with entries $d_{ij}^2$.

The centering matrix: Define the centering matrix:

$$H = I_n - \frac{1}{n}\mathbf{1}\mathbf{1}^T$$

where $I_n$ is the $n \times n$ identity matrix and $\mathbf{1}$ is a vector of ones. The matrix $H$ is symmetric and idempotent ($H^2 = H$). When applied to a data matrix, it centers the columns to have zero mean.

The double centering transformation:

$$B = -\frac{1}{2} H D^{(2)} H$$

Theorem (Torgerson, 1952): If $D^{(2)}$ is a squared Euclidean distance matrix derived from centered data $X$ (with $X^T\mathbf{1} = \mathbf{0}$), then:

$$B = XX^T$$

That is, double centering recovers the Gram matrix—the matrix of all pairwise inner products.

Derivation of Double Centering:

Let's prove why double centering works. Start with the squared distance:

$$d_{ij}^2 = |\mathbf{x}_i|^2 - 2\mathbf{x}_i^T\mathbf{x}_j + |\mathbf{x}_j|^2$$

Define $a_i = |\mathbf{x}_i|^2$. Then:

$$d_{ij}^2 = a_i + a_j - 2b_{ij}$$

where $b_{ij} = \mathbf{x}_i^T\mathbf{x}_j$ is the $(i,j)$ entry of the Gram matrix.

Row mean of squared distances: $$\bar{d}{i\cdot}^2 = \frac{1}{n}\sum_j d{ij}^2 = a_i + \bar{a} - \frac{2}{n}\sum_j b_{ij}$$

Since the data is centered, $\sum_j \mathbf{x}j = \mathbf{0}$, so $\sum_j b{ij} = \mathbf{x}_i^T \sum_j \mathbf{x}j = 0$. Thus: $$\bar{d}{i\cdot}^2 = a_i + \bar{a}$$

Column mean (by symmetry): $\bar{d}_{\cdot j}^2 = a_j + \bar{a}$

Grand mean: $\bar{d}_{\cdot\cdot}^2 = 2\bar{a}$

Apply double centering: $$-\frac{1}{2}(d_{ij}^2 - \bar{d}{i\cdot}^2 - \bar{d}{\cdot j}^2 + \bar{d}{\cdot\cdot}^2)$$ $$= -\frac{1}{2}[(a_i + a_j - 2b{ij}) - (a_i + \bar{a}) - (a_j + \bar{a}) + 2\bar{a}]$$ $$= -\frac{1}{2}[-2b_{ij}] = b_{ij}$$

This proves that double centering recovers the Gram matrix entry $b_{ij}$.

Once we have the Gram matrix $B = XX^T$, we can extract coordinates through eigendecomposition. This step reveals the intimate connection between Classical MDS and PCA.

The Gram matrix is symmetric and positive semi-definite (when $D$ is Euclidean), so it has a real eigendecomposition:

$$B = V \Lambda V^T$$

where:

• $V \in \mathbb{R}^{n \times n}$ is orthogonal, with columns being eigenvectors
• $\Lambda = \text{diag}(\lambda_1, \ldots, \lambda_n)$ with $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_n \geq 0$

Coordinate recovery: Since $B = XX^T$, we can write:

$$X = V \Lambda^{1/2}$$

where $\Lambda^{1/2} = \text{diag}(\sqrt{\lambda_1}, \ldots, \sqrt{\lambda_n})$.

Dimensionality reduction: To embed in $k$ dimensions, we take only the top $k$ eigenvector-eigenvalue pairs:

$$X_k = V_k \Lambda_k^{1/2}$$

where $V_k \in \mathbb{R}^{n \times k}$ contains the first $k$ eigenvectors and $\Lambda_k$ is the $k \times k$ diagonal matrix of corresponding eigenvalues.

Choosing the embedding dimensionality:

The eigenvalues tell us how much 'variance' each dimension captures. To choose $k$, we examine:

1. Scree plot: Plot eigenvalues in decreasing order. Look for an 'elbow' where eigenvalues drop sharply.

2. Proportion of variance explained: $$\text{PoV}(k) = \frac{\sum_{i=1}^k \lambda_i}{\sum_{i=1}^n \lambda_i}$$

A common heuristic is to choose $k$ such that PoV$(k) \geq 0.8$ or $0.9$.

3. Intrinsic dimensionality: If the data truly lies on a $k$-dimensional subspace, only $k$ eigenvalues will be non-zero (up to numerical precision).

Handling negative eigenvalues:

If the input dissimilarities are not exactly Euclidean, $B$ may have negative eigenvalues. Common strategies:

1. Ignore negative eigenvalues: Only use positive eigenvalues for embedding
2. Add a constant to off-diagonal distances: This can make the matrix Euclidean
3. Use absolute values: Replace $\sqrt{\lambda_i}$ with $\sqrt{|\lambda_i|}$ for all eigenvalues
4. Switch to metric MDS: Use optimization-based approaches that handle non-Euclidean inputs

One of the most beautiful results in multivariate statistics is the equivalence between Classical MDS and PCA when applied to Euclidean distances from centered data. Understanding this equivalence deepens our appreciation for both methods and reveals why they appear in different contexts despite solving related problems.

Setup:

• Let $X \in \mathbb{R}^{n \times p}$ be a centered data matrix ($X^T\mathbf{1} = \mathbf{0}$)
• Let $D$ be the Euclidean distance matrix: $d_{ij} = |\mathbf{x}_i - \mathbf{x}_j|_2$

PCA perspective:

• Compute covariance matrix: $C = \frac{1}{n}X^TX \in \mathbb{R}^{p \times p}$
• Eigendecomposition: $C = W\Lambda_C W^T$
• Project data: $Z_{\text{PCA}} = XW$

Classical MDS perspective:

• Compute Gram matrix via double centering: $B = XX^T \in \mathbb{R}^{n \times n}$
• Eigendecomposition: $B = V\Lambda_B V^T$
• Extract coordinates: $Z_{\text{MDS}} = V\Lambda_B^{1/2}$

The equivalence theorem:

$$Z_{\text{MDS}} = Z_{\text{PCA}} \quad \text{(up to orthogonal transformation)}$$

Why does this equivalence matter?

1. Computational choice:

   • When $n < p$ (few samples, many features): use Gram matrix approach ($n \times n$ eigendecomposition)
   • When $p < n$ (many samples, few features): use covariance matrix approach ($p \times p$ eigendecomposition)
   • This is the basis for 'kernel PCA' in disguise

2. Conceptual unification:

   • PCA view: Find directions of maximum variance
   • MDS view: Find configuration that best preserves pairwise distances
   • These are the same problem with different emphases

3. Extension to kernels:

   • MDS requires only pairwise distances
   • This enables the 'kernel trick': replace $X^TX$ with any positive semi-definite kernel matrix
   • Kernel PCA is essentially MDS applied to kernel-space distances

When do they differ?

• If data is not centered, MDS on raw Euclidean distances differs from PCA
• MDS can accept any dissimilarity matrix, not necessarily Euclidean
• PCA requires the feature matrix; MDS only needs distances

After computing a $k$-dimensional embedding, we need to assess how well it preserves the original distance structure. Several metrics quantify this.

1. Proportion of Variance Explained (PoV)

The ratio of eigenvalue sum determines how much geometric information is captured:

$$\text{PoV}(k) = \frac{\sum_{i=1}^k \lambda_i}{\sum_{i=1}^n |\lambda_i|}$$

For Euclidean inputs, all eigenvalues are non-negative, and this equals the fraction of total variance. Higher PoV indicates better preservation.

2. Residual Sum of Squares (RSS)

Compare original distances to embedded distances:

$$\text{RSS} = \sum_{i < j} (d_{ij} - \hat{d}_{ij})^2$$

where $\hat{d}_{ij} = |\mathbf{x}_i - \mathbf{x}_j|$ in the embedding. Lower is better.

3. Stress (in the MDS sense)

$$\text{Stress} = \sqrt{\frac{\sum_{i<j}(d_{ij} - \hat{d}{ij})^2}{\sum{i<j} d_{ij}^2}}$$

Normalized measure of distortion. Common interpretation:

• Stress < 0.05: Excellent
• Stress < 0.10: Good
• Stress < 0.20: Fair
• Stress > 0.20: Poor

4. Shepard Diagram

A Shepard diagram plots original distances $d_{ij}$ against embedded distances $\hat{d}_{ij}$. For a perfect embedding, all points lie on the diagonal. Deviations reveal systematic distortions:

• Points below diagonal: Embedded distances too small (compression)
• Points above diagonal: Embedded distances too large (expansion)
• Curved patterns: Nonlinear distortions

5. Procrustes Analysis

If you have ground truth coordinates (e.g., in simulation), Procrustes analysis finds the optimal rotation/reflection/scaling to align the embedding with truth, then measures residual error:

$$\text{Procrustes error} = \min_{R \in \mathcal{O}(k)} |X_{\text{true}} - X_{\text{MDS}} R|_F$$

where $\mathcal{O}(k)$ is the set of orthogonal matrices.

Classical MDS is computationally elegant but faces scalability challenges for large datasets. Let's analyze the computational requirements.

Time Complexity:

1. Distance matrix computation (if starting from features): $O(n^2 p)$
2. Squaring distances: $O(n^2)$
3. Double centering: $O(n^2)$
4. Eigendecomposition of $n \times n$ matrix: $O(n^3)$
5. Coordinate extraction: $O(nk)$

Total: $O(n^2 p + n^3)$

For $n = 10,000$ points, the eigendecomposition alone requires $\sim 10^{12}$ operations—potentially hours of computation.

Space Complexity:

• Distance matrix $D$: $O(n^2)$
• Gram matrix $B$: $O(n^2)$
• Eigenvectors $V$: $O(n^2)$

Total: $O(n^2)$

For $n = 50,000$ with float64: $50000^2 \times 8 \approx 20$ GB per matrix. This quickly exceeds RAM.

Optimization strategies:

1. Exploit sparsity in distance computation

• If only $k$-nearest neighbor distances are needed, use $k$-d trees or ball trees
• Approximate nearest neighbors (e.g., HNSW, Annoy) for very large $n$

2. Partial eigendecomposition

• If only $k$ components are needed, use iterative methods (Lanczos, randomized SVD)
• scipy.sparse.linalg.eigsh or sklearn.utils.extmath.randomized_svd
• Reduces complexity from $O(n^3)$ to $O(n^2 k)$ or better

3. Landmark MDS (L-MDS)

• Select $m \ll n$ landmark points
• Compute exact MDS on landmarks: $O(m^3)$
• Project remaining points using distances to landmarks: $O((n-m) \cdot m)$
• Total: $O(m^3 + nm)$

4. Nyström approximation

• Sample $m$ columns/rows of the Gram matrix
• Approximate full eigendecomposition from this sample
• Well-studied approximation guarantees

5. Out-of-sample extension

• Fit MDS on a representative subset
• Project new points without refitting

Classical Multidimensional Scaling provides a foundational approach to understanding and visualizing data from pairwise distances. Let's consolidate the key concepts:

What's next:

Classical MDS assumes that dissimilarities are Euclidean distances and optimizes a specific (implicit) objective. In the next page, we explore Metric MDS, which explicitly optimizes a stress function, allowing us to handle a broader class of dissimilarities and to weight distance preservation according to importance. Metric MDS introduces iterative optimization, providing a more flexible framework at the cost of computational efficiency.