DA: uncertainty is numerically unstable when n_landmarks >= n_combined (the GPSettings default on <5000 cells) #32

Summary

DifferentialAbundance with n_landmarks >= n_combined (which is what any GPSettings object routes kompot.da() into on a dataset under 5000 cells, since GPSettings.n_landmarks defaults to 5000) produces an uncertainty that is numerically unstable: the point estimates agree between the two orientations of the same comparison to 1e-11, and with the full GP to 1e-5, while the uncertainty and the z-scores differ between orientations by a median of 97% and 82% across every cell, and differ from the full-GP uncertainty by a median of 49%.

The inducing set here contains the training set, so mathematically the sparse and the full GP posterior variances must coincide. A 30x per-cell difference is therefore numerical, not a modelling choice: consistent with an ill-conditioned kernel over the whole union in the sparse uncertainty path. I have not investigated mellon's sparse uncertainty internals; this is a measurement with a stated hypothesis, not a diagnosis.

How it is reached

• kompot/anndata/differential_abundance.py:104 uses n_landmarks=None only when gp is None; any GPSettings carries n_landmarks=5000 (kompot/settings.py:87).
• DifferentialAbundance.fit then calls compute_landmarks(np.vstack([X1, X2]), gp_type="fixed", n_landmarks=5000, ...); mellon 1.7.1 parameters.py:280-288 returns the whole union as the landmark array when n_landmarks >= n, and both per-condition estimators are built with gp_type="fixed" on it.
• This is the path introduced by the 0.8.0 entry "Shared landmarks across conditions in DA" (CHANGELOG.md line 225 at c907da9). The k-means path (n_landmarks < n_combined) is well conditioned in the same measurement (orientation gap on uncertainty 1.7e-03 at 40 landmarks).

So the affected population is exactly the documented configuration path on small data: tutorials, tests, and any user who sets a seed via GPSettings(random_state=...).

Measurement

Synthetic data, seed 0: two 2-D Gaussian clouds of 400 cells each at x = -2.5 and +2.5, sd 0.45. DifferentialAbundance(n_landmarks=5000, random_state=0), fit (X1, X2) and, as the mirror, (X2, X1); both predicted on the same 800 stacked points. Relative discrepancy is |u - v| / max(|u|, |v|) per cell. Reference: n_landmarks=None (mellon full GP per condition). Tree c907da9, mellon 1.7.1, jax 0.10.0, numpy 2.4.4, CPU.

Uncertainty ranges: fixed-landmark run A [6.9e-03, 4.7e-01], run B [7.1e-03, 4.4e-01], full GP [1.1e-01, 4.3e-01]. No negative variances. Stored landmarks_info reports n_landmarks: 800.

The point-estimate rows localise the defect to the uncertainty path alone: densities and fold changes are reproduced across orientations and against the full GP; only the posterior variance is not. On this fixture the direction labels still exchange 800/800 because every cell clears both thresholds; on data near the thresholds, significance calls would depend on which condition is listed first.

Through the public API the same thing reads: kompot.da(adata, "cond", "A", "B", gp=GPSettings(random_state=0)) versus kompot.da(adata, "cond", "B", "A", gp=GPSettings(random_state=0)) gives ..._lfc agreement to 3.9e-11 and ..._lfc_zscore disagreement up to 0.84 relative; with gp=None every column is bit-identical.

Suggested checks

1. Reproduce with the snippet shape above; then vary n_landmarks across n_combined to see whether the instability is specific to n_landmarks >= n (union as landmarks) or grows continuously as the landmark count approaches n.
2. Compare mellon's sparse posterior-variance computation against the full GP when the inducing set equals the training set; a jitter or rank-reduction step is the likely place.
3. Decide whether n_landmarks >= n_combined should fall back to the full GP per condition rather than a fixed landmark set that is the union, or whether the anndata entry point should treat GPSettings.n_landmarks for DA the way it treats gp=None.

Found while reviewing PR #22; not caused by it.