CosmoForge: Optimal Spherical Analysis

• CosmoForge is a modular Python framework designed for optimal low‐ℓ analysis of spin-0 and spin-2 fields, integrating QML power-spectrum estimation with pixel-based likelihood evaluation.
• It smartly switches between harmonic and pixel bases to efficiently handle anisotropic noise and non-diagonal covariance matrices in large-scale cosmological data.
• The framework’s unified architecture, split into CosmoCore, QUBE, and PICSLike, ensures consistent covariance handling across various observables like CMB polarization, galaxy density, and weak-lensing shear.

CosmoForge is a public Python framework for optimal large-scale analysis of fields on the sphere, unifying Quadratic Maximum Likelihood (QML) power-spectrum estimation and pixel-based Gaussian likelihood evaluation for spin-0 and spin-2 fields with general, non-diagonal noise covariances (Galloni et al., 20 May 2026). It is designed primarily for low-$\ell$ cosmic microwave background analysis, but its formalism is stated to apply to any scalar or spin-2 observable on the sphere, including galaxy overdensity, weak-lensing shear, and stochastic gravitational-wave backgrounds. In parallel, several papers use “CosmoForge” as a broader label for integrated cosmological frameworks that connect simulation, emulation, reconstruction, or forward modelling across scales; this suggests a wider methodological family rather than a single software stack (Arnold et al., 2021).

1. Definition, scope, and scientific role

CosmoForge addresses a specific regime in which pseudo-$C_\ell$ methods are sub-optimal: large angular scales with masks and anisotropic or correlated noise. In that regime, the mode coupling induced by the sky cut and by realistic noise cannot be corrected exactly by approximate methods, while the small number of available modes makes information loss especially consequential. The framework is therefore positioned around two optimal estimators: QML for power-spectrum estimation and the pixel-space Gaussian likelihood for parameter inference (Galloni et al., 20 May 2026).

The framework is split into three installable packages: CosmoCore, which provides the infrastructure; QUBE, which provides Fisher-matrix and QML estimation; and PICSLike, which provides the pixel-space likelihood. A common interface exposes two interchangeable computation bases: a harmonic basis built on the Sherman–Morrison–Woodbury identity and a direct pixel-space basis. The code selects whichever basis is cheaper for the configuration at hand, while keeping the underlying covariance model shared between QML estimation and likelihood evaluation. This coupling is central to the framework’s stated purpose: one can estimate $C_\ell$ with QML and then perform pixel-likelihood inference using the same geometry, basis machinery, and covariance infrastructure (Galloni et al., 20 May 2026).

CosmoForge supports spin-0 and spin-2 fields in arbitrary combinations, including all auto- and cross-spectra. For spin-0, the paper cites temperature and galaxy density as examples; for spin-2, it cites CMB polarization and weak-lensing shear. It also supports general symmetric positive-definite noise covariance matrices, including anisotropic white noise, $1/f$ correlations, inter-pixel noise correlations, and inter-field correlations. This makes the framework broader than many earlier public QML implementations, which typically assumed diagonal noise or focused on only one computational formulation (Galloni et al., 20 May 2026).

2. Statistical formalism

The framework adopts the standard Gaussian pixel-space data model

$d = s + n,$

with signal covariance $S = \langle s s^\top \rangle$, noise covariance $N = \langle n n^\top \rangle$, and total covariance

$C = S + N.$

Given a theory parameter vector $\theta$, the pixel-space Gaussian likelihood is

$$\ln \mathcal{L}(\theta) = -\frac{1}{2}\left[d^\top C(\theta)^{-1} d + \ln|C(\theta)|\right] + \mathrm{const},$$

with $C_\ell$0. This likelihood is implemented in PICSLike (Galloni et al., 20 May 2026).

For QML, CosmoForge defines the covariance derivatives

$C_\ell$1

the Fisher matrix

$C_\ell$2

and the quadratic estimator components

