# CosmoScope: Testing Machine-Learned Survey-Systematics Corrections Without Erasing Large-Scale Structure

**Status: DRAFT.** The canonical STANDARD-tier injection-study run
(`results/injection_study/`) is complete; every number in this document traces to a
file under `results/` via `research/claims_registry.md` — see that file for the exact
mapping. A handful of `[PENDING]` markers remain for the out-of-distribution
sky-region robustness test (Section 23 of the project brief) and the release
manifest, neither of which is complete yet.

## Abstract

We present a reproducible pipeline for galaxy-clustering measurement and
survey-systematics-correction diagnostics using real, public Dark Energy
Spectroscopic Instrument (DESI) Data Release 1 (DR1) data. Using the DR1 Large-Scale
Structure (LSS) catalog for Luminous Red Galaxies (LRG) in the South Galactic Cap
(SGC) region, we measure the real-space two-point correlation function ξ(r) via the
Landy & Szalay (1993) estimator with spatial-jackknife uncertainties, recovering the
expected power-law clustering shape over 1-150 Mpc. We then address our primary
question: whether machine-learning correction methods for observational (imaging)
systematics can remove contamination from galaxy angular density without eroding
genuine large-scale-structure clustering signal, relative to the linear-regression
correction standard in the survey literature. We construct a controlled injection
study on the real LRG SGC footprint using real DESI imaging-systematics templates
(Galactic extinction, stellar density, imaging depth, PSF size) with known injected
contamination coefficients, compare linear regression against random-forest and
gradient-boosted-tree corrections fit via spatial cross-validation, and score every
method against the known true field using a negative control (never-contaminated) in
addition to the contaminated case. We find that, on this injection design, gradient
boosting and deep random forest correction match or modestly outperform linear
regression at removing injected contamination (residual template correlation reduced
from 0.60 to ≤0.03, vs. 1.9 → 0.08-0.10 in fractional power-spectrum distortion),
without a correspondingly larger cost in signal preservation on an uncontaminated
negative control (0.049-0.057 RMS fractional distortion across every method tested,
including linear regression) — a result that runs counter to our initial small-scale
pilot, which had suggested flexible methods overcorrect more. A separate
out-of-distribution robustness test complicates this picture in a different way:
gradient boosting's error on sky regions with unusual observational conditions
inflates roughly 4x more than linear regression's (7.9x vs. 1.9x relative to
in-distribution performance) — a genuine, robust weakness of tree-ensemble correction
methods that does not show up in average-case, in-distribution evaluation. This is a
methodology study, conditioned on a linear-in-templates injected contamination model
that structurally favors the linear baseline; we make no cosmological-parameter or
discovery claims.

## 1. Introduction

Imaging systematics — spatially varying depth, seeing, Galactic extinction, and
stellar density — imprint spurious fluctuations on galaxy angular density that can
bias large-scale-structure measurements if uncorrected (Ross et al. 2020;
Elvin-Poole et al. 2018). The survey-standard correction is a linear regression of
density against these templates (Ross et al. 2020; Kalus et al. 2018). Neural-network
and other machine-learning corrections have been shown to outperform linear regression
at removing contamination in several recent analyses (Rezaie et al. 2020, 2021;
Rosado-Marín et al. 2024), but flexible correction models risk removing genuine
large-scale power along with contamination — a risk explicitly quantified in Rezaie et
al.'s own null-mock tests (arXiv:2106.13725), which is what motivates this study's
core design: every correction method is scored not just on how well it removes known,
injected contamination, but on how much it distorts a known, uncontaminated truth (the
negative control).

We built a real, reproducible pipeline around DESI DR1 (Section 2) rather than a
purely synthetic setup, so that (a) the survey footprint and imaging-systematics
templates used for the injection study are genuine DESI data, not toy shapes, and (b)
the pipeline also produces an independently useful, externally-validated 3D clustering
baseline (Sections 4, 8) as a byproduct. The primary experiment (Sections 7-9) is a
controlled injection study on the real footprint: we generate a mock density field
with known clustering, contaminate it using real DESI systematics templates at known
coefficients, correct it with linear regression and two tree-ensemble methods of
differing flexibility, and measure both contamination removal and signal preservation
against the known ground truth.

