2.2 Classical and non-metric multidimensional scaling

Some datasets describe relations between objects rather than measurements of their features. We may know how far apart every pair of objects is without knowing a useful coordinate system for them. Multidimensional scaling asks how to turn those pairwise relations into positions in a low-dimensional Euclidean space. Classical MDS reconstructs inner products from squared distances; non-metric MDS instead fits the order of dissimilarities through a monotone transformation. Their objectives express different meanings of a good map.

Motivation

PCA starts with a table of feature vectors and finds directions in their feature space. MDS starts with a square table whose entries compare pairs of objects. A row of that table lists relations to the other objects, not a vector of measured features. Applying PCA directly to those rows would therefore solve a different problem from reproducing the supplied distances.

Let $D=(d_{ij})\in\mathbb R^{n\times n}$, with $n\ge2$. The target is a coordinate matrix $X\in\mathbb R^{n\times k}$: row $i$ stores the point $x_i\in\mathbb R^k$. We want the Euclidean distances between rows to express the pairwise information in $D$. Typical maps use one or two coordinates, but the algebra also identifies how many coordinates exact reconstruction needs.

A distance table cannot determine an absolute origin. Adding the same vector $c\in\mathbb R^k$ to every row preserves all differences between points. We therefore choose the centered representative, whose column means vanish. An origin is fixed by this choice; an orientation remains free.

Definitions

Definition

Dissimilarity, metric, and Euclidean distance tables

A dissimilarity table here is symmetric, nonnegative, and hollow: $d_{ij}=d_{ji}\ge0$ and $d_{ii}=0$. A metric additionally separates distinct objects and satisfies $d_{ij}\le d_{i\ell}+d_{\ell j}$ for every triple. A table is Euclidean if some points satisfy

$$d_{ij}=\|x_i-x_j\|_2.$$

Coincident points are allowed in this last definition; they produce zero distance between their labels. The triangle inequality alone does not characterize Euclidean distance tables.

For classical MDS, square each entry separately and define the centering matrix:

$$E=D^{(2)},\qquad E_{ij}=d_{ij}^2, \qquad H=I_n-\frac1n\mathbf1\mathbf1^T.$$

Thus $E$ is not the matrix product $DD$. The matrix $H$ is symmetric, $H^2=H$, and $H\mathbf1=0$. Multiplication $HX$ subtracts the coordinate mean from each row. Define row means and the grand mean of $E$ by

$$r_i=\frac1n\sum_j E_{ij},\qquad g=\frac1{n^2}\sum_{i,j}E_{ij}.$$

These are means of squared distances, rather than means of distances. Because $E$ is symmetric, its column means equal its row means. Double centering subtracts both sets of means and adds back the grand mean once.

Theorem and proposition

Theorem

The centered Gram matrix and exact Euclidean reconstruction

For a symmetric, nonnegative, hollow table $D$, set

$$B=-\frac12HEH,\qquad b_{ij}=-\frac12(E_{ij}-r_i-r_j+g).$$

Then $D$ is Euclidean exactly when $B$ is positive semidefinite. When this holds, its positive eigendecomposition gives centered coordinates

$$B=V_+\Lambda_+V_+^T,\qquad X=V_+\Lambda_+^{1/2},\qquad XX^T=B.$$

The smallest exact embedding dimension is $\operatorname{rank}(B)$, at most $n-1$. A zero-rank table describes coincident points.

Proposition

Positive spectral truncation and classical strain

Order the eigenvalues of $B$ from largest to smallest. Among centered positive semidefinite Gram matrices $G$ of rank at most $k$, the minimum of $\|B-G\|_F^2$ is attained by retaining the largest at most $k$ positive eigenvalues. If $I_k$ indexes those retained values, then

$$G_k=\sum_{j\in I_k}\lambda_jv_jv_j^T, \qquad X_k=V_{I_k}\Lambda_{I_k}^{1/2}, \qquad \|B-G_k\|_F^2=\sum_{j\notin I_k}\lambda_j^2.$$

Negative eigenvalues contribute unavoidable error to this positive Gram approximation. Repeated eigenvalues can give several equally good coordinate orientations or retained subspaces.

The discrepancy in the second statement is often called classical strain. It compares inner products, so it is distinct from the sum of squared errors in the pairwise distances themselves. A low-dimensional classical map optimizes this Gram criterion; interpreting its coordinates still requires examining the distances and the discarded spectrum.

Deriving double centering

Suppose centered coordinates exist. Their Gram matrix is $B=XX^T$, with $b_{ij}=\langle x_i,x_j\rangle$. The distance identity gives

