Sparse SGD for larger graphs

Use sparse SGD

g <- suitesparse.graph('dwt_307')
fit <- metric.mds(edges=g$edges,n=g$n,edge.weights=g$edge_weights,
  dim=3,approximation='sparse',backend='sgd',seed=1,
  sparse.control=list(n.pivots=32),max.iter=100)
ivue::plot3D.plain(fit$coords, layers=list(ivue::layer3D.edges(g$edges)), controls = FALSE)

Full MDS remains the default. Sparse mode currently supports SGD only, with inverse-squared stress weights. Its defaults are random initialization, up to 200 pivots, and 30 epochs: provisional calibration choices. The 32 pivots and 100 epochs above are explicit illustrative choices. Supply connected graphs with positive lengths; disconnected graphs are rejected without adding edges.

Sparse mode accepts supplied starting coordinates but does not accept classical initialization, uniform weights, multiple starts, dense prepared objects, or dense diagnostics. Observation-based sparse preparations are supported as described below. These restrictions avoid an unexpected all-pairs calculation. Its sparse_proxy_stress describes the fitted sparse constraints, not full stress. When a supplied edge is longer than an alternative route, sparse and full MDS can use different targets for that edge. That distinction does not arise for the unit and Euclidean chord lengths used in this experiment.

Prepare directly from observations

prepare.sparse.mds() constructs the pivots and weights internally, without an all-pairs distance matrix. The default method="euclidean" uses original-space Euclidean pivot distances, published distance-proportional randomized pivot selection, regional endpoint weights and symmetric nearest-neighbor local pairs. Neighbor search uses the smallest ordinary PCA projection reaching 90% variance; all edge lengths and fitting targets still use the original observations.

prep <- prepare.sparse.mds(X, method="euclidean", n.pivots=200,
  q=10, neighbor.variance.target=0.90, seed=1)
fit <- metric.mds(prepared=prep, approximation="sparse", backend="sgd",
  dim=3, max.iter=100, seed=1,
  sgd.control=list(scheduler="zheng", schedule.epsilon=0.1))

Use pivot.selection="farthest" for the Euclidean or L-infinity farthest-first alternative. method="uniform" uses uniform landmarks without regional multipliers or local edges, but retains inverse-squared distance weights. method="geodesic" instead uses shortest-path targets on a connected neighbor graph. Its exact component-MST repair evaluates distances as needed and never stores their full matrix; worst-case repair time can nevertheless be quadratic.

For maximum absolute coordinate differences, use method="linf". It computes exact neighbors in the original coordinates; equal distances are ordered by row index. Its default neighbor.variance.target is NULL, and explicit PCA targets are rejected because rotation changes this metric. No normalization is implicit. For nonnegative abundance rows with positive maxima, max normalization is:

stopifnot(all(is.finite(X)), all(X >= 0), all(apply(X, 1, max) > 0))
Xmax <- X / apply(X, 1, max)
# Consolidate any duplicate normalized rows and retain their sample mapping.
prep <- prepare.sparse.mds(Xmax, method="linf", n.pivots=200, q=10, seed=1)

Exact L-infinity search has quadratic worst-case work but does not store an all-pairs matrix. Its returned graph$edges table holds the symmetric union of nearest-neighbor pairs and their L-infinity lengths; graph$metadata records the search rule. It needs no dgraphs dependency.

The Euclidean and geodesic methods require the accompanying dgraphs support for search.data and streaming component repair. Input rows must be distinct. Requested pivot counts are capped at nrow(X)-1; keep the count small relative to the data size for useful storage savings. The preparation can be reused across dimensions and fit seeds. Its default schedule is "zheng", using the published weight-based endpoints; the legacy graph interface keeps its existing default schedule. Completing the epoch budget does not certify convergence. See ?prepare.sparse.mds for returned audit information and controls.

The October 2026 source-tree methods note under vignettes/articles/sparse-mds-methods/ documents formulas and a separate 288-fit quadratic-form validation. The existing experiment below concerns the legacy graph interface and its recorded settings.

Quality, time, and memory on existing examples

The saved experiment compares full SGD, sparse SGD at several pivot counts, and weighted GRIP. Each setting uses seeds 1–3. SGD uses 100 epochs with the same random initial coordinates for paired runs; GRIP uses its own initialization. Full MDS retains its best checkpoint and fits scale; sparse SGD returns its terminal iterate without fitting scale. These are comparisons of the supplied methods, not equal work per epoch.

The structural graph below connects this comparison to the graph gallery. The larger saddle illustrates approximation at 4,096 vertices. Full SGD was not attempted above 512 vertices in this bounded experiment; that cutoff is not a measured capacity limit.

Errors are independently evaluated against graph shortest paths: all pairs for dwt_307, and all distances from 32 independently sampled source vertices for the larger saddle. Original error uses returned coordinates. Scale-adjusted error additionally fits one positive scalar, a diagnostic added after inspecting the original run; no layouts were refitted.

Total seconds includes graph preparation and fitting. Peak MiB includes the entire isolated R process, including package loading. Tables report medians over returned runs and show both the returned and attempted counts.

Complete results and reproduction

The experiment appendix preserves every example, including the smaller saddle, paraboloid, and helix, all outcomes, numerical protocol, resource limits, and reproduction commands. A returned fit does not establish convergence. Compare independently evaluated error rather than sparse proxy values from different pivot counts.

Explicit distance constraints

Sparse SGD also accepts distances already computed by another method. Supply n and a constraints list instead of graph edges. Targets are used as given; no shortest-path calculation or automatic pivot selection is performed. This is useful for exact complete-graph Fermat distances computed at selected pairs without storing a complete graph.

constraints <- list(pairs = rbind(c(1, 2), c(2, 3), c(1, 3)),
                    targets = c(1, 1, 2))
fit <- metric.mds(n = 3, approximation = "sparse", constraints = constraints,
                  seed = 1, max.iter = 30)
fit$metadata$constraint_source
#> [1] "explicit_distances"

Unordered pairs must be unique and form a connected network. Targets must be positive and finite; inverse-squared weighting cannot represent zero distances. Optional count_i and count_j give nonnegative endpoint multiplicities; both default to one. They are divided by squared target distance for updates. The defaults give symmetric sparse stress. Supplying asymmetric multiplicities supports the pivot-region weighting used by the graph-based sparse method, but the caller must compute those multiplicities consistently with its metric. The returned metadata preserves the targets and multiplicities.

For Fermat geometry, local constraints also need shortest-path targets: powered direct-edge lengths can exceed them. The dgraphs Fermat vignette includes an implicit-distance example and a helper constructing exact local and pivot constraints. The graph-based interface remains unchanged and continues to use supplied edge lengths for its local targets. Sparse fitting does not certify unmeasured distances or recover full stress; evaluate on independent pairs.