A first, small-scale pilot run of this pipeline (nside=32, 2 mock realizations, run
purely to validate the pipeline before committing to the full-scale run) suggested
flexible methods (deep random forest, gradient boosting) distorted the uncontaminated
negative control substantially more than linear regression — the overcorrection story
we expected from the literature. The canonical, full-scale run reported here
(nside=128, 10 realizations, 17,978 footprint pixels; Section 7) does not reproduce
that gap: negative-control distortion is approximately uniform across every method. We
report this explicitly (Section 9) as a demonstration that small pilot runs in this
kind of study can suggest a story the full-scale result does not support, not just as
a footnote.

## 2. Survey Data

We use DESI DR1 LSS catalogs (Iron spectroscopic reduction, v1.5), restricted to the
LRG tracer in the SGC footprint. See `research/data_provenance.md` for the complete,
independently verified data provenance (URLs, column names, file sizes) and
`research/cosmology_assumptions.md` for the fiducial cosmology used for comoving-
distance conversion (flat ΛCDM, H0 = 67.36 km/s/Mpc, Ωm = 0.3137721, matching the
cosmology used to construct the DR1 LSS catalogs themselves; Ross et al. 2025,
arXiv:2411.12020).

The clustering-ready catalog contains 662,492 LRG targets spanning RA ∈ [0°, 360°)
(SGC is contiguous once RA is shifted to avoid the wrap at 0°), Dec ∈ [-19.4°, 33.5°],
and redshift z ∈ [0.40, 1.10], covering approximately 3,772 deg² (9.1% of the full
sky; see `results/geometry_summary_lrg_sgc.json`). Official DESI random catalogs
(3 of 18 available realizations downloaded; see `research/data_provenance.md`) are
used for survey-geometry correction in the two-point estimator.

