Cross-Validating Polynomial Shift Surfaces
Fit error always falls as polynomial degree rises, so it cannot choose the degree — only error on points the fit has never seen can. This guide, part of polynomial shift algorithms for regional adjustments, implements k-fold cross-validation for a shift surface: how to fold a spatial control set without leaking information between folds, what statistic to select on, and how to report the result so a reviewer can see the degree was chosen rather than assumed.
Why Spatial Folds Are Different
Standard k-fold cross-validation assigns observations to folds at random, which assumes they are independent. Control points are not: two monuments two hundred metres apart see nearly the same distortion surface, so a random split puts near-duplicates on both sides of the divide and the held-out error comes out optimistically small. The model looks better at extrapolation than it is, which is precisely the property being tested.
Figure — random folds leak neighbours across the divide; blocked folds do not.
Spatially blocked folds fix it. Partition the area into contiguous blocks and assign whole blocks to folds, so a held-out point has no near neighbour in the training set. The resulting error is larger and honest, and the gap between the random-fold and blocked-fold estimates is itself informative: a wide gap means the control is clustered and the surface is being interpolated between clusters rather than modelled.
Complete Runnable Implementation
from __future__ import annotations
from dataclasses import dataclass
import numpy as np
def design_matrix(x: np.ndarray, y: np.ndarray, degree: int) -> np.ndarray:
"""Polynomial terms up to `degree` in two variables, centred inputs assumed."""
cols = [x ** i * y ** j
for total in range(degree + 1)
for i in range(total + 1)
for j in [total - i]]
return np.column_stack(cols)
def spatial_folds(x: np.ndarray, y: np.ndarray, k: int,
seed: int = 0) -> list[np.ndarray]:
"""Contiguous spatial blocks assigned to k folds.
Blocks rather than random points: two monuments a few hundred metres apart see
almost the same surface, so splitting them across the divide leaks the answer
into the training set and makes the held-out error look better than it is.
"""
n = len(x)
if k < 2 or k > n:
raise ValueError(f"k must be between 2 and {n}")
side = int(np.ceil(np.sqrt(k * 2)))
xi = np.clip(((x - x.min()) / max(np.ptp(x), 1e-12) * side).astype(int), 0, side - 1)
yi = np.clip(((y - y.min()) / max(np.ptp(y), 1e-12) * side).astype(int), 0, side - 1)
block = xi * side + yi
rng = np.random.default_rng(seed)
unique = rng.permutation(np.unique(block))
assignment = {b: i % k for i, b in enumerate(unique)}
fold_of = np.array([assignment[b] for b in block])
return [np.flatnonzero(fold_of == f) for f in range(k)]
@dataclass(frozen=True)
class CrossValidation:
degree: int
fit_rmse_m: float
holdout_rmse_m: float
holdout_p95_m: float
n_terms: int
n_points: int
def describe(self) -> str:
return (f"degree {self.degree}: fit {self.fit_rmse_m * 1000:.1f} mm, "
f"hold-out {self.holdout_rmse_m * 1000:.1f} mm "
f"({self.n_terms} terms, {self.n_points} points)")
def cross_validate(x: np.ndarray, y: np.ndarray, shift: np.ndarray,
degree: int, k: int = 5) -> CrossValidation:
"""K-fold spatially blocked cross-validation of one polynomial degree."""
xc, yc = x - x.mean(), y - y.mean() # centre: conditioning, not taste
a_full = design_matrix(xc, yc, degree)
if a_full.shape[1] > len(x) // 2:
raise ValueError(
f"degree {degree} needs {a_full.shape[1]} terms for {len(x)} points; "
f"the fit would be describing the control, not the surface"
)
coef, *_ = np.linalg.lstsq(a_full, shift, rcond=None)
fit_rmse = float(np.sqrt(np.mean((a_full @ coef - shift) ** 2)))
errors: list[float] = []
for hold in spatial_folds(xc, yc, k):
train = np.setdiff1d(np.arange(len(x)), hold)
c, *_ = np.linalg.lstsq(design_matrix(xc[train], yc[train], degree),
shift[train], rcond=None)
pred = design_matrix(xc[hold], yc[hold], degree) @ c
errors.extend(np.abs(pred - shift[hold]).tolist())
err = np.array(errors)
return CrossValidation(degree, fit_rmse, float(np.sqrt(np.mean(err ** 2))),
float(np.percentile(err, 95)), a_full.shape[1], len(x))
def choose_degree(x: np.ndarray, y: np.ndarray, shift: np.ndarray,
degrees=(1, 2, 3, 4)) -> CrossValidation:
"""Pick the degree with the lowest hold-out RMSE, ties going to the simpler."""
results = [cross_validate(x, y, shift, d) for d in degrees]
return min(results, key=lambda r: (round(r.holdout_rmse_m, 6), r.degree))Parameter Reference
| Name | Type | Units | Note |
|---|---|---|---|
x, y |
np.ndarray |
m or degrees | Centred internally for conditioning |
shift |
np.ndarray |
m or arc-seconds | One component at a time |
k |
int |
— | 5 is a reasonable default for tens of points |
degree |
int |
— | Refused when terms exceed half the point count |
holdout_p95_m |
float |
m | Reported with RMSE; the tail is what fails |
Worked Example
import numpy as np
rng = np.random.default_rng(5)
x = rng.uniform(-1.0, 1.0, 60)
y = rng.uniform(-1.0, 1.0, 60)
truth = 0.20 + 0.03 * x + 0.02 * y + 0.015 * x * y - 0.010 * x ** 2
shift = truth + rng.normal(0.0, 0.004, 60)
for d in (1, 2, 3, 4):
print(cross_validate(x, y, shift, d).describe())
# degree 1: fit 8.0 mm, hold-out 8.7 mm (3 terms, 60 points)
# degree 2: fit 3.9 mm, hold-out 4.5 mm (6 terms, 60 points)
# degree 3: fit 3.7 mm, hold-out 4.8 mm (10 terms, 60 points)
# degree 4: fit 3.5 mm, hold-out 5.9 mm (15 terms, 60 points)The fit column falls monotonically and says nothing; the hold-out column bottoms out at degree two — which is the degree the data was generated at — and then climbs. That climb is the model fitting the noise, and it is invisible in the fit statistic.
Figure — fit error falls forever; hold-out error is the one with a minimum.
Validation Check
def assert_folds_are_spatial(x: np.ndarray, y: np.ndarray, k: int = 5) -> None:
"""Held-out points must not have a training neighbour on top of them."""
folds = spatial_folds(x, y, k)
for hold in folds:
train = np.setdiff1d(np.arange(len(x)), hold)
d = np.hypot(x[hold][:, None] - x[train][None, :],
y[hold][:, None] - y[train][None, :])
assert d.min() > 0.0, "a held-out point coincides with a training point"Common Mistakes
Random folds on spatial data. The commonest and the most flattering: it reports a hold-out error close to the fit error, and the conclusion is that a high degree generalises well when it has merely been shown its own neighbours.
Selecting on RMSE alone. RMSE hides the tail, and a surface that is excellent almost everywhere and metres out in one corner has a respectable RMSE. Report the 95th percentile alongside it and look at where the worst folds are.
Choosing a degree the control cannot support. Fifteen terms fitted to twenty points will fit them beautifully and behave arbitrarily between them. The guard above refuses when the term count exceeds half the point count, which is conservative and still permissive compared with the practical guidance in choosing polynomial degree for regional fits.
Reporting the Selection
A degree chosen by cross-validation is defensible; a degree chosen by cross-validation and then reported as a bare number is indistinguishable from one chosen by habit. Four items make the selection reviewable.
Figure — four items that make a chosen degree reviewable rather than asserted.
The table of degrees tried, with fit and hold-out error for each — this is the evidence, and it is three lines of output. The fold strategy, stated as spatial or random, because the two produce materially different numbers and a reader cannot tell which was used from the result. The control count and extent, since a degree that is right for sixty points over a county is wrong for fifteen over the same area. And the extent of validity, which is the convex hull of the control and is the boundary beyond which the surface must not be used.
Publishing that alongside the coefficients also settles a question that otherwise recurs on every project: whether a higher degree was tried. It usually was, and the fact that it was rejected on hold-out error is the interesting part of the record.
Frequently Asked Questions
How many folds should I use?
Five is a reasonable default. More folds mean more training data per fit and a less stable estimate per fold; fewer mean the training set is noticeably smaller than the real one, which biases the estimate pessimistically. With very few control points, leave-one-block-out is the honest limit — and a control set that small is an argument against a polynomial rather than for a particular k.
Should each shift component be cross-validated separately?
Yes. The latitude and longitude shift surfaces have different structure, and nothing requires the same degree to be best for both. Fitting them at a common degree is a defensible simplification if it is recorded; assuming it without checking is not.
What if the best degree is one?
Then the distortion between the frames is essentially a plane over this area, and a similarity or affine transformation is the more interpretable model with the same fit — as implementing affine transformations for local grids sets out. A degree-one polynomial is an affine transformation written awkwardly.
Does cross-validation say anything about extrapolation?
Only indirectly, and it is worth being explicit about the limit: every held-out point still sits inside the convex hull of the full control set. Cross-validation measures interpolation quality. Beyond the hull a polynomial diverges regardless of what any fold said, which is why the model should be clamped to its published extent.