$C_\ell$3

The resulting QML estimator is

$C_\ell$4

with covariance $C_\ell$5. The paper states that this estimator is unbiased and saturates the Cramér–Rao bound when the covariance model is correct (Galloni et al., 20 May 2026).

Cross-spectra are implemented in two ways. The first is a full joint QML treatment in which multiple maps are stacked into a single data vector with a joint covariance, yielding minimum-variance estimates of all auto- and cross-spectra. The second is a cross-QML estimator,

$C_\ell$6

with mode-mixing matrix

$C_\ell$7

The latter is unbiased and noise-bias-free when the two datasets have independent noise, but the paper explicitly notes that it is not minimum-variance because it does not exploit the off-diagonal signal correlations in the full joint covariance (Galloni et al., 20 May 2026).

The spin structure is handled through a field-collection abstraction. In harmonic space, the signal covariance becomes diagonal in $C_\ell$8 but retains block structure in field indices. For a single spin-2 field, the paper writes

$C_\ell$9

while for a joint $C_\ell$0 analysis it uses the corresponding $C_\ell$1 block. This makes the framework applicable to $C_\ell$2, $C_\ell$3, $C_\ell$4, $C_\ell$5, $C_\ell$6, and $C_\ell$7, as well as to arbitrary multi-field combinations (Galloni et al., 20 May 2026).

3. Computational architecture and algorithmic design

The framework’s central algorithmic distinction is between two computation bases. In the harmonic basis, the signal covariance is factorized as

$C_\ell$8

so that the total covariance is

$C_\ell$9

Applying the Sherman–Morrison–Woodbury identity with

$1/f$0

gives

$1/f$1

The matrix determinant lemma then yields

$1/f$2

In the pixel basis, by contrast, $1/f$3 is built directly in pixel space, $1/f$4 is formed explicitly, and Cholesky factorization supplies both $1/f$5 and $1/f$6 (Galloni et al., 20 May 2026).

A prominent implementation detail is numerical stabilization of the projected inverse

$1/f$7

The paper notes that the naive SMW form,

$1/f$8

can suffer catastrophic cancellation in high signal-to-noise regimes. CosmoForge therefore uses the algebraically equivalent stable form

$1/f$9

implemented through linear solves rather than subtraction of nearly equal matrices. The stated consequence is improved numerical stability and positive-definite Fisher matrices even in cosmic-variance-dominated regimes (Galloni et al., 20 May 2026).

The Fisher stage is further accelerated by exploiting the sparsity of the harmonic derivative matrices $d = s + n,$0, where $d = s + n,$1. Using cyclicity of the trace,

$d = s + n,$2

and evaluating only the non-zero entries of $d = s + n,$3, the framework reduces the Fisher cost to $d = s + n,$4 for arbitrary noise covariances. The paper identifies this as an exact algorithmic optimization, not an approximation (Galloni et al., 20 May 2026).

Several additional optimizations are built in. CosmoForge supports automatic field block-diagonalization when field groups decouple, reducing cubic costs in the SMW kernel. It uses Numba JIT compilation for hot kernels, including spin-0 Legendre polynomials, spin-2 Wigner $d = s + n,$5-matrix recurrences, and pixel-space signal-matrix construction. It uses MPI parallelization in the QML Fisher stage, in spectrum estimation, and in the pixel-likelihood scan. It also implements multipole switching, following Planck low-$d = s + n,$6, so that only a chosen multipole interval is varied while the rest of the spectrum is absorbed into an effective noise term. This reduces the effective number of modes and hence the cost of the SMW inversion (Galloni et al., 20 May 2026).

4. Outputs, validation, and performance

CosmoForge provides native multipole binning and three output normalizations through a single code path. The paper emphasizes that native binning operates directly on the sufficient statistics, rather than averaging already-estimated unbinned $d = s + n,$7, and is therefore preferred when the Fisher matrix has substantial off-diagonal structure (Galloni et al., 20 May 2026).