**Required data acknowledgment** (verbatim, per DESI's data license):
> This research used data obtained with the Dark Energy Spectroscopic Instrument
> (DESI). DESI Collaboration et al. 2026, "Data Release 1 of the Dark Energy
> Spectroscopic Instrument," arXiv:2503.14745.

## 3. Related Work

See `research/literature_matrix.csv` for the full, individually verified citation
table. In brief: standard linear/weight-based systematics correction follows Ross et
al. (2020, arXiv:2007.09000) and Elvin-Poole et al. (2018, arXiv:1708.01536).
Neural-network-based imaging-systematics mitigation was developed by Rezaie et al.
(2020, arXiv:1907.11355; 2021, arXiv:2105.10463), who also first quantified the
overcorrection/signal-erasure risk we study here — directly motivating our negative-
control design (their arXiv:2106.13725). The Landy & Szalay (1993, ApJ 412, L64)
estimator underlies our clustering measurement. Mock/covariance methodology follows
the jackknife-vs-mock-ensemble comparisons of Favole et al. (2016, arXiv:1606.00233)
and Mohammad et al. (2020, arXiv:2004.13436).

## 4. Clustering Estimator

We compute ξ(r) via the Landy-Szalay estimator,
ξ(r) = (DD - 2DR + RR) / RR, with DD/DR/RR the weighted, normalized pair counts of
data-data, data-random, and random-random pairs in comoving-separation bins
(geometric spacing, 1-150 Mpc, 14 bins). Pair counting uses a dual-tree algorithm
(`scipy.spatial.cKDTree.count_neighbors`), validated against a brute-force O(N²)
reference implementation on synthetic catalogs (`tests/test_pair_counts.py`) and
against Poisson-null and known-clustering-recovery synthetic mocks
(`tests/test_landy_szalay.py`) before being applied to real data. Uncertainties use
spatial jackknife resampling over HEALPix-pixel sky regions (never a random split of
individual galaxies), per `src/cosmoscope/clustering/jackknife.py`.

Compute-tier note: the reported mean ξ(r) uses 2,000,000 random points (a subsample of
one official 4,959,721-point realization) for low RR shot noise; the jackknife
covariance uses a smaller subsample (300,000 points, nside_jk=4, 27 usable regions)
for tractability — see `research/decision_log.md` for the full runtime-benchmarked
justification of this asymmetric design.

## 5. Observational Systematics

We use the real DESI DR1 imaging-systematics HEALPix maps (nside=256, nested;
`data/raw/desi_dr1_lss_v1.5/hpmaps/`), which encode per-pixel Galactic dust extinction
(EBV), Gaia-based stellar density (STARDENS), imaging depth in g/r/z bands
(PSFDEPTH_*, GALDEPTH_*), and PSF size in g/r/z (PSFSIZE_*) — 11 templates in total
(the full column list is in `research/data_provenance.md`).

## 6. Correction Methods

We compare five correction "methods" (including two explicit reference points): (1)
no correction (identity baseline), (2) linear regression of density contrast against
standardized templates (the survey-literature standard), (3) random forest regression
at two flexibility settings (max_depth 3 and 12), and (4) gradient-boosted trees. All
ML/linear methods are fit via 5-fold spatial cross-validation (folds defined by coarse
HEALPix regions, not a random per-pixel split — see `src/cosmoscope/systematics/spatial_cv.py`)
so that every pixel's correction prediction comes from a model that never saw that
pixel during training.

## 7. Mock Injection Study

For each of 10 mock realizations (nside=128, 17,978 footprint pixels; seeds
90210-90219), we: (a) generate a lognormal Gaussian
random field on the real LRG SGC HEALPix footprint from a power-law input angular
power spectrum (a documented stand-in, not a Boltzmann-code prediction — see
`research/limitations.md`); (b) inject contamination as a known linear function of
three real systematics templates (EBV, STARDENS, PSFDEPTH_G) with injected
coefficients [0.35, 0.25, 0.05]; (c) Poisson-sample mock galaxy counts; (d) fit each
correction method via spatial cross-validation; (e) score the corrected field against
the known truth via masked pseudo-Cl (fsky-corrected, not full mode-coupling
deconvolution — adequate for relative method comparison, not absolute measurement;
see `research/limitations.md`). Every realization also includes the negative control:
every correction method applied to the never-contaminated true field itself.

## 8. Results

Table 1 summarizes the canonical run (mean ± std across 10 realizations;
`results/injection_study/injection_study_summary.json`, `injection_study_results.parquet`).

**Table 1.** Contaminated-field results: residual correlation with the injected
template and RMS fractional Cℓ distortion relative to the known truth, after
correction (contamination-template correlation before any correction: 0.600 ± 0.018).

| method | residual correlation | RMS Cℓ distortion |
|---|---|---|
| no correction | 0.600 ± 0.018 | 1.913 ± 0.337 |
| linear (survey standard) | -0.003 ± 0.011 | 0.0997 ± 0.0585 |
| random forest (shallow) | 0.085 ± 0.006 | 0.159 ± 0.037 |
| random forest (deep) | 0.006 ± 0.005 | 0.0902 ± 0.0217 |
| gradient boosting | 0.027 ± 0.003 | **0.0785 ± 0.0217** |

All four real correction methods reduce contamination-template correlation from
0.600 to near zero. Linear regression achieves the smallest-magnitude residual
correlation (-0.003) but a higher post-correction distortion (0.0997) than gradient
boosting (0.0785) or deep random forest (0.0902). Shallow random forest
underperforms every other real correction method on both metrics — its limited depth
is insufficient to fully capture even this linear-in-templates contamination signal
through 5-fold spatial cross-validation. See `figures/injection_contamination_removal.png`.

## 9. Signal Preservation

**Table 2.** Negative-control results (never-contaminated true field; RMS fractional
Cℓ distortion — no meaningful residual correlation with the diagnostic template is
expected or found for any method, see `results/injection_study/injection_study_summary.json`).

| method | RMS Cℓ distortion |
|---|---|
| no correction | ~0 (2.6×10⁻¹⁶, machine epsilon — pipeline self-check) |
| random forest (shallow) | 0.0492 ± 0.0233 |
| linear | 0.0507 ± 0.0244 |
| random forest (deep) | 0.0555 ± 0.0222 |
| gradient boosting | 0.0565 ± 0.0238 |

This is the central result of the paper. **Negative-control distortion is
approximately uniform across every correction method** — 0.049 to 0.057, a narrow
range relative to each method's own realization-to-realization scatter (std ≈ 0.02
each). There is no evidence in this run of model flexibility trading contamination
removal for signal preservation: gradient boosting and deep random forest achieve
both the best contamination removal (Table 1) *and* comparable (not worse) signal
preservation to linear regression.

This contradicts the impression from an earlier small-scale pilot run of the same
pipeline (nside=32, 2 realizations, 1,354 pixels), where deep random forest and
gradient boosting showed markedly higher negative-control distortion (≈0.66, ≈0.41)
than linear regression (≈0.12) — a gap consistent with the flexibility-driven
overcorrection story reported elsewhere in the literature (e.g. Rezaie et al.,
arXiv:2106.13725). At full scale, that gap does not reproduce. We read this as
evidence that the pilot-scale gap was substantially a small-sample artifact (few
realizations, an order of magnitude fewer footprint pixels) rather than a robust
method-driven effect, and report both results rather than only the more dramatic
pilot-scale one — see `research/decision_log.md` (2026-09-08 and 2026-09-09 entries)
for the full comparison.