$$E_{ij}=b_{ii}+b_{jj}-2b_{ij}.$$

Let $a=(b_{11},\ldots,b_{nn})^T$. In matrix form this identity reads $E=a\mathbf1^T+\mathbf1a^T-2B$. Multiplying on both sides by $H$ removes the first two terms because $H\mathbf1=0$. Centering gives $HX=X$, so $HBH=B$. Hence $HEH=-2B$, proving the formula using only the distance table. The entrywise mean formula is precisely the expansion of this matrix product.

Proof

Why positive semidefiniteness is also sufficient

If coordinates exist, then $z^TBz=\|X^Tz\|_2^2\ge0$ for every $z$, so $B$ is positive semidefinite. Conversely, assume that the constructed $B$ is positive semidefinite and factor it spectrally as $XX^T$. Since $B\mathbf1=0$, every positive-eigenvalue eigenvector is perpendicular to $\mathbf1$; the resulting coordinate columns have mean zero.

The hollow diagonal of $E$ implies $b_{ii}=r_i-g/2$. Therefore

$$b_{ii}+b_{jj}-2b_{ij} =(r_i-g/2)+(r_j-g/2)+(E_{ij}-r_i-r_j+g) =E_{ij}.$$

The recovered points reproduce the squared distances, and nonnegativity of $D$ recovers the distances themselves. A Gram matrix of coordinates in $k$ dimensions has rank at most $k$, while the factor above uses exactly one column per positive eigenvalue. This proves the minimum-dimension claim.

For the truncation result, orthogonal invariance of the Frobenius norm allows the comparison in an eigenbasis of $B$. The spectral low-rank approximation principle assigns available nonnegative eigenvalues to the largest positive directions; in each selected direction the best value is $\lambda_j$. Unselected directions contribute their squared eigenvalues. Every retained positive direction is centered because it is orthogonal to $\mathbf1$.

Centered coordinates have a remaining geometric ambiguity. For an orthogonal $k\times k$ matrix $Q$, the matrix $XQ$ has the same Gram matrix and the same pairwise distances. Rotations, reflections, and changes of eigenvector sign consequently produce equally valid maps. A meaningful comparison uses relative positions and distances, rather than the sign of a plotted axis.

Worked examples

Worked example

Example 1: Recover three collinear points

Take

$$D=\begin{bmatrix}0&1&3\\1&0&2\\3&2&0\end{bmatrix}, \qquad E=\begin{bmatrix}0&1&9\\1&0&4\\9&4&0\end{bmatrix}.$$

The squared-distance row means are $10/3$, $5/3$, and $13/3$, with grand mean $28/9$. Double centering gives

$$B=\frac19\begin{bmatrix}16&4&-20\\4&1&-5\\-20&-5&25\end{bmatrix}.$$

This matrix equals $uu^T/9$ for $u=(-4,-1,5)^T$. Its only positive eigenvalue is $\|u\|_2^2/9=14/3$. Consequently one coordinate suffices:

$$X=\begin{bmatrix}-4/3\\-1/3\\5/3\end{bmatrix}.$$

The coordinate mean is zero, and the three pairwise distances are $1$, $2$, and $3$ in their prescribed positions. Adding $4/3$ to each coordinate gives $0$, $1$, and $3$; reflection gives another valid arrangement. These changes illustrate the translation and orientation information that a distance table does not contain.

The negative off-diagonal entries of $B$ are ordinary negative inner products. They do not contradict positive semidefiniteness, which is a condition on quadratic forms or eigenvalues rather than the signs of individual entries.

Worked example

Example 2: A metric that is not Euclidean

Let $D$ record shortest-path distances around a cycle of four unit edges:

$$D=\begin{bmatrix}0&1&2&1\\1&0&1&2\\2&1&0&1\\1&2&1&0\end{bmatrix}, \qquad B=\frac14\begin{bmatrix}3&1&-5&1\\1&3&1&-5\\-5&1&3&1\\1&-5&1&3\end{bmatrix}.$$

The triangle inequality holds: an opposite pair is distance $2$, reached by two adjacent edges of length $1$, and no alternate path is shorter than the shortest path. Nevertheless the eigenvalues of $B$ are $2$, $2$, $0$, and $-1$. For example, $(1,-1,1,-1)^T$ is an eigenvector with eigenvalue $-1$. The exact Euclidean criterion therefore fails in every dimension.

Retaining the two positive directions gives coordinates