Output mode	Definition	Interpretation
Deconvolved	$d = s + n,$8	Estimates of the true $d = s + n,$9
Decorrelated	$S = \langle s s^\top \rangle$0	Unit covariance, mixed expectation value
Convolved	Return $S = \langle s s^\top \rangle$1 with window $S = \langle s s^\top \rangle$2	Theory is compared after convolution

Validation is performed against the Planck low-$S = \langle s s^\top \rangle$3 Fortran reference implementation. In a $S = \langle s s^\top \rangle$4-mode-only configuration with $S = \langle s s^\top \rangle$5 and $S = \langle s s^\top \rangle$6, the paper reports bit-exact agreement for the noise covariance $S = \langle s s^\top \rangle$7, and matrix-level agreement for $S = \langle s s^\top \rangle$8, $S = \langle s s^\top \rangle$9, $N = \langle n n^\top \rangle$0, $N = \langle n n^\top \rangle$1, $N = \langle n n^\top \rangle$2, $N = \langle n n^\top \rangle$3, and $N = \langle n n^\top \rangle$4. Typical relative differences are reported as $N = \langle n n^\top \rangle$5 for $N = \langle n n^\top \rangle$6, $N = \langle n n^\top \rangle$7 for Fisher and QML products, and $N = \langle n n^\top \rangle$8 for the likelihood, with the differences attributed to floating-point accumulation order rather than systematic bias. Harmonic and pixel bases are stated to agree internally to $N = \langle n n^\top \rangle$9 for the Fisher matrix in both spin-0 and multi-field cases (Galloni et al., 20 May 2026).

Monte Carlo tests are carried out at $C = S + N.$0 for $C = S + N.$1 over $C = S + N.$2 using 10 000 Gaussian polarization realizations with realistic beam and mask. The recovered mean $C = S + N.$3 spectrum matches the input, residuals are much smaller than the standard error on the mean, and empirical standard deviations match $C = S + N.$4. The paper stresses that this is a stringent test because $C = S + N.$5 is subdominant to $C = S + N.$6, so any $C = S + N.$7 leakage would appear as a positive bias (Galloni et al., 20 May 2026).

The framework is not intended for arbitrarily high resolutions or full-sky high-$C = S + N.$8 analysis without binning. The paper states that practical QML analyses are targeted to low and intermediate multipoles, sky fractions from $C = S + N.$9 to $\theta$0, and $\theta$1, depending on sky fraction and basis choice. Benchmark figures illustrate the intended regime: for $\theta$2, pixel basis is preferred, and a $\theta$3 analysis at $\theta$4 with about 15 600 pixels on a 48-core node completes full Fisher plus 10 QML simulations in $\theta$5 minutes; for $\theta$6 and $\theta$7 at $\theta$8, Fisher plus 10 simulations takes $\theta$9 s, with QML cost per simulation of about 4 s. The sparse Fisher traces are reported to reduce Fisher-stage runtime by factors $$\ln \mathcal{L}(\theta) = -\frac{1}{2}\left[d^\top C(\theta)^{-1} d + \ln|C(\theta)|\right] + \mathrm{const},$$0 relative to naive dense harmonic implementations (Galloni et al., 20 May 2026).

