Higher-dimensional metric MDS for 3D visualization

Three routes to a 3D drawing

1. Direct 3D: minimize metric-MDS stress with three coordinates per vertex.
2. Higher-dimensional MDS + PCA: fit, for example, ten coordinates per vertex; center them and retain the first three principal components.
3. Higher-dimensional MDS + PCA + 3D refinement: use those three PCA coordinates as the starting configuration for another metric-MDS fit. Every vertex can then move freely in 3D.

PCA chooses a three-dimensional linear subspace that retains coordinate variance. It does not minimize error against the original distance targets. Free 3D refinement optimizes those original targets again; the result need not remain a linear projection of the higher-dimensional coordinates. This is different from the optional edge.kk() stage in the synthetic graph gallery, which fits directly connected vertices rather than all finite pairs.

An optimized higher-dimensional fit can have low stress while its PCA projection has much higher stress. Evaluate the actual three-column result, not just the higher-dimensional fit. Also compare a classical-MDS start followed by 3D metric MDS: a useful initialization need not require a higher-dimensional iterative fit.

A runnable 10D example

We use the bundled 200-point paraboloid with Euclidean edge lengths on a symmetric six-nearest-neighbor graph. Its targets are shortest-path distances. Preparing the graph once lets every route use exactly the same target matrix.

library(grip)
bundle <- readRDS(system.file('extdata','high-dimensional-mds','examples.rds',package='grip'))
example <- bundle$cases[['paraboloid_17_6']]
prepared <- prepare.graph.geodesic.mds(
  edges=example$edges, n=example$n, edge.weights=example$weights)
D <- prepared$distance_matrix
stopifnot(isTRUE(all.equal(unname(D),unname(example$distances),tolerance=1e-10)))

The following is a short usage demonstration, with 20 SGD passes per fitted stage and one start. It is not a convergence study or an equal-time benchmark. The 10D route performs an additional fitting stage. To expose the initialization, we supply random matrices explicitly; the direct start uses the first three columns of the same draw. Internal scale adjustment differs between dimensions.

set.seed(17)
start10 <- matrix(rnorm(example$n * 10),example$n,10)
fit <- function(dimension, initial) {
  # Record schedule termination below; suppress its expected warning here.
  withCallingHandlers(
    metric.mds(prepared=prepared, dim=dimension, init=initial,
      backend='sgd', approximation='full', pair.weights='uniform',
      n.init=1, max.iter=20, seed=17, diagnostics=FALSE,
      sgd.control=list(checkpoint.every=1L)),
    warning=function(w) {
      if (grepl('iteration_limit',conditionMessage(w),fixed=TRUE))
        invokeRestart('muffleWarning')
    })
}
direct <- fit(3,start10[,1:3])
high <- fit(10,start10)
pca3 <- prcomp(high$coords,center=TRUE,scale.=FALSE,rank.=3)$x[,1:3,drop=FALSE]
refined <- fit(3,pca3)
classical.start <- fit(3,'classical')

Use scale.=FALSE: rescaling the ten coordinate axes to unit variance would change their Euclidean geometry before projection. project.3d() is a drawing helper, not PCA; it uses the first three coordinate columns.

We score the returned coordinates without an additional fitted scale:

\[ S_{\mathrm{norm}}(Z)= \frac{\sum_{i<j:\ d_{ij}<\infty}(\lVert z_i-z_j\rVert-d_{ij})^2} {\sum_{i<j:\ d_{ij}<\infty}d_{ij}^{2}}. \]

Lower is better; zero is an exact realization. This dimensionless normalized stress is the square of the normalized error reported in the backend comparison. The metric-MDS solver itself fits an overall scale when returning coordinates. Here we evaluate those coordinates as returned, and the PCA coordinates as projected, without another rescaling during evaluation.

normalized.stress <- function(z,D) {
  target <- as.vector(as.dist(D))
  keep <- is.finite(target)
  sum((as.vector(dist(z))[keep]-target[keep])^2)/sum(target[keep]^2)
}
layouts <- list('Direct 3D'=direct$coords, '10D before projection'=high$coords,
  '10D + PCA'=pca3, '10D + PCA + 3D refinement'=refined$coords,
  'Classical + 3D refinement'=classical.start$coords)
knitr::kable(data.frame(Route=names(layouts),
  Dimensions=vapply(layouts,ncol,0L),
  `Normalized stress`=vapply(layouts,normalized.stress,0,D=D),
  check.names=FALSE),digits=7,row.names=FALSE)

iteration_limit means that the prescribed schedule ended, not that an optimum was reached. These fresh demonstration fits are separate from the saved results below; their starts and implementation wrapper are not identical.

Saved examples and comparison conditions

