Spatially constrained MDS¶

The default backend uses weighted metric multidimensional scaling (MDS) with a geographic coordinate constraint. It preserves the central idea of classical MDS—representing a dissimilarity matrix as Euclidean distances—but optimizes raw metric stress so pair weights and fixed reference coordinates can enter the objective directly.

Classical versus metric MDS

“Classical MDS” normally means the Torgerson–Gower double-centering and eigendecomposition solution, which minimizes strain. SpatialDeform's backend="mds" instead minimizes weighted stress with regularized SMACOF. This is the constrained extension described below. The backend="isomap" implementation uses classical scaling after computing graph shortest-path distances.

Objective¶

Let:

\(G_i\in\mathbb{R}^2\) be the original projected coordinate of node \(i\);
\(X_i\in\mathbb{R}^2\) be its deformed coordinate;
\(T_{ij}\) be the shortest-path travel cost between nodes \(i\) and \(j\);
\(s>0\) convert cost units into map units;
\(\delta_{ij}=sT_{ij}\) be the target map distance;
\(g_{ij}=\lVert G_i-G_j\rVert_2\) be original geographic distance; and
\(w_{ij}=f(g_{ij})\) be a non-negative spatial pair weight.

For \(n\) nodes, the implementation minimizes:

\[ \boxed{ \mathcal{L}(X)= \underbrace{\sum_{i<j}w_{ij} \left(\lVert X_i-X_j\rVert_2-\delta_{ij}\right)^2}_{ \text{weighted travel-cost stress}} + \underbrace{\alpha\sum_{i=1}^{n}\lVert X_i-G_i\rVert_2^2}_{ \text{geographic coordinate constraint}} } \]

The first term asks Euclidean distance on the deformed map to represent network travel cost. The second penalizes displacement from observed coordinates, preserving recognizability and fixing translation, rotation, and orientation.

The effective coordinate-constraint coefficient is normalized to the weight matrix:

\[ \alpha=\lambda_{geo}\,\bar w, \qquad \bar w=\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}, \]

where \(\lambda_{geo}\) is the public geo_weight parameter. This normalization keeps its practical scale more stable when switching between uniform, Gaussian, and inverse pair weights.

Spatial pair weights¶

The weights depend only on original coordinates, not on the changing deformed coordinates:

\[ w_{ij}= \begin{cases} 1, & \text{uniform},\\[3pt] \exp\!\left[-\frac{1}{2}\left(\frac{g_{ij}}{b}\right)^2\right], & \text{Gaussian},\\[6pt] \left(1+\frac{g_{ij}}{b}\right)^{-1}, & \text{inverse}, \end{cases} \qquad w_{ii}=0. \]

Here \(b\) is spatial_bandwidth in working-CRS units. If it is omitted, the median positive \(g_{ij}\) is used. Gaussian weighting concentrates the loss on nearby node pairs; inverse weighting gives distant pairs a longer-tailed influence.

Value	Weight	Cartographic effect
`uniform`	\(1\)	every pair contributes equally
`gaussian`	\(\exp[-\tfrac12(g_{ij}/b)^2]\)	strongest local emphasis
`inverse`	\((1+g_{ij}/b)^{-1}\)	local emphasis with a longer tail

local_model = SpatialDeformer(
    backend="mds",
    spatial_weight="gaussian",
    spatial_bandwidth=5_000,  # metres in a metric CRS
    geo_weight=0.1,
)

Cost-to-map scale¶

With scale="auto", \(s\) is resolved before optimization by an unweighted least-squares fit from travel costs to original geographic distances:

\[ s^*=\underset{s}{\operatorname{argmin}} \sum_{i<j}(g_{ij}-sT_{ij})^2 = \frac{\sum_{i<j}g_{ij}T_{ij}} {\sum_{i<j}T_{ij}^2}. \]

This removes arbitrary units such as seconds versus minutes. It is deliberately separate from spatial pair weighting. Reuse a fixed scale_ to retain global contraction and expansion across scenarios; see Multiple scenarios.

SMACOF update¶

Define the weighted graph Laplacian

\[ V=\operatorname{diag}(W\mathbf{1})-W \]

and, for the current configuration \(X^{(k)}\), define the SMACOF majorization matrix \(B(X^{(k)})\) by

\[ b_{ij}= \begin{cases} -w_{ij}\dfrac{\delta_{ij}} {\lVert X_i^{(k)}-X_j^{(k)}\rVert_2}, & i\ne j \text{ and the denominator is non-zero},\\[8pt] 0, & i\ne j\text{ and the denominator is zero},\\[3pt] -\sum_{\ell\ne i}b_{i\ell}, & i=j. \end{cases} \]

Each iteration solves

\[ X^{(k+1)}= \left(V+(\alpha+\varepsilon)I\right)^{-1} \left[B(X^{(k)})X^{(k)}+\alpha G\right], \]

where \(\varepsilon\) is a tiny numerical ridge used only to make the linear system stable. The optimizer starts at \(X^{(0)}=G\). When geo_weight=0, the updated configuration is recentered to the geographic centroid because an unanchored distance objective is translation-invariant.

For uniform pair weights, \(V=nI-\mathbf{1}\mathbf{1}^{\mathsf T}\). SpatialDeform therefore applies the equivalent Sherman–Morrison rank-one solve instead of a dense Cholesky factorization. See Performance.

Iterations stop when relative pair-stress change falls below tol or when max_iter is reached.

Geographic regularization¶

geo_weight controls the balance:

flexible = SpatialDeformer(geo_weight=0.01)
balanced = SpatialDeformer(geo_weight=0.1)
conservative = SpatialDeformer(geo_weight=1.0)

Smaller values prioritize travel-cost fidelity and allow greater distortion.
Larger values preserve geography but leave more residual stress.
0 removes the geographic penalty; the implementation still centers the solution, but rotation and orientation are less meaningful cartographically.

Tune this parameter according to the intended interpretation, not only the smallest stress.

Reading stress¶

stress_ reports only the data-fit term from the boxed objective:

\[ \operatorname{stress}(X)= \sum_{i<j}w_{ij} \left(\lVert X_i-X_j\rVert_2-\delta_{ij}\right)^2. \]

The geographic penalty is not included in this reported diagnostic. normalized_stress_ is

\[ \operatorname{normalized\_stress}(X)= \frac{\operatorname{stress}(X)} {\sum_{i<j}w_{ij}\delta_{ij}^2}, \]

which is more useful when comparing hyperparameters on the same network.

The estimator’s score() returns negative normalized stress so that higher is better, following scikit-learn model-selection convention.

for weight in (0.01, 0.1, 1.0):
    model = SpatialDeformer(geo_weight=weight).fit(edges)
    print(weight, model.normalized_stress_, model.score())

Stress alone does not measure the coordinate penalty, geographic readability, line crossing, or topological ambiguity. Always inspect the result.