= Solution
<Classical multidimensional scaling> begins with a symmetric <dissimilarity matrix> $D$, rather than necessarily with measured coordinates. Its aim is to find a low-dimensional Euclidean configuration representing those dissimilarities. Square the entries, not the matrix product, and define
$$
J=I-\frac1n\mathbf1\mathbf1^T,\qquad B=-\frac12J D^{(2)}J.
$$
For centered coordinates $X$, squared <Euclidean distances> satisfy $D_{ij}^2=B_{ii}+B_{jj}-2B_{ij}$ with $B=XX^T$. Multiplication on both sides by $J$ annihilates the two single-index terms, proving the double-centering formula.
If $B$ is <positive semidefinite>, use its <spectral decomposition> $B=U\Lambda U^T$ to construct
$$
\boxed{X_q=U_q\Lambda_q^{1/2}.}
$$
Keeping all positive <eigenvalues> reproduces the original <Euclidean distances> exactly; their number is the minimal embedding dimension. Keeping fewer gives a low-rank approximation to the <Gram matrix>. It minimizes squared Gram-matrix error, often called strain, rather than necessarily the sum of squared errors in raw distances. Negative <eigenvalues> signal that the dissimilarities cannot be represented exactly by Euclidean coordinates; discarding them produces an approximation whose adequacy should be inspected.
In the upper-right sketch, the supplied distances come from four points in a rectangle, and the two-dimensional reconstruction recovers that shape. The overall translation, rotation and reflection of the configuration are not identifiable from distances. \b[Classical scaling reconstructs centered geometry from pairwise distances.] If the input distances are computed from centered multivariate observations, classical scaling and <principal component analysis> give the same score configuration up to an orthogonal transformation: their matrices $ZZ^T$ and $Z^TZ$ share the nonzero squared singular values.
Back to article page