The saved study used full, uniformly weighted SGD-MDS on nine 200-point clouds (plane, paraboloid, saddle; three data seeds), each with five distance targets, plus six SuiteSparse graphs. It retained dimensions 4–20 and optimization starts 17, 29, and 43. Here we show four target datasets and dimensions 4, 6, 8, 10, 15, and 20. The examples were selected after examining the study to teach contrasting outcomes; they are not an unbiased success-rate sample. Every displayed fit uses start 17, without choosing its appearance or score. Tables retain all three starts; classical MDS is deterministic and appears once.

Synthetic fits allowed 20 passes per stage; SuiteSparse fits allowed four. The PCA-refined route therefore includes both the higher-dimensional stage and a second 3D stage. Route seconds include the preceding fitting stages, their initialization and projection, but exclude shared target preparation and subsequent evaluation/rendering. Equal pass counts across dimensions do not mean equal time or arithmetic work.

To compare cost, the study also continued direct 3D SGD and retained its best scored checkpoint within each refined route’s charged time allowance. The dimension tables compare against that stronger direct baseline. The allowance includes some archival overhead absent from direct continuation, favoring the direct method. Times describe these short recorded runs on one machine, not isolated solver speed or performance on an idle machine. Neither route is claimed to have converged.

The viewers show the saved initial direct stage, whose coordinates exist. The time-matched direct checkpoints have saved scores but not saved coordinates; their scores appear only in tables. A displayed direct drawing must therefore not be mistaken for the time-matched comparator.

All views preserve coordinate scale. Display alignment permits translation, rotation and reflection using corresponding vertex identities. It does not deform the drawing or affect scores. Gray reference points, where present, show the generating cloud; similarity to that cloud is a different question from fitting graph-path distances. Blue denotes direct/classical fits, orange PCA projections, and green refined fits. Use individual selections before overlaying all layouts.

One paraboloid, two distance targets

The 200 points have data seed 17 and satisfy \(z=x^2+y^2\), with parameters sampled uniformly from the unit disk. This samples the parameter domain rather than surface area uniformly.

SuiteSparse: an improvement and a counterexample

These are unweighted graph conversions of HB/illc1033 and HB/lock_700 from the SuiteSparse Matrix Collection, not embeddings of their original numeric matrix values. HB/illc1033 uses distinct row and column vertices with an edge for each nonzero matrix entry. HB/lock_700 uses symmetrized off-diagonal support. Every edge has unit traversal length. The package’s suitesparse.graph() catalogue retains the original source attribution and conversion metadata.

Disconnected components were fitted separately, and pairs in different components are excluded from stress. All vertices are retained. HB/lock_700’s nine isolates are placed in a separate display row; their positions relative to the main component have no distance-based interpretation. For legibility, each dense graph view draws at most 2,500 edges, selected at evenly spaced indices in the saved edge order. Fitting and scoring use the full graph.

Choosing a workflow

Start with direct 3D MDS and a classical-initialized 3D control. Try several random starts before attributing a difference to dimension. If higher-dimensional fitting is promising, compare PCA-only and refined results at a few dimensions and include the total route cost. Ten is an example, not a recommended optimum; the original study’s preselected summary dimension was eight.

Across all six study graphs, refinement at dimension eight lowered mean time-matched stress by more than 1% on three graphs. The short schedules and small number of starts limit that finding. Choosing the best dimension after examining these tables is an exploratory choice, not independent validation. Stress also does not determine visual usefulness: inspect neighborhoods, overlap and recognizable structure in the interactive views.

Full MDS requires dense all-pairs distances and quadratic work/storage. Larger problems may need sparse MDS. Sparse fitting has a different pair set and objective; its proxy score cannot be substituted for the full uniform stress used here. Higher-dimensional sparse fitting is supported, but its performance is not tested by these saved comparisons.

Reproduction and saved evidence

inst/extdata/high-dimensional-mds/examples.rds contains the four graph fixtures, two small synthetic distance matrices, the displayed 3D coordinates, and selected three-start scores/timings. Its source-hash table records the original study files. The saved study used a timing-instrumented SGD wrapper with component handling; the short public-API recipe above does not reproduce its exact starts, checkpoints or timings.

From a source checkout, the maintainer can extract this subset from the retained study archive, verifying its hashes and independently rescoring every included 3D layout, with:

Rscript data-raw/build-high-dimensional-mds.R /path/to/high_dimensional_mds_to_3d
make high-dimensional-vignette

The extraction runs no optimization. Rebuilding the vignette runs only the small public-API demonstration; the saved tables and 3D views require no external archive, server, or gflowui installation. The bundle’s provenance document describes selection, timing and attribution. See also the function guide for dimension support and drawing conventions.