Landmark Multi-Dimensional Scaling

• LMDS is a computational method that approximates classical MDS by using a small set of landmark points to reduce the full eigendecomposition cost.
• It employs the Nyström approximation to extend the embedding from landmarks to non-landmark points while preserving geometric properties.
• The method balances efficiency and accuracy, with landmark selection strategies, such as random sampling or k-means, directly impacting performance.

Landmark multi-dimensional scaling (LMDS), also called landmark MDS, is a landmark-based approximation to classical multidimensional scaling in which only $m\ll n$ “landmark” points are used to avoid the full $n\times n$ eigendecomposition. In the survey literature, it is presented as a Nyström approximation for big-data embedding; in the localization literature, it is described as a leading method for embedding new points into an existing coordinate system from observed distance information, and as Gower’s method from the 1960s viewed through the Nyström lens (Ghojogh et al., 2020, Ouyang et al., 14 Sep 2025). Its defining feature is that landmarks are embedded by classical MDS on a small block, while the remaining points are placed by an out-of-sample extension derived from cross-distances to the landmarks.

1. Classical MDS as the starting point

Classical MDS begins from $n$ points $x_1,\dots,x_n\in\mathbb R^d$ and their squared pairwise-distance matrix

$$D_{ij}=\|x_i-x_j\|^2,\quad D\in\mathbb R^{n\times n}.$$

With the centering matrix

$H=I_n-\tfrac1n\,\mathbf1\mathbf1^T,$

the associated Gram matrix is

$K=-\tfrac12\,H\,D\,H\in\mathbb R^{n\times n}.$

When the entries of $D$ arise from Euclidean distances, $K$ is symmetric and positive semidefinite. Classical MDS then computes

$K=V\,\Lambda\,V^T,\quad \Lambda=\diag(\lambda_1,\dots,\lambda_n),\;\lambda_1\ge\cdots\ge\lambda_n\ge0,$

and forms a $n\times n$0-dimensional embedding

$n\times n$1

Equivalently, the embedded point $n\times n$2 has coordinates $n\times n$3 (Ghojogh et al., 2020).

The computational cost of this construction is the main motivation for LMDS. Forming $n\times n$4 costs $n\times n$5, forming $n\times n$6 costs $n\times n$7, the full eigendecomposition costs $n\times n$8, and storing $n\times n$9 requires $n$0 memory. LMDS keeps the classical MDS geometry but replaces the full spectral step by a low-rank landmark approximation.

2. Nyström construction of LMDS

LMDS proceeds by selecting $n$1 landmarks and permuting indices so that the first $n$2 correspond to the landmark set $n$3, with the remaining $n$4 points in $n$5. The Gram matrix is partitioned as

$n$6

Only $n$7 and $n$8 are computed explicitly. Under the assumption that $n$9 has approximate rank $x_1,\dots,x_n\in\mathbb R^d$0, the Nyström approximation is

$x_1,\dots,x_n\in\mathbb R^d$1

where $x_1,\dots,x_n\in\mathbb R^d$2 is the Moore–Penrose inverse. The stated intuition is that once the landmark–landmark block and the landmark off-diagonal blocks are known, the remainder is determined in order to preserve positive semidefiniteness (Ghojogh et al., 2020).

The spectral computation is then confined to the small block: $x_1,\dots,x_n\in\mathbb R^d$3 Keeping the top $x_1,\dots,x_n\in\mathbb R^d$4 modes gives $x_1,\dots,x_n\in\mathbb R^d$5 and $x_1,\dots,x_n\in\mathbb R^d$6. The landmark coordinates are

$x_1,\dots,x_n\in\mathbb R^d$7

while the non-landmarks are embedded by the Nyström out-of-sample formula

$x_1,\dots,x_n\in\mathbb R^d$8

Collecting $x_1,\dots,x_n\in\mathbb R^d$9 yields exactly the same embedding one would obtain by applying full classical MDS to $$D_{ij}=\|x_i-x_j\|^2,\quad D\in\mathbb R^{n\times n}.$$0 (Ghojogh et al., 2020).

3. Algorithmic workflow, complexity, and landmark choice

The standard LMDS workflow is explicit. One chooses landmarks $$D_{ij}=\|x_i-x_j\|^2,\quad D\in\mathbb R^{n\times n}.$$1 with $$D_{ij}=\|x_i-x_j\|^2,\quad D\in\mathbb R^{n\times n}.$$2, sets $$D_{ij}=\|x_i-x_j\|^2,\quad D\in\mathbb R^{n\times n}.$$3, computes distances from every point to every landmark, forms the landmark distance block $$D_{ij}=\|x_i-x_j\|^2,\quad D\in\mathbb R^{n\times n}.$$4 and cross-distance block $$D_{ij}=\|x_i-x_j\|^2,\quad D\in\mathbb R^{n\times n}.$$5, double-centers them to obtain