A comparison with pseudo-$$\ln \mathcal{L}(\theta) = -\frac{1}{2}\left[d^\top C(\theta)^{-1} d + \ln|C(\theta)|\right] + \mathrm{const},$$1 analysis using NaMaster quantifies the gain from optimal methods. For $$\ln \mathcal{L}(\theta) = -\frac{1}{2}\left[d^\top C(\theta)^{-1} d + \ln|C(\theta)|\right] + \mathrm{const},$$2, $$\ln \mathcal{L}(\theta) = -\frac{1}{2}\left[d^\top C(\theta)^{-1} d + \ln|C(\theta)|\right] + \mathrm{const},$$3, and $$\ln \mathcal{L}(\theta) = -\frac{1}{2}\left[d^\top C(\theta)^{-1} d + \ln|C(\theta)|\right] + \mathrm{const},$$4, QML errors are reported to be 35–60% smaller than pseudo-$$\ln \mathcal{L}(\theta) = -\frac{1}{2}\left[d^\top C(\theta)^{-1} d + \ln|C(\theta)|\right] + \mathrm{const},$$5 at $$\ln \mathcal{L}(\theta) = -\frac{1}{2}\left[d^\top C(\theta)^{-1} d + \ln|C(\theta)|\right] + \mathrm{const},$$6. Even at $$\ln \mathcal{L}(\theta) = -\frac{1}{2}\left[d^\top C(\theta)^{-1} d + \ln|C(\theta)|\right] + \mathrm{const},$$7, pseudo-$$\ln \mathcal{L}(\theta) = -\frac{1}{2}\left[d^\top C(\theta)^{-1} d + \ln|C(\theta)|\right] + \mathrm{const},$$8 error bars remain up to $$\ln \mathcal{L}(\theta) = -\frac{1}{2}\left[d^\top C(\theta)^{-1} d + \ln|C(\theta)|\right] + \mathrm{const},$$9 larger at the lowest multipoles, and QML decorrelated spectra have nearly diagonal correlation matrices while pseudo-$C_\ell$00 retains significant off-diagonal correlations (Galloni et al., 20 May 2026).

5. Related “CosmoForge-type” frameworks

Several contemporaneous papers explicitly describe their systems as prototypes or templates for a “CosmoForge”-type framework, but in domains distinct from low-$C_\ell$01 spherical likelihood analysis. These usages do not denote the same software as CosmoForge proper; rather, they indicate a broader architectural idea: cosmological pipelines that connect large-scale initial conditions or simulations to optimized inference or realistic observables.

Paper	Domain	Stated “CosmoForge”-type role
"FORGE -- the f(R) gravity cosmic emulator project I: Introduction and matter power spectrum emulator" (Arnold et al., 2021)	Modified-gravity emulation	A prototype for a “CosmoForge”-type framework turning a carefully designed suite of simulations into a Gaussian-process emulator for nonlinear $C_\ell$02
"FORECAST: a flexible software to forward model cosmological hydrodynamical simulations mimicking real observations" (Fortuni et al., 2023)	Survey forward modelling	Organized as ingredients one would want in a “CosmoForge”-type framework linking hydrodynamical simulations to realistic images and catalogs
"Introducing cosmosTNG: simulating galaxy formation with constrained realizations of the COSMOS field" (Byrohl et al., 2024)	Constrained hydrodynamical simulation	Presented as a complete end-to-end prototype for a flexible cosmological simulation or modeling framework that ingests survey data, reconstructs initial conditions, runs hydrodynamical simulations, and produces synthetic survey products
"Fast Lightcones for Combined Cosmological Probes" (Sgier et al., 2020)	Multi-probe lightcones and covariance	Treated as an architectural template for a CosmoForge-like framework for self-consistent full-sky maps and non-Gaussian multi-probe covariance
"FORGE'd in the Early Universe: The Effect of Protostellar Outflows on Pop III Accretion" (Meziani et al., 3 Feb 2026)	Cosmological RMHD star formation	Described as a prototype of a “CosmoForge” framework embedding STARFORGE-like star-formation physics into FIRE-like cosmological zoom-ins

