Separate legacy 2D rescaling from dimension-generic normalisation

[t-uda]

Problem

EllipseCloud.rescale() is currently restricted to two-dimensional ellipse clouds. The original proposal attempted to generalise it using the volume-equivalent radius

$$ r_{\mathrm{vol}}(\Sigma) = \det(\Sigma)^{1/(2d)}. $$

However, this does not reduce to the existing two-dimensional implementation.

For aggregated semi-axis lengths (A \ge B), the current implementation uses

$$ s_{\mathrm{legacy}}=\frac{A^2}{B}, $$

whereas the volume-equivalent radius is (\sqrt{AB}). These coincide only in the isotropic case (A=B). Replacing the current formula would therefore silently change all existing two-dimensional results.

Interpretation of the current rule

For an ellipse with semi-axes (a\ge b), the quantity

$$\frac{a^2}{b} $$

is the radius of curvature at a vertex on the minor axis. Thus, the existing implementation should be regarded as a curvature-based two-dimensional calibration, rather than a volume-based normalisation.

This interpretation is particularly relevant when ellipses are estimated from local covariance matrices: the ratio (a^2/b) can remove the leading dependence on neighbourhood width and retain information related to the curvature radius of an underlying sampled curve.

Why there is no automatic arbitrary-dimensional extension

In dimensions greater than two, several inequivalent scalar size measures are available:

• volume-equivalent radius;
• maximum or minimum principal curvature radius;
• a mean-curvature-based radius;
• a Gauss–Kronecker-curvature radius.

For an ellipsoid with semi-axes

$$ s_1\ge \cdots \ge s_d, $$

the geometric mean of the principal curvature radii at a vertex on the shortest axis is

$$ R_K = \frac{ {( \prod_{j\lt d} s_j )}^{2/(d-1)} }{ s_d }. $$

This reduces to (s_1^2/s_2) when (d=2), but it is naturally associated with a hypersurface interpretation. It is not necessarily appropriate for lower-dimensional manifolds embedded in a higher-dimensional ambient space, such as curves in (\mathbb R^3).

Consequently, ambient dimension alone is insufficient to determine a canonical generalisation.

Proposed API direction

Separate the existing behaviour from genuinely dimension-generic normalisations.

1. Preserve the current two-dimensional behaviour

The existing rule remains available without numerical changes:

cloud.rescale(method="median", measure="legacy")

For two-dimensional data, measure="legacy" computes the current per-axis aggregate followed by

scale = major**2 / minor

The default behaviour for existing two-dimensional calls must remain backward compatible.

2. Add an explicit volume-based measure

Provide a dimension-generic option:

cloud.rescale(method="median", measure="volume")

For each covariance matrix (\Sigma_i), define

$$ r_i=\det(\Sigma_i)^{1/(2d)}, $$

and aggregate the scalar values (r_i) using the requested statistic:

$$ s=\mathrm{median}_i r_i $$

or the corresponding arithmetic mean.

This option represents volume-equivalent radius normalisation and must not be presented as equivalent to the existing two-dimensional rule.

3. Defer curvature-based higher-dimensional generalisation

Do not introduce (R_K), or another curvature-derived scalar, as an automatic default in this issue. Such an option requires a documented geometric assumption, including the intended intrinsic dimension of the sampled structure.

A separate follow-up issue may investigate curvature-based normalisation for hypersurface data.

Diagnostics

RescaleDiagnostics must report enough information to identify the applied convention. At minimum, diagnostics should contain:

• the scalar scale factor;
• the aggregation method;
• the selected measure;
• the relevant pre- and post-rescaling summaries.

Legacy per-axis summaries remain appropriate for measure="legacy". Volume-based rescaling should instead report summaries of the scalar values (\det(\Sigma_i)^{1/(2d)}).

Positive-definiteness requirements

All implemented measures assume positive-definite covariance matrices. In higher dimensions, local covariance matrices are necessarily singular when the neighbourhood contains too few independent samples.

The implementation must therefore detect non-positive eigenvalues and raise a clear error. Covariance regularisation is outside the scope of this issue and should not be applied silently.

Acceptance criteria

• Existing two-dimensional calls retain exactly the current numerical behaviour.
• The formula (A^2/B) is documented explicitly.
• A dimension-generic volume-based option is available only under an explicit measure name.
• The implementation does not claim that volume-based scaling generalises the existing rule.
• Diagnostics distinguish legacy and volume-based measures.
• Singular or non-positive-definite covariance matrices produce an informative error.
• Tests cover isotropic and anisotropic 2D cases, at least one 3D volume-based case, and singular covariance input.