$$D_{ij}=\|x_i-x_j\|^2,\quad D\in\mathbb R^{n\times n}.$$6

eigendecomposes $$D_{ij}=\|x_i-x_j\|^2,\quad D\in\mathbb R^{n\times n}.$$7, computes $$D_{ij}=\|x_i-x_j\|^2,\quad D\in\mathbb R^{n\times n}.$$8, computes $$D_{ij}=\|x_i-x_j\|^2,\quad D\in\mathbb R^{n\times n}.$$9, and concatenates the result as $H=I_n-\tfrac1n\,\mathbf1\mathbf1^T,$0 (Ghojogh et al., 2020).

The complexity profile is the principal practical advantage. Computing all distances to landmarks costs $H=I_n-\tfrac1n\,\mathbf1\mathbf1^T,$1, centering and forming kernel subblocks costs $H=I_n-\tfrac1n\,\mathbf1\mathbf1^T,$2, eigendecomposition of the $H=I_n-\tfrac1n\,\mathbf1\mathbf1^T,$3 block costs $H=I_n-\tfrac1n\,\mathbf1\mathbf1^T,$4, and extension to the remaining points costs $H=I_n-\tfrac1n\,\mathbf1\mathbf1^T,$5. The total time is therefore roughly

$H=I_n-\tfrac1n\,\mathbf1\mathbf1^T,$6

with memory $H=I_n-\tfrac1n\,\mathbf1\mathbf1^T,$7 rather than $H=I_n-\tfrac1n\,\mathbf1\mathbf1^T,$8. This is the sense in which LMDS replaces an intractable $H=I_n-\tfrac1n\,\mathbf1\mathbf1^T,$9 eigenproblem by an $K=-\tfrac12\,H\,D\,H\in\mathbb R^{n\times n}.$0 one plus matrix multiplications.

Landmark choice controls the accuracy–cost trade-off. The survey explicitly lists three strategies: uniform random sampling, $K=-\tfrac12\,H\,D\,H\in\mathbb R^{n\times n}.$1-means centers, and a farthest-point heuristic. Larger $K=-\tfrac12\,H\,D\,H\in\mathbb R^{n\times n}.$2 improves the rank-$K=-\tfrac12\,H\,D\,H\in\mathbb R^{n\times n}.$3 approximation but increases time and memory cost; the reported practical rule is that $K=-\tfrac12\,H\,D\,H\in\mathbb R^{n\times n}.$4–$K=-\tfrac12\,H\,D\,H\in\mathbb R^{n\times n}.$5 often suffices (Ghojogh et al., 2020).

4. Gower’s formulation and the geometry of LMDS

In the anchor-based formulation, let $K=-\tfrac12\,H\,D\,H\in\mathbb R^{n\times n}.$6 landmarks have squared-distance matrix

$K=-\tfrac12\,H\,D\,H\in\mathbb R^{n\times n}.$7

and let

$K=-\tfrac12\,H\,D\,H\in\mathbb R^{n\times n}.$8

If $K=-\tfrac12\,H\,D\,H\in\mathbb R^{n\times n}.$9 has rank $D$0 and eigendecomposition

$D$1

then the MDS embedding of the landmarks is

$D$2

with $D$3 and $D$4. Defining

$D$5

one has

$D$6

For a new point with observed squared distances $D$7 to the anchors, the Nyström viewpoint sets

$D$8

and yields the LMDS coordinate

$D$9

Gower’s classical derivation starts from the linear system

$K$0

which leads to the same expression. The explicit conclusion is that LMDS $K$1 Gower $K$2 Anderson–Robinson (with the length term dropped), and that LMDS is exactly the Nyström extension of classical MDS (Ouyang et al., 14 Sep 2025).

This formulation also makes the geometry explicit. LMDS can be viewed as minimizing the “angle-only” objective

$K$3

Its normal equations are

$K$4

The interpretation given in the localization literature is that LMDS tries to match inner products $K$5, hence “angle preservation,” but has no length term. Consequently, when the $K$6 contain large radial errors, LMDS may mis-scale the new point (Ouyang et al., 14 Sep 2025).

5. Weighted variants, least-squares relations, and global optimization

The 2025 analysis places LMDS inside a one-parameter family of objectives obtained by augmenting the angle-only criterion with a length-preservation term. With

$K$7

the family is

$K$8

This objective interpolates among several range-based embedding methods and makes explicit what each method minimizes (Ouyang et al., 14 Sep 2025).