$$X_2=\begin{bmatrix}1&0\\0&1\\-1&0\\0&-1\end{bmatrix}.$$

Opposite distances remain $2$, while adjacent distances become $\sqrt2$. The minimum squared Gram discrepancy is $(-1)^2=1$. The negative direction records a genuine obstruction to the requested geometry; it is not an extra Euclidean coordinate whose square root could be added. The approximation replaces that geometry by a positive Gram matrix and visibly changes some distances.

Worked example

Example 3: Exact ordinal structure without exact input distances

For three objects, prescribe dissimilarities $d_{12}=1$, $d_{23}=2$, and $d_{13}=8$. They violate the triangle inequality because $8>1+2$. Their strict order can nevertheless be represented by the one-dimensional coordinates $0$, $1$, and $3$, whose corresponding distances are $1$, $2$, and $3$.

A monotone transformation can take

$$f(1)=1,\qquad f(2)=2,\qquad f(8)=3.$$

Join these values with increasing line segments if a function between the observed dissimilarities is required. Every transformed dissimilarity equals its fitted Euclidean distance. The normalized stress defined below is zero, with positive denominator $1^2+2^2+3^2=14$.

The map captures the ranking perfectly, even though the input value $8$ is not reproduced. This is exactly the information an ordinal objective chooses to retain. The numerical gaps between successive dissimilarities are free to change, so zero ordinal stress is a different conclusion from exact classical distance reconstruction.

A numerical classical MDS calculation

The following calculation implements Example 1 with symmetric eigenvalue decomposition. The explicit elementwise square distinguishes the input distances from the Gram construction. A small positive threshold separates the zero eigenvalues in this example from floating-point roundoff.

import numpy as np

D = np.array([[0, 1, 3], [1, 0, 2], [3, 2, 0]], dtype=float)
n = D.shape[0]
H = np.eye(n) - np.ones((n, n)) / n
E = D ** 2
B = -0.5 * H @ E @ H
eigenvalues, eigenvectors = np.linalg.eigh(B)
order = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[order]
eigenvectors = eigenvectors[:, order]
selected = np.flatnonzero(eigenvalues > 1e-10)[:1]
X = eigenvectors[:, selected] * np.sqrt(eigenvalues[selected])
fitted = np.linalg.norm(X[:, None, :] - X[None, :, :], axis=2)
print(np.round(eigenvalues, 6))
print(np.round(fitted, 6))

The eigenvalues are approximately $14/3$, $0$, and $0$; the reconstructed distance matrix equals $D$ up to numerical precision. A returned coordinate column may have the opposite sign from the hand calculation. When inspecting another dataset, examine the negative spectrum before truncating it: the size of a substantive negative eigenvalue measures information discarded by the positive Gram approximation.

The connection with PCA

Suppose the feature measurements are available after all. Let $A\in\mathbb R^{n\times d}$ be their centered matrix, and form $D$ from Euclidean distances between its rows. The derivation above then yields

$$B=AA^T=U\Sigma^2U^T,\qquad X_k=U_k\Sigma_k=AV_k.$$

Here the last expression is precisely the leading PCA score matrix. Classical MDS and PCA consequently produce the same reduced configuration when they start with this same Euclidean geometry. One starts from distances and the other from features; the common object is the centered inner-product structure. Signs and rotations within repeated eigenspaces can change the displayed coordinate axes while preserving that structure.

The eigenvalues in this Gram matrix are squared singular values, rather than sample covariance eigenvalues. Sample covariance divides by $n-1$. Taking square roots of those covariance eigenvalues as MDS coordinate scales would shrink every distance by a factor of $\sqrt{n-1}$. Retaining the correct normalization matters when a map must reproduce given numerical distances.

The same comparison helps interpret dimension reduction. Discarding a positive eigenvalue removes a direction that the complete Euclidean data actually possess. Discarding a negative eigenvalue removes a direction that cannot belong to any real Euclidean Gram factor in the first place. These two sources of approximation enter the strain sum together, but represent different geometric reasons for a mismatch. Inspecting their contributions separately explains whether an inaccurate map mainly reflects the chosen dimension or the non-Euclidean structure of the input relations.

Non-metric MDS and stress

Definition

Monotone disparities and normalized stress

For fitted distances $\widehat d_{ij}=\|x_i-x_j\|_2$, choose a nondecreasing function $f$ of the observed dissimilarities. The values $f(d_{ij})$ are called disparities. The normalized stress used here is