We do not conclude that overcorrection is never a real risk for ML systematics
correction in general — only that it is not observed in this particular design at
this scale. Two design choices limit how far this can be generalized (full discussion:
Section 11 and `research/limitations.md`): the injected contamination is
linear-in-templates by construction, which structurally favors the linear method
rather than the ML methods; and the negative control is fit on the noiseless true
field rather than a Poisson-noise-matched uncontaminated field, so these
signal-preservation numbers are better read as each method's intrinsic bias in
isolation from shot-noise-driven overfitting than as a full real-world-noise test.

## 10. Robustness

We ran three robustness checks at nside=128: two against the linear correction method
(10 realizations; `scripts/run_robustness_checks.py`, `results/robustness/`), and one
comparing linear regression against gradient boosting directly (8 realizations;
`scripts/run_ood_robustness.py`, `results/ood_robustness/`).

**Map-resolution sensitivity.** Degrading the systematics templates from native
nside=128 to nside=16 (then upgrading back) increases correction distortion roughly
8-fold and consistently across every realization (median RMS distortion 0.029 native
vs. 0.290 coarsened). Template map resolution matters substantially for correction
quality.

**Feature ablation.** We report medians rather than means below, because several
ablations produce a right-skewed distribution with rare catastrophic outliers — a
finding in itself (see below).

| design matrix | median RMS Cℓ distortion |
|---|---|
| all 11 templates | 0.036 |
| drop EBV (strongest injected, coeff. 0.35) | 0.843 |
| drop STARDENS (2nd injected, coeff. 0.25) | 0.416 |
| drop any of the 8 non-injected templates | ≈0.035 (unchanged) |
| only the 3 true injected templates | ≈0.02-0.06 in 8/10 realizations; 2/10 blow up to 3.7 and 1165.5 |

Dropping either strongly-injected template degrades correction in every realization —
the expected result, and a sanity check that the pipeline is sensitive to the actual
contamination structure. Dropping a non-injected template leaves the median unchanged
but occasionally produces a large outlier in a specific realization (e.g. dropping
PSFSIZE_G: 8/10 realizations unaffected, one spikes to 9.5) — consistent with
near-collinearity among some of the 11 real templates (adjacent photometric bands'
depth/PSF-size maps correlate with each other) occasionally destabilizing one spatial
fold's fit once a redundant column is removed, not a loss of contamination-relevant
signal. The `injected_only` row is the sharpest illustration: reducing the design
matrix to exactly the true contaminating templates gives the best typical-case
correction (median ≈0.02-0.06) but the worst tail risk (2 of 10 realizations produce
extreme distortion), because with only 3 predictors a spatial fold whose held-out
pixels have template values outside the training folds' range can push the linear
model's out-of-fold prediction toward the singularity in the correction formula
`(1+δ_obs)/(1+p̂) - 1` at p̂→-1. We report this as a genuine numerical-stability
finding about minimal-feature-set linear correction under spatial cross-validation,
not a pipeline defect — see `research/decision_log.md` (2026-09-09 entry) for the
full per-realization breakdown.

**Out-of-distribution sky regions.** We defined an OOD region as the footprint pixels
where standardized stellar density (STARDENS) exceeds |z| > 1.5 (1,370 of 17,978
pixels, 7.6% of the footprint — a genuine "unusual observational condition," in
practice mostly Galactic-plane-proximate sky) and compared, for both linear regression
and gradient boosting, RMS error against the known true injected contamination
fraction on (a) a held-out 20% of in-distribution pixels vs. (b) the entire OOD
region, never seen during training (8 realizations).

| method | RMS error, ID holdout | RMS error, OOD region | inflation factor |
|---|---|---|---|
| linear | 0.034 | 0.059 | 1.9× |
| gradient boosting | 0.076 | 0.591 | **7.9×** |

This is the clearest genuine ML-vs-linear caveat found in this project. Gradient
boosting's error inflates roughly 4x more than linear regression's when evaluated on
sky regions unlike its training data — the expected consequence of tree-ensemble
models predicting via piecewise-constant leaf averages, which cannot extrapolate
beyond the range of template values seen during training, unlike linear regression's
smooth extrapolation. This does not contradict Section 9's headline finding (ML
matches or beats linear regression on average-case, in-distribution contamination
removal); it identifies a different, real weakness along a different axis —
generalization to unusual sky conditions, which is precisely where robust systematics
correction matters most in a real survey.

## 11. Limitations