Method	Weight	Characterization
LMDS / Gower	$K$9	Enforces $K=V\,\Lambda\,V^T,\quad \Lambda=\diag(\lambda_1,\dots,\lambda_n),\;\lambda_1\ge\cdots\ge\lambda_n\ge0,$0 exactly
T-LMDS	$K=V\,\Lambda\,V^T,\quad \Lambda=\diag(\lambda_1,\dots,\lambda_n),\;\lambda_1\ge\cdots\ge\lambda_n\ge0,$1	Total LMDS
LS	$K=V\,\Lambda\,V^T,\quad \Lambda=\diag(\lambda_1,\dots,\lambda_n),\;\lambda_1\ge\cdots\ge\lambda_n\ge0,$2	Least-squares embedding in MDS coordinates

A central result is that the classical least-squares estimator

$K=V\,\Lambda\,V^T,\quad \Lambda=\diag(\lambda_1,\dots,\lambda_n),\;\lambda_1\ge\cdots\ge\lambda_n\ge0,$3

is algebraically equivalent, after centering and rotation into MDS coordinates, to minimizing $K=V\,\Lambda\,V^T,\quad \Lambda=\diag(\lambda_1,\dots,\lambda_n),\;\lambda_1\ge\cdots\ge\lambda_n\ge0,$4 with $K=V\,\Lambda\,V^T,\quad \Lambda=\diag(\lambda_1,\dots,\lambda_n),\;\lambda_1\ge\cdots\ge\lambda_n\ge0,$5. The stated interpretation is therefore that LS is a regularized solution of LMDS, or, equivalently, LS $K=V\,\Lambda\,V^T,\quad \Lambda=\diag(\lambda_1,\dots,\lambda_n),\;\lambda_1\ge\cdots\ge\lambda_n\ge0,$6 T-LMDS regularized toward length preservation with weight $K=V\,\Lambda\,V^T,\quad \Lambda=\diag(\lambda_1,\dots,\lambda_n),\;\lambda_1\ge\cdots\ge\lambda_n\ge0,$7. Smaller $K=V\,\Lambda\,V^T,\quad \Lambda=\diag(\lambda_1,\dots,\lambda_n),\;\lambda_1\ge\cdots\ge\lambda_n\ge0,$8 places less emphasis on angle preservation, while larger $K=V\,\Lambda\,V^T,\quad \Lambda=\diag(\lambda_1,\dots,\lambda_n),\;\lambda_1\ge\cdots\ge\lambda_n\ge0,$9 approaches LMDS.

Although each $n\times n$00 is a nonconvex 4th-order polynomial, the same work shows that it can be globally optimized through a generalized trust-region method. By introducing $n\times n$01, the problem is recast as a quadratic objective under a quadratic equality constraint. A KKT analysis gives conditions

$n\times n$02

and the solution is obtained either from a closed-form hard-case treatment or by bisection on a scalar dual equation. The reported convergence is global and linear in $n\times n$03. The same trust-region machinery is extended to multiple-source localization with anchors, where block-coordinate updates solve a sequence of related generalized trust-region subproblems (Ouyang et al., 14 Sep 2025).

6. Scope, applications, and distinction from Bayesian landmark methods

LMDS is positioned in the survey literature as a method for large-scale visualization of high-dimensional data, fast approximate embedding when new points arrive, and preprocessing for clustering or classification on manifolds. The same survey places landmark MDS and landmark Isomap under the broader use of Nyström approximation for big-data embedding. In the localization literature, LMDS is treated as a method for embedding unknown points into a coordinate system determined by known landmarks or anchors, with both single-source and multiple-source settings receiving explicit formulations (Ghojogh et al., 2020, Ouyang et al., 14 Sep 2025).

It is important to distinguish classical LMDS from landmark-based Bayesian constructions. A direct comparison is given in the sparse Bayesian multidimensional scaling literature: classical LMDS takes $n\times n$04 landmark points, forms an $n\times n$05 dissimilarity submatrix, performs classical MDS by spectral decomposition in $n\times n$06, and then “triangulates” the remaining $n\times n$07 points by solving $n\times n$08 sets of linear equations using distances to the $n\times n$09 landmarks in $n\times n$10. By contrast, landmark sparse Bayesian MDS (L-sBMDS) is fully Bayesian, posits a joint likelihood for all $n\times n$11 points on a subsampled set of dissimilarities, provides posterior samples of all $n\times n$12 rather than only a point estimate, and performs inference for non-landmarks within MCMC rather than by post-hoc triangulation (Sheth et al., 2024).

This comparison clarifies the methodological boundary around LMDS. In its classical form, LMDS is a deterministic Nyström extension built on the geometry of classical MDS and a landmark-induced low-rank approximation. Later work shows that the same formula is Gower’s method in MDS coordinates and that it occupies one endpoint of a broader family balancing angle and length preservation. A plausible implication is that, in practice, the choice among LMDS, T-LMDS, and LS is best understood not as a choice among unrelated heuristics, but as a choice among explicit objective functions with different geometric emphases.