$$S(X,f)=\left( \frac{\sum_{i<j}(f(d_{ij})-\widehat d_{ij})^2} {\sum_{i<j}\widehat d_{ij}^2} \right)^{1/2}, \qquad \sum_{i<j}\widehat d_{ij}^2>0.$$

Equal input dissimilarities receive equal disparities under this function convention. Distinct inputs may receive tied disparities because nondecreasing functions can have flat portions.

With $X$ fixed, the denominator is constant. Sort the observed pairs by their dissimilarities and fit their current Euclidean distances with a nondecreasing sequence. This is isotonic regression. If successive fitted values violate the ordering, pooling a block replaces its values by their average in the unweighted least-squares step. Tied inputs are treated as one disparity level.

With disparities fixed, change the coordinates to reduce the chosen stress. The coordinate problem is nonconvex, so initialization and the progress of this step affect the attained configuration. Alternating the two steps describes an optimization scheme; a stopping criterion terminates that scheme without establishing a global optimum. In particular, small stress measures closeness to monotone disparities, rather than certifying every strict order between fitted distances.

The positive-denominator condition excludes a completely collapsed configuration, where every point is identical and this ratio is undefined. For $a>0$, replacing $(X,f)$ by $(aX,af)$ leaves stress unchanged. A numerical scheme may therefore fix an overall scale as well as center its coordinates. These choices remove irrelevant location and scale freedoms while leaving the ordinal comparison intact.

Common mistakes

Common mistake

Applying an appropriate criterion to the wrong object

Square the entries of $D$ before double centering. Test positive semidefiniteness of the resulting Gram matrix, rather than requiring its entries to be nonnegative. A valid metric can still have an indefinite centered Gram matrix. Classical truncation, raw-distance fitting, and ordinal stress compare different objects, so their numerical errors are not interchangeable.

Summary

Classical MDS turns squared pairwise distances into a centered Gram matrix. Its spectrum decides exact Euclidean embeddability and supplies coordinates or a positive low-rank approximation. Non-metric MDS instead compares fitted distances with a monotone transformation of dissimilarities. Centering fixes translation; orthogonal changes preserve classical geometry, and a joint positive scaling of coordinates and disparities preserves normalized ordinal stress.

Exercises

Checkpoint

Exercise 1: Two objects have distance $4$. Compute their centered Gram matrix, its positive eigenvalue, and one centered one-dimensional configuration.

Use the hollow matrix with off-diagonal entries $4$ and square those entries before centering.

Checkpoint

Exercise 2: A centered Gram matrix has eigenvalues $5$, $2$, $0$, and $-1$. Can it reproduce Euclidean distances exactly? What is its minimum squared Gram discrepancy using one positive direction?

Separate the Euclidean criterion from the rank-one approximation criterion.

Checkpoint

Exercise 3: Keep input dissimilarities $d_{12}=1$, $d_{23}=2$, $d_{13}=8$, but fix coordinates at $0$, $2$, $3$. Find the isotonic disparities and the resulting normalized stress.

In increasing dissimilarity order, the current distances are $2$, $1$, $3$. Pool the violating adjacent block.

Solutions

Solution · Exercise 1 solution

The mean squared-distance row value and grand mean both equal $8$, giving

$$B=\begin{bmatrix}4&-4\\-4&4\end{bmatrix}.$$

Its eigenvalues are $8$ and $0$. Coordinates $(-2,2)^T$ have mean zero, Gram matrix $B$, and distance $4$. The positive rank is one.

Solution · Exercise 2 solution

The negative eigenvalue prevents exact Euclidean reconstruction in any dimension. The rank-one positive approximation retains $5$. Its squared Gram discrepancy is $2^2+0^2+(-1)^2=5$. The retained coordinate column is $\sqrt5$ times a unit eigenvector for the eigenvalue $5$.

Solution · Exercise 3 solution

The first two distances violate the required nondecreasing order. Their least-squares pooled value is $(2+1)/2=3/2$, while the third disparity remains $3$. Thus

$$f(1)=f(2)=\frac32,\qquad f(8)=3,\qquad S=\sqrt{\frac{(3/2-2)^2+(3/2-1)^2+(3-3)^2}{2^2+1^2+3^2}} =\frac1{\sqrt{28}}.$$

The flat fitted block is permitted by nondecreasing isotonic regression. It does not remove the distance inversion while the coordinates remain fixed; the subsequent coordinate step addresses that remaining mismatch.

Related notes

SVD and principal components introduces centering and spectral coordinates from measured features. Centered kernels and kernel PCA starts from inner products in a feature space. Fisher discriminants and class scatter uses class labels to choose projection directions. The course reading list collects the course books.