These examples span several methodological layers. The modified-gravity FORGE project uses 200 dark-matter-only simulations over a Latin hypercube in $C_\ell$03 and a Gaussian-process emulator for the nonlinear matter power spectrum, reaching accuracy better than $C_\ell$04 for the majority of nodes up to $C_\ell$05 (Arnold et al., 2021). FORECAST performs particle-based forward modelling of hydrodynamical simulations into realistic synthetic images and catalogs, including light-cone construction, SSP-based flux computation, IGM attenuation, dust attenuation, PSF convolution, and survey-like noise (Fortuni et al., 2023). cosmosTNG couples observationally constrained initial conditions from CLAMATO and zCOSMOS, reconstructed with TARDIS, to IllustrisTNG hydrodynamical evolution in eight realizations of a COSMOS subvolume at $C_\ell$06 (Byrohl et al., 2024). UFalcon produces self-consistent full-sky maps for weak lensing shear, galaxy overdensity including RSD, CMB lensing, and ISW temperature anisotropies from the same underlying density field, with post-processing in about one walltime-hour corresponding to about 100 CPU-hours (Sgier et al., 2020). FORGE’d in FIRE combines FIRE-3 cosmological zoom-ins, FORGE hyper-refinement, and STARFORGE-grade RMHD star-formation physics to follow a single metal-free Pop III protostar from Mpc scales to sub-100 au scales (Meziani et al., 3 Feb 2026).

Taken together, these papers suggest a broader usage in which “CosmoForge” denotes an integrated framework joining cosmological structure, physical modelling, and observational or inferential products. A plausible implication is that the low-$C_\ell$07 CosmoForge code (Galloni et al., 20 May 2026) belongs to the same methodological trend: replacing isolated analysis tools with modular, end-to-end systems whose internal covariance and numerical assumptions remain consistent across stages.

6. Limitations, extensibility, and prospective significance

CosmoForge’s intended regime is explicit. Its computational complexity is dominated by $C_\ell$08 in pixel basis and by $C_\ell$09 plus $C_\ell$10 in harmonic basis, so it is aimed at low and intermediate multipoles rather than full high-$C_\ell$11 CMB analysis. The paper therefore positions pseudo-$C_\ell$12 methods such as NaMaster as the practical workhorses for high-$C_\ell$13 applications, while CosmoForge targets the large-scale regime where exact optimality matters most (Galloni et al., 20 May 2026).

Within that regime, the framework is designed to be extensible. It is publicly available on GitHub and on PyPI as cosmocore, qube-qml, picslike, and the meta-package cosmoforge. Because the internal abstractions are built around generic spin-0 and spin-2 fields, the same machinery can be applied to galaxy clustering, weak-lensing shear, 21 cm intensity mapping, stochastic gravitational-wave background anisotropies, and cross-correlations among such fields, provided the inputs are maps on the sphere, an angular power-spectrum description of the signal, and a noise covariance model (Galloni et al., 20 May 2026).

The broader “CosmoForge-type” literature points in complementary directions for future extension. FORGE indicates how a simulation suite can be turned into an emulator for fast inference in extended-gravity parameter spaces (Arnold et al., 2021). FORECAST shows how hydro outputs can be propagated into realistic synthetic survey images and catalogs for end-to-end validation of photometric pipelines (Fortuni et al., 2023). cosmosTNG demonstrates an observations-to-initial-conditions-to-hydrodynamics workflow for specific survey fields (Byrohl et al., 2024). UFalcon provides a template for rapid lightcone construction and simulation-based non-Gaussian covariance estimation for combined probes (Sgier et al., 2020). FORGE’d in FIRE illustrates how cosmological initial conditions can be connected to AU-scale RMHD feedback in first-star formation (Meziani et al., 3 Feb 2026). These adjacent developments do not alter the software identity of CosmoForge proper, but they clarify the conceptual space into which it fits.

In that sense, CosmoForge occupies a sharply defined place: it is an optimal, modular, validated framework for low-$C_\ell$14 spherical-data analysis, centered on exact covariance handling, interchangeable harmonic and pixel formulations, and the unification of QML estimation with pixel-based likelihood inference (Galloni et al., 20 May 2026). The wider literature suggests that the same design logic—shared covariance infrastructure, modular computation backends, and consistency across inference stages—is becoming characteristic of cosmological analysis frameworks more generally.