Out of scope

• Choosing a canonical curvature scale for arbitrary intrinsic dimensions.
• Estimating intrinsic dimension.
• Regularising singular local covariance matrices.
• Changing the filtration or tangency solver itself.

[t-uda]

Triage note (2026-07): sequencing decision from the backlog triage —

1. PR Add RescaleDiagnostics to return detailed rescaling information #119 (rescale diagnostics, feat: add diagnostics for EllipseCloud.rescale() #105) merges first; it computes RescaleDiagnostics summaries from the current 2D-only internals.
2. This issue's n-dim generalisation must therefore also make the diagnostics dimension-generic (same method, same PR): pre/post summaries computed from det(Σ)^(1/2d)-based scales rather than 2D ell_scales.
3. 良い感じのグローバルスケーリング可視化 #125 (global-scaling visualisation) has an open design question about whether it belongs in rescale() — decide it after this lands; do not block on it.

The reference implementation in the body (median geometric mean of semi-axes, verified on 3D galaxy data) is the agreed approach. Implementing agents: read contrib/agent-skill/SKILL.md before writing tests/experiments, and base the PR on the release branch per CONTRIBUTING.md.

[t-uda]

The proposed formula does not reduce to the current 2D implementation. They coincide only when every ellipse is a circle.

Current 2D implementation

scales     = sqrt(eigvalsh(cov))       # (N, 2), ascending -> [minor, major]
ell_scales = median(scales, axis=0)    # per-axis median -> [B, A]
ell_scale  = A**2 / B

Proposed: s = median_i (a_i b_i)^{1/2}.

These differ in two independent ways: the median is taken per axis and then combined, rather than combined per ellipse and then medianed; and the combination rule is A²/B, not √(AB). For a cloud of identical ellipses with semi-axes a ≥ b, the two agree iff √(ab) = a²/b, i.e. iff a = b.

A quick numerical check (N=5000, b ~ U(0.3, 1.0), a = b · U(1, 4)):

current 2D implementation : 3.344713
proposed formula at d=2   : 0.972858     (ratio 3.44)
circles only              : identical

So merging the proposal as written would silently change every existing 2D result, including the experiments behind the published paper.

A geometric reading of the current rule. For an ellipse with semi-axes a ≥ b, a²/b is the radius of curvature at the vertex on the minor axis — the osculating circle at the flattest point of the boundary. That is a natural scale for a tangency-based filtration, and it is a different notion from the volume-equivalent radius √(ab) that the proposal uses. If the current rule was arrived at empirically as the best match against the Euclidean side, the curvature reading may be why it works: contact events are governed by curvature, not by enclosed volume.

A generalisation that does reduce exactly. For an ellipsoid with semi-axes s_1 ≥ s_2 ≥ … ≥ s_d, the principal curvature radii at the vertex on the shortest axis are s_j² / s_d for j < d. Their geometric mean is

R = (∏_{j<d} s_j)^{2/(d-1)} / s_d

Applying it to the per-axis medians S_j (keeping the current order of operations):

• d = 2 → S_1² / S_2, i.e. exactly the current implementation.
• d = 3 → S_1 S_2 / S_3.

def rescale(self, *, method="median") -> float:
    d = self.n_dim
    scales = numpy.sqrt(numpy.linalg.eigvalsh(self.cov))      # (N, d), ascending
    agg = numpy.median if method == "median" else numpy.average
    S = agg(scales, axis=0)[::-1]                             # descending
    s = float(numpy.prod(S[:-1]) ** (2.0 / (d - 1)) / S[-1])
    self.cov /= s ** 2
    self.coef *= s ** 2
    return s

Both candidates are defensible; they just mean different things (curvature radius vs volume-equivalent radius). The one above is the one that preserves current behaviour. Whichever is adopted, it would help to state the closed-form formula in the rescale docstring, since the choice is not recoverable from the code without working out the algebra.

[t-uda]

The issue body has been rewritten following the mathematical correction above.

The original proposal incorrectly stated that the volume-equivalent radius (\det(\Sigma)^{1/(2d)}) reduces to the existing two-dimensional formula. The revised scope now separates:

1. the backward-compatible two-dimensional rule (A^2/B);
2. an explicit, dimension-generic volume-based normalisation; and
3. possible higher-dimensional curvature-based normalisations, which require additional geometric assumptions and are deferred.

The previous statement that the determinant-based formula was the agreed generalisation should therefore be regarded as superseded.