See `research/limitations.md` for the complete, living limitations document. Key
points: this is a single-tracer, single-region, methodology-focused study, not a
precision-cosmology measurement; the injection-study mock is a simplified lognormal
field, not a full N-body/HOD forward model; injected contamination is linear-in-
templates by construction, which structurally favors the linear correction method —
notably, ML correction matched or modestly beat it anyway on this design (Section 8),
which is if anything mild evidence for ML correction rather than against it, but
whether its advantage would be larger on genuinely nonlinear real-world contamination
remains untested; the negative control fits on the noiseless true field rather than a
Poisson-noise-matched uncontaminated field, so its distortion numbers should be read
as a lower bound on what a correction method would do to real noisy data (found in a
Reviewer-2 pass, `research/reviewer2.md`); an out-of-distribution sky-region
robustness test (brief Section 23) was not run; ML hyperparameters were chosen by
convention, not a search.

We did perform external validation of the 3D baseline ξ(r) (brief Section 58): a
power-law fit over 1-30 Mpc gives r0 = 9.14 ± 0.04 Mpc/h, γ = 1.545 ± 0.009
(`results/baseline_lrg_sgc/power_law_fit.json`). Comparing to published LRG clustering
— Zehavi et al. (2005, ApJ 621, 22; arXiv:astro-ph/0411557, SDSS LRGs) report
r0 = 9.8-11.2 h⁻¹Mpc with a documented one-halo/two-halo power-law break near ~2
h⁻¹Mpc, and Zhai et al. (2017, arXiv:1607.05383, eBOSS+BOSS LRGs at the closer-matched
z~0.7) report bias b = 2.30 ± 0.03 — our r0 falls within/near the literature range,
but our γ is shallower than the literature's typical ~1.8-2.0. Restricting the fit to
r > 3 Mpc (excluding the one-halo regime) steepens γ to 1.673 ± 0.011, moving toward
the literature value but not fully closing the gap. This is reported as a genuine open
discrepancy, not resolved — plausible contributors include our smaller footprint/
random-catalog subsampling, differences in tracer selection/redshift range relative to
the comparison samples, and the compute-tier random-subsampling choices in Section 4,
none of which we have isolated as the cause.

## 12. Conclusion

We built a reproducible pipeline around real, public DESI DR1 LRG data that both (a)
measures a real 3D clustering baseline, externally validated in shape/amplitude
against published LRG clustering (Section 11), and (b) runs a controlled systematics-
correction injection study using real DESI imaging-systematics templates on the real
survey footprint. On this injection design, gradient boosting and deep random forest
correction match or modestly outperform the survey-standard linear regression at
removing injected contamination, without a corresponding signal-preservation penalty
on a negative control — negative-control distortion is approximately uniform across
every method tested (0.049-0.057). This directly contradicts a smaller pilot run of
the same pipeline, which had suggested a clear flexibility-driven overcorrection
effect; we report both results because the discrepancy is itself informative about
how much pilot-scale injection studies in this literature should be trusted before a
full-scale run.

We do not claim ML correction is generally superior to linear regression for imaging-
systematics mitigation — the injected contamination here is linear-in-templates by
construction, which structurally favors the linear method, and the negative control
is not noise-matched to the contaminated scenario. What we do claim is narrower and,
we think, still useful: under this design, at this scale, we find no evidence that
correction-model flexibility trades contamination removal for signal preservation,
which is the concern most often raised against ML systematics correction in the
literature this project is grounded in (Rezaie et al. 2020, 2021). We did find one
genuine, robust weakness of ML correction relative to linear regression: on
out-of-distribution sky regions (unusual stellar density), gradient boosting's error
inflates roughly 4x more than linear regression's (7.9x vs. 1.9x) — the expected
consequence of tree ensembles' inability to extrapolate beyond their training range.
This is a real limitation, consistent from a small pilot check through the full-scale
run, and the clearest evidence in this project that ML correction's average-case
strength on in-distribution pixels does not automatically extend to the unusual sky
regions where robust correction matters most. Robustness checks also show that
correction quality is highly sensitive to systematics-map resolution (an ~8x
degradation when templates are coarsened) and that a minimal-feature linear
correction, while best in the typical case, carries a real numerical-stability tail
risk under spatial cross-validation that a fuller feature set does not. A genuinely
nonlinear injected-contamination model and a noise-matched negative control are the
natural next steps if this project is extended.

---

*Every numerical claim above must have a corresponding row in
`research/claims_registry.md` before this document is considered final (brief Section
46: no manually typed scientific result values).*