Power Spectrum Estimation

The Spectra class implements Quadratic Maximum Likelihood (QML) estimation for power spectrum recovery from partial-sky observations of spin-0 and spin-2 fields. This method provides unbiased, minimum-variance power spectrum estimates while properly accounting for noise, masking, and correlations. It applies to any Gaussian field on the sphere, such as CMB temperature and polarization, galaxy density, or 21 cm intensity.

Mathematical Background

QML Theory

QML power spectrum estimation maximizes the likelihood function for Gaussian random fields. For each multipole \(\ell\), the QML estimator is:

\[\hat{C}_\ell = \frac{\mathbf{x}^T \mathbf{E}_\ell \mathbf{x}}{\text{Tr}[\mathbf{E}_\ell \mathbf{S}]}\]

where: - \(\mathbf{x}\) is the observed data vector - \(\mathbf{E}_\ell = \mathbf{C}^{-1} \frac{\partial \mathbf{C}}{\partial C_\ell} \mathbf{C}^{-1}\) is the estimator matrix - \(\mathbf{S}\) is the total covariance matrix - \(\mathbf{C} = \mathbf{S} + \mathbf{N}\) includes signal and noise

This estimator is unbiased: \(\langle \hat{C}_\ell \rangle = C_\ell\) and achieves the minimum variance bound given by the Fisher information.

Bias Correction

The normalized QML estimator includes bias correction for finite-size effects:

\[q_\ell = \frac{\mathbf{x}^T \mathbf{E}_\ell \mathbf{x}}{\text{Tr}[\mathbf{E}_\ell \mathbf{M}]}\]

where \(\mathbf{M}\) is the appropriate normalization matrix that accounts for the specific covariance structure.

Class Documentation

class qube.spectra.Spectra(params_file=None, fisher=None, compression=None, **kwargs)[source]

Bases: Core, MPISharedMemoryMixin

Quadratic Maximum Likelihood (QML) power spectrum estimator.

This class implements the QML method for optimal power spectrum estimation from observations of spin-0 and spin-2 fields on the sphere (e.g. CMB temperature and polarization, cosmic shear). The QML estimator provides unbiased, minimum-variance estimates of angular power spectra C_l while properly accounting for instrumental noise, sky cuts, and pixel correlations.

The QML estimator computes power spectrum amplitudes as:

q̂_l = (1/2) * x^T * E_l * x

where E_l = C^(-1) * ∂C/∂q_l * C^(-1) is the quadratic estimator matrix, x is the data vector, and C is the total covariance matrix.

Parameters:

  • params_file (str, optional) – Path to YAML parameter file containing analysis configuration.

  • fisher (Fisher, optional) – Pre-computed Fisher matrix instance. If provided, reuses computed components (covariance matrices, geometry, etc.) for efficiency.

  • **kwargs (dict) – Additional keyword arguments passed to the Core parent class.

comm

MPI communicator for parallel computation.

Type:

MPI.Comm

rank

MPI process rank (0 for master process).

Type:

int

size

Total number of MPI processes.

Type:

int

fisher_instance

Fisher matrix computation instance, either provided or computed.

Type:

Fisher

maps1, maps2

Input map data for primary and secondary fields (cross-correlation).

Type:

numpy.ndarray

qml_results

QML power spectrum estimates for all simulations and multipoles.

Type:

numpy.ndarray

qml_noise_bias

Noise bias estimates for auto-correlation spectra.

Type:

numpy.ndarray

invfisher

Inverse of the beam-smoothed Fisher matrix.

Type:

numpy.ndarray

inv_cov1, inv_cov2

Inverted covariance matrices for primary and secondary datasets.

Type:

numpy.ndarray

Examples

Basic QML power spectrum estimation:

>>> from cosmoforge.qube import Spectra
>>> spectra = Spectra("config/qml_analysis.yaml")
>>> spectra.run()
>>> power_spectra = spectra.get_power_spectra()
>>> noise_bias = spectra.get_noise_bias()

Using pre-computed Fisher matrix:

>>> from cosmoforge.qube import Fisher, Spectra
>>> fisher = Fisher("config/fisher_config.yaml")
>>> fisher.run()
>>> spectra = Spectra("config/qml_config.yaml", fisher=fisher)
>>> spectra.run()  # Reuses Fisher components for efficiency

MPI parallel execution:

>>> # Run with: mpirun -n 8 python qml_analysis.py
>>> spectra = Spectra("config/high_res_config.yaml")
>>> spectra.run()  # Distributes computation across processes

Notes

The QML method offers several advantages over pseudo-C_l estimators: 1. Optimal statistical properties (minimum variance, unbiased) 2. Exact treatment of sky cuts and inhomogeneous noise 3. Proper error propagation through Fisher matrix 4. Natural handling of mode coupling

However, it is computationally expensive, scaling as O(N_pix^3) due to matrix inversions and O(N_ell * N_pix^2) for quadratic form evaluations.

For auto-correlation analyses, noise bias is computed and can be subtracted: noise_bias_l = (1/2) * Tr[N * E_l] where N is the noise covariance matrix.

Cross-correlation analyses are naturally noise-bias free when using independent realizations of noise between maps.

The implementation supports both temperature-only and polarization analyses with proper treatment of spin-2 fields and E/B mode decomposition.

See also

Fisher

Fisher information matrix computation

cosmocore.Core

Base class providing fundamental analysis infrastructure

__init__(params_file=None, fisher=None, compression=None, **kwargs)[source]

Initialize QML power spectrum estimation class.

Parameters:

  • params_file (str, optional) – Path to YAML configuration file.

  • fisher (Fisher, optional) – Pre-computed Fisher instance. If provided, reuses computed components (covariance matrices, geometry, field collections) for efficiency.

  • compression (dict, optional) – Computation basis configuration (method, epsilon, basis, mode_fraction).

  • **kwargs (dict) – Additional arguments passed to Core.

Raises:

  • TypeError – If fisher is not a Fisher instance.

  • ValueError – If fisher doesn’t contain a valid Fisher matrix.

property lmax_signal: int

Maximum multipole for signal/derivative matrix computation.

This defaults to 4*nside to match the Fortran reference implementation. The derivative matrices are computed up to this lmax, while the output power spectra use params.lmax.

Returns:

Maximum multipole for signal covariance and derivative computation.

Return type:

int

setup_maps()[source]

Read observational map data from FITS files.

Loads maps with proper pixel selection, field extraction, and calibration. For cross-correlation (do_cross=True), loads both primary and secondary maps. Output shape: (n_active_pixels, n_simulations).

Raises:
setup_fisher_inversion()[source]

Prepare inverted Fisher matrix for QML estimation.

The Fisher matrix is already beam-smoothed (beam window functions absorbed into derivatives). This method computes F⁻¹, F⁻¹/², and writes results to output files.

Raises:

  • ValueError – If Fisher matrix is not available or singular.

  • LinAlgError – If Fisher matrix inversion fails.

setup_qml_computation()[source]

Initialize arrays and variables for QML power spectrum computation.

This method allocates memory for QML estimation results and prepares the computational infrastructure for the main QML calculation phase. Arrays are sized based on the number of multipoles, spectra, and Monte Carlo simulations.

Notes

Initializes the following arrays:

  • qml_results: Storage for QML power spectrum estimates with shape (n_simulations, n_total_parameters) where n_total_parameters = n_spectra * (lmax - 1)

  • qml_noise_bias: Storage for noise bias estimates (auto-correlation only) with shape (n_total_parameters,)

The qml_results array stores the raw quadratic estimates: q̂_l^(sim) = (1/2) * x^(sim)T * E_l * x^(sim)

for each simulation and multipole. These raw estimates will be combined using the inverse Fisher matrix to produce the final optimally-weighted power spectrum estimates.

For auto-correlation analyses, noise bias terms are computed as: bias_l = (1/2) * Tr[N * E_l] where N is the noise covariance matrix and E_l is the quadratic estimator matrix for multipole l.

Memory allocation scales as O(n_sims * n_ell * n_spec) for results and O(n_ell * n_spec) for bias terms, which is typically much smaller than the O(n_pix^2) covariance matrix storage.

Examples

This method is called automatically during computation:

>>> spectra = Spectra("config.yaml")
>>> spectra.setup_qml_computation()
>>> print(f"QML results shape: {spectra.qml_results.shape}")
>>> print(f"Total parameters: {spectra.qml_results.shape[1]}")
compute_qml_spectra()[source]

Execute parallel QML power spectrum computation.

Distributes multipole computations across MPI processes (round-robin). For each multipole: computes signal derivative, E-operator, and quadratic estimates for all simulations. Selects compressed or traditional method based on basis_manager availability.

compute()[source]

Execute QML power spectrum computation.

Calls setup_qml_computation() to initialize arrays, then compute_qml_spectra() for the main parallel computation.

run()[source]

Execute the complete QML power spectrum analysis pipeline.

Pipeline phases: 1. Master setup: fields, geometry, covariance, spectra, beams

(reuses Fisher components if provided)

  1. QML setup: maps, Fisher inversion

  2. MPI broadcast of shared data to workers

  3. Parallel QML computation

  4. MPI synchronization

Raises:
get_power_spectra(mode='deconvolved', *, as_dict=False)[source]

Retrieve power spectrum estimates in specified normalization mode.

This method returns power spectrum estimates with three different normalization prescriptions, allowing users to choose the output format best suited to their analysis needs.

Parameters:

  • mode (str, optional) –

    Normalization mode for output (default: “deconvolved”):

    • ”deconvolved”: F⁻¹y - estimates of true C_ℓ with correlated errors. This is the standard QML output that attempts to recover the true underlying spectrum by inverting the window function.

    • ”decorrelated”: F⁻¹/²y - uncorrelated bandpower estimates with identity covariance matrix. Useful when independent error bars are needed.

    • ”convolved”: Raw y estimates with window matrix W for theory comparison. Instead of deconvolving the window, compare with window-convolved theory: <y> = W @ C_true.

  • as_dict (bool, optional) – If False (default, back-compat), return the flat numpy array described below. If True, return a dict keyed by physical labels ("TT", "EE", "TE", …) whose values have shape (n_simulations, n_bins). In "convolved" mode, the dict applies to the raw estimates y only; the window matrix W and convolve_theory callable are unchanged (use Fisher.get_bandpower_slices() on self.fisher_instance to navigate W).

Return type:

ndarray | tuple[ndarray, ndarray, Callable] | dict | None

Returns:

  • numpy.ndarray or tuple or dict or None – For “deconvolved” or “decorrelated” modes with as_dict=False:

    Array of shape (n_simulations, n_parameters) where n_parameters = n_spectra * n_bins.

    For “convolved” mode with as_dict=False:

    Tuple of (y, W, convolve_theory_func) where:

    • y: Raw estimates of shape (n_simulations, n_parameters)

    • W: Window matrix of shape (n_parameters, n_parameters)

    • convolve_theory_func: Callable that applies W @ theory

    With as_dict=True: a dict for the per-spectrum arrays (and, in convolved mode, a 3-tuple (y_dict, W, convolve)).

    Returns None for worker processes (rank != 0) or if computation has not completed.

  • **Layout (flat-array shape)**

  • The columns of the returned array are organised as n_spectra

  • contiguous blocks of n_bins columns each. Spectrum order

  • follows Fisher.spectra_list on

  • self.fisher_instance — auto-spectra first (per component,

  • in the order each component emits them), then cross-spectra in

  • (i, j) order with i < j. For TQU (spins=[0, 2],

  • labels=["T", "E", "B"]) under SymmetryMode.SYMMETRIC

  • with nbins=10, n_parameters = 60 and the columns are:: – cols 0..9 : TT cols 10..19 : EE cols 20..29 : BB cols 30..39 : EB cols 40..49 : TE cols 50..59 : TB

  • Use Fisher.get_bandpower_slices() (on

  • self.fisher_instance) or pass as_dict=True to look up

  • per-spectrum slices without computing these offsets by hand.

Raises:

ValueError – If mode is not one of “deconvolved”, “decorrelated”, “convolved”.

Notes

Mode Comparison:

Mode | Formula | Covariance | Use Case |

|------|———|------------|———-| | deconvolved | F⁻¹y | F⁻¹ (correlated) | Standard analysis | | decorrelated | F⁻¹/²y | I (identity) | Independent errors | | convolved | y | F | Theory comparison |

Deconvolved Mode (default): The standard QML output that inverts the window function to recover estimates of the true C_ℓ. Errors are correlated between multipoles.

Decorrelated Mode: Produces bandpower estimates with uncorrelated errors (unit variance). Useful for plotting error bars or when independent measurements are required. Note that information content is preserved but redistributed.

Convolved Mode: Returns raw QML estimates without deconvolution, along with the window matrix for comparing with convolved theoretical spectra. This avoids numerical issues from inverting poorly-conditioned window matrices.

Output Convention: When output_convention is set to "Dl" in the configuration, all returned spectra are converted from C_ℓ to D_ℓ = ℓ(ℓ+1)/(2π) C_ℓ. In convolved mode, the window matrix is also transformed so that the returned convolve_theory_func expects D_ℓ input.

Examples

Default (flat array) — back-compat:

>>> spectra = Spectra("config.yaml")
>>> spectra.run()
>>> cl_deconv = spectra.get_power_spectra()        # shape (n_sims, 60)
>>> ee_slice = cl_deconv[:, 10:20]                 # bins 0..9 of EE

Label-keyed dict (recommended for human-readable code):

>>> cl = spectra.get_power_spectra(as_dict=True)
>>> cl["EE"].shape                                 # (n_sims, n_bins)
>>> cl["TE"]                                       # cross-spectrum

Decorrelated bandpowers:

>>> cl_decorr = spectra.get_power_spectra(mode="decorrelated")
>>> # Errors are all 1.0 by construction

Convolved mode for theory comparison:

>>> y, W, convolve = spectra.get_power_spectra(mode="convolved")
>>> cl_theory_convolved = convolve(cl_theory)
>>> # Compare: y should match cl_theory_convolved

>>> # With as_dict=True, the y component becomes a dict; W and
>>> # the convolve callable are unchanged.
>>> y_dict, W, convolve = spectra.get_power_spectra(
...     mode="convolved", as_dict=True
... )

See also

get_covariance

Get covariance matrix for specified mode

get_error_bars

Get 1-sigma error bars for specified mode

convolve_theory

Apply window matrix to theoretical spectrum

Fisher.get_bandpower_slices

Slice mapping for navigating the flat layout.

get_effective_ells()[source]

Return effective multipole for each bin.

For unbinned (delta_ell=1), returns integer ells from 2 to lmax. For binned, returns bin midpoints.

Returns:

Effective multipoles, shape (nbins,). None for workers.

Return type:

np.ndarray or None

get_noise_bias()[source]

Retrieve noise bias estimates for auto-correlation spectra.

Returns:

Shape (n_params,). Returns None for workers, cross-correlation analyses (do_cross=True), or if computation incomplete.

Return type:

np.ndarray or None

Notes

Noise bias: (1/2) * Tr[N * E_l]. Cross-correlations are bias-free.

convolve_theory_for_inference(cl_theory)[source]

Apply the QML bandpower window to a per-ℓ theory spectrum.

Returns the expected binned bandpower for the given theory, in the same convention as get_power_spectra(mode="deconvolved"). Use this in a Gaussian likelihood for parameter inference:

mu_b(θ) = convolve_theory_for_inference(C_ℓ(θ)) -2 ln L = (d - mu)ᵀ F_b (d - mu)

Parameters:

cl_theory (np.ndarray) –

Per-ℓ theory C_ℓ. Two formats accepted: - shape (n_ell,) for ℓ=lmin..lmax (length lmax-lmin+1) - shape (lmax+1,) starting at ℓ=0 (entries below lmin

are ignored)

Must be unbeamed physical C_ℓ — beam² is already absorbed into the window.

Returns:

cl_binned – Expected binned bandpower for the given theory. If the output convention is “Dl”, the result is converted to D_ℓ using the bin effective ells. Returns None on worker processes (rank != 0).

Return type:

np.ndarray of shape (n_bins,), or None

Notes

This delegates to Fisher.get_bandpower_window_function() on the underlying Fisher instance. The first call triggers a per-ℓ Fisher computation (cached for subsequent calls).

Scope: single-spectrum only. For multi-spectrum likelihoods, use get_power_spectra() (mode="convolved") — the returned convolve_theory_func accepts a label-keyed dict and handles the full n_spectra * n_bins window matrix.

See Fisher.get_bandpower_window_function() for details on the buffer approach and the multi-spectrum limitation.

Examples

>>> spectra = Spectra("config.yaml", fisher=fisher,
...                   compression={"method": "harmonic"})
>>> spectra.set_binning(bins)
>>> spectra.run()
>>> cl_th = compute_camb_cl(theta)        # ell=0..lmax
>>> mu_b = spectra.convolve_theory_for_inference(cl_th)
>>> data = spectra.get_power_spectra(mode="deconvolved").mean(0)
>>> F_b = spectra.fisher_instance.get_fisher_matrix()
>>> neg2lnL = (data - mu_b) @ F_b @ (data - mu_b)

See also

Fisher.get_bandpower_window_function

Underlying window matrix.

get_power_spectra

Get the data bandpower (deconvolved mode).

get_covariance

Get the covariance matrix (= F_b⁻¹).

get_covariance(mode='deconvolved')[source]

Get covariance matrix for power spectrum estimates in specified mode.

Returns the covariance matrix appropriate for the normalization mode, which quantifies the statistical uncertainties and correlations between different multipole estimates.

Parameters:

mode (str, optional) –

Normalization mode (default: “deconvolved”):

  • ”deconvolved”: Returns F⁻¹, the inverse Fisher matrix with correlated errors between multipoles.

  • ”decorrelated”: Returns identity matrix I, as decorrelated estimates have unit variance by construction.

  • ”convolved”: Returns F, the normalized Fisher matrix, which is the covariance of raw y estimates.

Returns:

Covariance matrix of shape (n_params, n_params). Returns None for worker processes (rank != 0) or if matrices are not available.

Return type:

numpy.ndarray or None

Raises:

ValueError – If mode is not one of “deconvolved”, “decorrelated”, “convolved”.

Notes

Mode-specific covariances:

Mode | Covariance | Interpretation |

|------|————|----------------| | deconvolved | F⁻¹ | Correlated errors on C_ℓ estimates | | decorrelated | I | Unit variance (uncorrelated) | | convolved | F | Covariance of raw y estimates |

The covariance matrix is essential for: - Computing chi-square statistics - Parameter estimation from power spectra - Constructing likelihood functions - Understanding multipole correlations

Examples

Get covariance for chi-square computation:

>>> cov = spectra.get_covariance(mode="deconvolved")
>>> cl_data = spectra.get_power_spectra(mode="deconvolved")
>>> residuals = np.mean(cl_data, axis=0) - cl_theory
>>> chi2 = residuals @ np.linalg.inv(cov) @ residuals

Decorrelated mode has identity covariance:

>>> cov_decorr = spectra.get_covariance(mode="decorrelated")
>>> # cov_decorr is identity matrix

See also

get_power_spectra

Get power spectrum estimates

get_error_bars

Get 1-sigma error bars (sqrt of diagonal)

get_error_bars(mode='deconvolved')[source]

Get 1-sigma error bars for power spectrum estimates.

Returns the marginal standard deviations for each multipole estimate, computed from the diagonal of the covariance matrix.

Parameters:

mode (str, optional) –

Normalization mode (default: “deconvolved”):

  • ”deconvolved”: sqrt(diag(F⁻¹)) - standard QML errors

  • ”decorrelated”: all ones (unit variance by construction)

  • ”convolved”: sqrt(diag(F)) - errors on raw y estimates

Returns:

1D array of error bars with shape (n_params,). Returns None for worker processes (rank != 0) or if covariance not available.

Return type:

numpy.ndarray or None

Raises:

ValueError – If mode is not one of “deconvolved”, “decorrelated”, “convolved”.

Notes

The error bars are computed as: σ_i = sqrt(Cov_ii)

where Cov is the mode-appropriate covariance matrix. These represent marginal uncertainties on individual estimates, marginalized over all other multipoles.

Important: For deconvolved and convolved modes, errors are correlated between multipoles. The full covariance matrix (from get_covariance) should be used for proper statistical analysis.

For decorrelated mode, error bars are all 1.0 by construction, making them suitable for simple error bar plots.

Examples

Get error bars for plotting:

>>> errors = spectra.get_error_bars(mode="deconvolved")
>>> cl = np.mean(spectra.get_power_spectra(), axis=0)
>>> plt.errorbar(ell, cl, yerr=errors)

Decorrelated mode has unit errors:

>>> errors_decorr = spectra.get_error_bars(mode="decorrelated")
>>> # All elements are 1.0

See also

get_covariance

Get full covariance matrix

get_power_spectra

Get power spectrum estimates

convolve_theory(cl_theory)[source]

Apply window matrix to theoretical power spectrum.

Convolves a theoretical power spectrum with the QML window matrix, producing the expected value of the raw QML estimates y for that theory: <y> = W @ C_theory.

Parameters:

cl_theory (numpy.ndarray) – Theoretical power spectrum values. Should be a 1D array with shape (n_params,) matching the QML output dimensions. When output_convention="Dl", this should be D_ℓ values.

Returns:

Window-convolved theoretical spectrum with shape (n_params,), in the same convention as the input. Returns None if window matrix is not available.

Return type:

numpy.ndarray or None

Notes

This method is useful for comparing QML estimates with theoretical predictions when using the “convolved” normalization mode. Instead of deconvolving the window function from the data (which can be numerically unstable), the window is applied to the theory.

The relationship is: <y> = W @ C_true

where y are the raw QML estimates and W is the window matrix. When output_convention="Dl", the window matrix is internally transformed so the input and output are both in D_ℓ convention.

Examples

Compare convolved estimates with theory:

>>> y, W, _ = spectra.get_power_spectra(mode="convolved")
>>> cl_theory_convolved = spectra.convolve_theory(cl_theory)
>>> residuals = np.mean(y, axis=0) - cl_theory_convolved

See also

get_power_spectra

Get power spectra (convolved mode returns window)

write_power_spectra(mode='deconvolved', filename=None, include_errors=True)[source]

Write power spectra to file in specified normalization mode.

Outputs power spectrum estimates and optionally error bars to text files. For convolved mode, also writes the window matrix.

Parameters:

  • mode (str, optional) – Normalization mode (default: “deconvolved”): “deconvolved”, “decorrelated”, or “convolved”.

  • filename (str or None, optional) – Output filename. If None, auto-generates based on mode: “{outclfile_base}_{mode}.{ext}”

  • include_errors (bool, optional) – If True, writes error bars to a separate file (default: True). For convolved mode, this parameter is ignored.

Return type:

None

Notes

File outputs by mode:

  • deconvolved/decorrelated: Writes mean spectrum across simulations and optionally error bars from diagonal of covariance.

  • convolved: Writes mean raw y estimates and the window matrix W to separate files (filename and filename_window).

Files are only written on the master process (rank 0).

Examples

Write all modes:

>>> spectra.write_power_spectra(mode="deconvolved")
>>> spectra.write_power_spectra(mode="decorrelated")
>>> spectra.write_power_spectra(mode="convolved")

Custom filename:

>>> spectra.write_power_spectra(
...     mode="deconvolved",
...     filename="my_results.dat"
... )

See also

get_power_spectra

Get power spectrum estimates

get_error_bars

Get 1-sigma error bars

Key Methods

Initialization

Spectra.__init__(params_file=None, fisher=None, compression=None, **kwargs)[source]

Initialize QML power spectrum estimation class.

Parameters:

  • params_file (str, optional) – Path to YAML configuration file.

  • fisher (Fisher, optional) – Pre-computed Fisher instance. If provided, reuses computed components (covariance matrices, geometry, field collections) for efficiency.

  • compression (dict, optional) – Computation basis configuration (method, epsilon, basis, mode_fraction).

  • **kwargs (dict) – Additional arguments passed to Core.

Raises:

  • TypeError – If fisher is not a Fisher instance.

  • ValueError – If fisher doesn’t contain a valid Fisher matrix.

Spectra.setup_geometry()

Set up geometry including active pixels and pointing vectors.

Returns:

Tuple containing (active_pixels, pointing_vectors) where: - active_pixels: Object array of active pixel indices per field component - pointing_vectors: Tuple of unit vector arrays for each component

Return type:

tuple of numpy.ndarray

Raises:

ValueError – If fields have not been set up before calling this method.

Notes

This method: 1. Extracts active pixel counts from each field 2. Computes 3D pointing vectors for each active pixel 3. Computes theta (colatitude) and phi (longitude) for computation basis 4. Sets pointing vectors in the field collection 5. Optionally outputs geometry data to file

The pointing vectors are unit vectors in 3D space pointing to pixel centers, used for spherical harmonic transforms and coordinate operations.

Core Computation

Spectra.setup_maps()[source]

Read observational map data from FITS files.

Loads maps with proper pixel selection, field extraction, and calibration. For cross-correlation (do_cross=True), loads both primary and secondary maps. Output shape: (n_active_pixels, n_simulations).

Raises:
Spectra.setup_fisher_inversion()[source]

Prepare inverted Fisher matrix for QML estimation.

The Fisher matrix is already beam-smoothed (beam window functions absorbed into derivatives). This method computes F⁻¹, F⁻¹/², and writes results to output files.

Raises:

  • ValueError – If Fisher matrix is not available or singular.

  • LinAlgError – If Fisher matrix inversion fails.

Spectra.setup_qml_computation()[source]

Initialize arrays and variables for QML power spectrum computation.

This method allocates memory for QML estimation results and prepares the computational infrastructure for the main QML calculation phase. Arrays are sized based on the number of multipoles, spectra, and Monte Carlo simulations.

Notes

Initializes the following arrays:

  • qml_results: Storage for QML power spectrum estimates with shape (n_simulations, n_total_parameters) where n_total_parameters = n_spectra * (lmax - 1)

  • qml_noise_bias: Storage for noise bias estimates (auto-correlation only) with shape (n_total_parameters,)

The qml_results array stores the raw quadratic estimates: q̂_l^(sim) = (1/2) * x^(sim)T * E_l * x^(sim)

for each simulation and multipole. These raw estimates will be combined using the inverse Fisher matrix to produce the final optimally-weighted power spectrum estimates.

For auto-correlation analyses, noise bias terms are computed as: bias_l = (1/2) * Tr[N * E_l] where N is the noise covariance matrix and E_l is the quadratic estimator matrix for multipole l.

Memory allocation scales as O(n_sims * n_ell * n_spec) for results and O(n_ell * n_spec) for bias terms, which is typically much smaller than the O(n_pix^2) covariance matrix storage.

Examples

This method is called automatically during computation:

>>> spectra = Spectra("config.yaml")
>>> spectra.setup_qml_computation()
>>> print(f"QML results shape: {spectra.qml_results.shape}")
>>> print(f"Total parameters: {spectra.qml_results.shape[1]}")
Spectra.compute_qml_spectra()[source]

Execute parallel QML power spectrum computation.

Distributes multipole computations across MPI processes (round-robin). For each multipole: computes signal derivative, E-operator, and quadratic estimates for all simulations. Selects compressed or traditional method based on basis_manager availability.

Spectra.compute()[source]

Execute QML power spectrum computation.

Calls setup_qml_computation() to initialize arrays, then compute_qml_spectra() for the main parallel computation.

Spectra.run()[source]

Execute the complete QML power spectrum analysis pipeline.

Pipeline phases: 1. Master setup: fields, geometry, covariance, spectra, beams

(reuses Fisher components if provided)

  1. QML setup: maps, Fisher inversion

  2. MPI broadcast of shared data to workers

  3. Parallel QML computation

  4. MPI synchronization

Raises:

Results Processing

Spectra.get_power_spectra(mode='deconvolved', *, as_dict=False)[source]

Retrieve power spectrum estimates in specified normalization mode.

This method returns power spectrum estimates with three different normalization prescriptions, allowing users to choose the output format best suited to their analysis needs.

Parameters:

  • mode (str, optional) –

    Normalization mode for output (default: “deconvolved”):

    • ”deconvolved”: F⁻¹y - estimates of true C_ℓ with correlated errors. This is the standard QML output that attempts to recover the true underlying spectrum by inverting the window function.

    • ”decorrelated”: F⁻¹/²y - uncorrelated bandpower estimates with identity covariance matrix. Useful when independent error bars are needed.

    • ”convolved”: Raw y estimates with window matrix W for theory comparison. Instead of deconvolving the window, compare with window-convolved theory: <y> = W @ C_true.

  • as_dict (bool, optional) – If False (default, back-compat), return the flat numpy array described below. If True, return a dict keyed by physical labels ("TT", "EE", "TE", …) whose values have shape (n_simulations, n_bins). In "convolved" mode, the dict applies to the raw estimates y only; the window matrix W and convolve_theory callable are unchanged (use Fisher.get_bandpower_slices() on self.fisher_instance to navigate W).

Return type:

ndarray | tuple[ndarray, ndarray, Callable] | dict | None

Returns:

  • numpy.ndarray or tuple or dict or None – For “deconvolved” or “decorrelated” modes with as_dict=False:

    Array of shape (n_simulations, n_parameters) where n_parameters = n_spectra * n_bins.

    For “convolved” mode with as_dict=False:

    Tuple of (y, W, convolve_theory_func) where:

    • y: Raw estimates of shape (n_simulations, n_parameters)

    • W: Window matrix of shape (n_parameters, n_parameters)

    • convolve_theory_func: Callable that applies W @ theory

    With as_dict=True: a dict for the per-spectrum arrays (and, in convolved mode, a 3-tuple (y_dict, W, convolve)).

    Returns None for worker processes (rank != 0) or if computation has not completed.

  • **Layout (flat-array shape)**

  • The columns of the returned array are organised as n_spectra

  • contiguous blocks of n_bins columns each. Spectrum order

  • follows Fisher.spectra_list on

  • self.fisher_instance — auto-spectra first (per component,

  • in the order each component emits them), then cross-spectra in

  • (i, j) order with i < j. For TQU (spins=[0, 2],

  • labels=["T", "E", "B"]) under SymmetryMode.SYMMETRIC

  • with nbins=10, n_parameters = 60 and the columns are:: – cols 0..9 : TT cols 10..19 : EE cols 20..29 : BB cols 30..39 : EB cols 40..49 : TE cols 50..59 : TB

  • Use Fisher.get_bandpower_slices() (on

  • self.fisher_instance) or pass as_dict=True to look up

  • per-spectrum slices without computing these offsets by hand.

Raises:

ValueError – If mode is not one of “deconvolved”, “decorrelated”, “convolved”.

Notes

Mode Comparison:

Mode | Formula | Covariance | Use Case |

|------|———|------------|———-| | deconvolved | F⁻¹y | F⁻¹ (correlated) | Standard analysis | | decorrelated | F⁻¹/²y | I (identity) | Independent errors | | convolved | y | F | Theory comparison |

Deconvolved Mode (default): The standard QML output that inverts the window function to recover estimates of the true C_ℓ. Errors are correlated between multipoles.

Decorrelated Mode: Produces bandpower estimates with uncorrelated errors (unit variance). Useful for plotting error bars or when independent measurements are required. Note that information content is preserved but redistributed.

Convolved Mode: Returns raw QML estimates without deconvolution, along with the window matrix for comparing with convolved theoretical spectra. This avoids numerical issues from inverting poorly-conditioned window matrices.

Output Convention: When output_convention is set to "Dl" in the configuration, all returned spectra are converted from C_ℓ to D_ℓ = ℓ(ℓ+1)/(2π) C_ℓ. In convolved mode, the window matrix is also transformed so that the returned convolve_theory_func expects D_ℓ input.

Examples

Default (flat array) — back-compat:

>>> spectra = Spectra("config.yaml")
>>> spectra.run()
>>> cl_deconv = spectra.get_power_spectra()        # shape (n_sims, 60)
>>> ee_slice = cl_deconv[:, 10:20]                 # bins 0..9 of EE

Label-keyed dict (recommended for human-readable code):

>>> cl = spectra.get_power_spectra(as_dict=True)
>>> cl["EE"].shape                                 # (n_sims, n_bins)
>>> cl["TE"]                                       # cross-spectrum

Decorrelated bandpowers:

>>> cl_decorr = spectra.get_power_spectra(mode="decorrelated")
>>> # Errors are all 1.0 by construction

Convolved mode for theory comparison:

>>> y, W, convolve = spectra.get_power_spectra(mode="convolved")
>>> cl_theory_convolved = convolve(cl_theory)
>>> # Compare: y should match cl_theory_convolved

>>> # With as_dict=True, the y component becomes a dict; W and
>>> # the convolve callable are unchanged.
>>> y_dict, W, convolve = spectra.get_power_spectra(
...     mode="convolved", as_dict=True
... )

See also

get_covariance

Get covariance matrix for specified mode

get_error_bars

Get 1-sigma error bars for specified mode

convolve_theory

Apply window matrix to theoretical spectrum

Fisher.get_bandpower_slices

Slice mapping for navigating the flat layout.

Spectra.get_noise_bias()[source]

Retrieve noise bias estimates for auto-correlation spectra.

Returns:

Shape (n_params,). Returns None for workers, cross-correlation analyses (do_cross=True), or if computation incomplete.

Return type:

np.ndarray or None

Notes

Noise bias: (1/2) * Tr[N * E_l]. Cross-correlations are bias-free.

Spectra.get_covariance(mode='deconvolved')[source]

Get covariance matrix for power spectrum estimates in specified mode.

Returns the covariance matrix appropriate for the normalization mode, which quantifies the statistical uncertainties and correlations between different multipole estimates.

Parameters:

mode (str, optional) –

Normalization mode (default: “deconvolved”):

  • ”deconvolved”: Returns F⁻¹, the inverse Fisher matrix with correlated errors between multipoles.

  • ”decorrelated”: Returns identity matrix I, as decorrelated estimates have unit variance by construction.

  • ”convolved”: Returns F, the normalized Fisher matrix, which is the covariance of raw y estimates.

Returns:

Covariance matrix of shape (n_params, n_params). Returns None for worker processes (rank != 0) or if matrices are not available.

Return type:

numpy.ndarray or None

Raises:

ValueError – If mode is not one of “deconvolved”, “decorrelated”, “convolved”.

Notes

Mode-specific covariances:

Mode | Covariance | Interpretation |

|------|————|----------------| | deconvolved | F⁻¹ | Correlated errors on C_ℓ estimates | | decorrelated | I | Unit variance (uncorrelated) | | convolved | F | Covariance of raw y estimates |

The covariance matrix is essential for: - Computing chi-square statistics - Parameter estimation from power spectra - Constructing likelihood functions - Understanding multipole correlations

Examples

Get covariance for chi-square computation:

>>> cov = spectra.get_covariance(mode="deconvolved")
>>> cl_data = spectra.get_power_spectra(mode="deconvolved")
>>> residuals = np.mean(cl_data, axis=0) - cl_theory
>>> chi2 = residuals @ np.linalg.inv(cov) @ residuals

Decorrelated mode has identity covariance:

>>> cov_decorr = spectra.get_covariance(mode="decorrelated")
>>> # cov_decorr is identity matrix

See also

get_power_spectra

Get power spectrum estimates

get_error_bars

Get 1-sigma error bars (sqrt of diagonal)

Spectra.get_error_bars(mode='deconvolved')[source]

Get 1-sigma error bars for power spectrum estimates.

Returns the marginal standard deviations for each multipole estimate, computed from the diagonal of the covariance matrix.

Parameters:

mode (str, optional) –

Normalization mode (default: “deconvolved”):

  • ”deconvolved”: sqrt(diag(F⁻¹)) - standard QML errors

  • ”decorrelated”: all ones (unit variance by construction)

  • ”convolved”: sqrt(diag(F)) - errors on raw y estimates

Returns:

1D array of error bars with shape (n_params,). Returns None for worker processes (rank != 0) or if covariance not available.

Return type:

numpy.ndarray or None

Raises:

ValueError – If mode is not one of “deconvolved”, “decorrelated”, “convolved”.

Notes

The error bars are computed as: σ_i = sqrt(Cov_ii)

where Cov is the mode-appropriate covariance matrix. These represent marginal uncertainties on individual estimates, marginalized over all other multipoles.

Important: For deconvolved and convolved modes, errors are correlated between multipoles. The full covariance matrix (from get_covariance) should be used for proper statistical analysis.

For decorrelated mode, error bars are all 1.0 by construction, making them suitable for simple error bar plots.

Examples

Get error bars for plotting:

>>> errors = spectra.get_error_bars(mode="deconvolved")
>>> cl = np.mean(spectra.get_power_spectra(), axis=0)
>>> plt.errorbar(ell, cl, yerr=errors)

Decorrelated mode has unit errors:

>>> errors_decorr = spectra.get_error_bars(mode="decorrelated")
>>> # All elements are 1.0

See also

get_covariance

Get full covariance matrix

get_power_spectra

Get power spectrum estimates

Spectra.convolve_theory(cl_theory)[source]

Apply window matrix to theoretical power spectrum.

Convolves a theoretical power spectrum with the QML window matrix, producing the expected value of the raw QML estimates y for that theory: <y> = W @ C_theory.

Parameters:

cl_theory (numpy.ndarray) – Theoretical power spectrum values. Should be a 1D array with shape (n_params,) matching the QML output dimensions. When output_convention="Dl", this should be D_ℓ values.

Returns:

Window-convolved theoretical spectrum with shape (n_params,), in the same convention as the input. Returns None if window matrix is not available.

Return type:

numpy.ndarray or None

Notes

This method is useful for comparing QML estimates with theoretical predictions when using the “convolved” normalization mode. Instead of deconvolving the window function from the data (which can be numerically unstable), the window is applied to the theory.

The relationship is: <y> = W @ C_true

where y are the raw QML estimates and W is the window matrix. When output_convention="Dl", the window matrix is internally transformed so the input and output are both in D_ℓ convention.

Examples

Compare convolved estimates with theory:

>>> y, W, _ = spectra.get_power_spectra(mode="convolved")
>>> cl_theory_convolved = spectra.convolve_theory(cl_theory)
>>> residuals = np.mean(y, axis=0) - cl_theory_convolved

See also

get_power_spectra

Get power spectra (convolved mode returns window)

Normalization Modes

Power spectrum estimates can be retrieved in three normalization modes, each suited to different analysis needs. The raw QML estimates \(\mathbf{y}\) are related to the true power spectrum \(C_\ell\) through the Fisher matrix \(\mathbf{F}\):

\[\langle \mathbf{y} \rangle = \mathbf{F} \, \mathbf{C}_\ell\]

The three modes differ in how \(\mathbf{y}\) is post-processed before being returned to the user.

Mode

Formula

Covariance

Use case

"deconvolved"

\(\mathbf{F}^{-1} \mathbf{y}\)

\(\mathbf{F}^{-1}\) (correlated)

Standard analysis: estimates of the true \(C_\ell\)

"decorrelated"

\(\mathbf{F}^{-1/2} \mathbf{y}\)

\(\mathbf{I}\) (identity)

Independent error bars for visualization

"convolved"

\(\mathbf{y}\)

\(\mathbf{F}\)

Theory comparison via \(\langle \mathbf{y} \rangle = \mathbf{W} C_\ell^{\rm th}\)

Deconvolved (default)

The standard QML output. The window function is inverted so the returned estimates are direct estimates of the true \(C_\ell\). Errors are correlated between multipoles; use get_covariance() for the full covariance matrix.

cl = spectra.get_power_spectra()  # or mode="deconvolved"
cov = spectra.get_covariance()
errors = spectra.get_error_bars()

Decorrelated

Applies \(\mathbf{F}^{-1/2}\) so that the output covariance is the identity matrix. Error bars are unity by construction, making this mode useful for quick visualization. Information content is preserved but redistributed.

cl_decorr = spectra.get_power_spectra(mode="decorrelated")
errors = spectra.get_error_bars(mode="decorrelated")  # all 1.0

Convolved

Returns the raw estimates together with the window matrix \(\mathbf{W}\) and a convenience function to convolve a theoretical spectrum. This avoids inverting a potentially ill-conditioned Fisher matrix.

y, W, convolve = spectra.get_power_spectra(mode="convolved")
cl_theory_convolved = convolve(cl_theory)
# Compare: y ≈ cl_theory_convolved

Chi-squared equivalence

All three modes yield the same chi-squared value for any given data realization:

\[\chi^2_{\rm deconv} = (\hat{C} - C^{\rm th})^T \mathbf{F} (\hat{C} - C^{\rm th}) = (\mathbf{y} - \mathbf{W} C^{\rm th})^T \mathbf{F}^{-1} (\mathbf{y} - \mathbf{W} C^{\rm th}) = \| \mathbf{q} - \mathbf{F}^{1/2} C^{\rm th} \|^2\]

Usage Examples

Basic QML Estimation

from cosmoforge.qube import Spectra

# Initialize QML analysis
spectra = Spectra("config/qml_config.yaml")

# Run complete estimation pipeline
spectra.run()

# Get power spectrum estimates (master process)
if spectra.rank == 0:
    cl = spectra.get_power_spectra()
    errors = spectra.get_error_bars()

Cross-Correlation Analysis

# Configure for cross-correlation between two datasets
config = {
    'do_cross': True,
    'covmatfile1': 'dataset1_noise.bin',
    'covmatfile2': 'dataset2_noise.bin',
    'mapfiles1': ['map1_T.fits', 'map1_Q.fits', 'map1_U.fits'],
    'mapfiles2': ['map2_T.fits', 'map2_Q.fits', 'map2_U.fits'],
}

spectra = Spectra(config)
spectra.run()

MPI Parallel Execution

# Run QML estimation with 16 processes
mpirun -n 16 python -c "
from cosmoforge.qube import Spectra
spectra = Spectra('qml_config.yaml')
spectra.run()
"

Error Analysis

# Compute error estimates from Fisher matrix
from cosmoforge.qube import Fisher, Spectra

# First run Fisher analysis for error forecasting
fisher = Fisher("fisher_config.yaml")
fisher.run()

# Then run QML estimation
spectra = Spectra("qml_config.yaml")
spectra.run()

if spectra.rank == 0:
    # Compare with Fisher forecasts
    fisher_errors = fisher.get_error_bars()
    qml_estimates = spectra.get_power_spectra()

Computational Performance

Algorithm Complexity

The QML algorithm scales as:

  • Matrix inversion: \(O(N_{pix}^3)\) one-time cost

  • Estimator matrices: \(O(N_{ell} \times N_{pix}^2)\)

  • Quadratic forms: \(O(N_{ell} \times N_{pix}^2)\) per data realization

where \(N_{pix}\) is the number of active pixels and \(N_{ell}\) is the number of multipole bins.

Memory Requirements

Peak memory usage is dominated by:

  • Covariance matrix: \(N_{pix}^2 \times 8\) bytes

  • Estimator matrices: \(N_{ell} \times N_{pix}^2 \times 8\) bytes

  • Working arrays: Additional \(\sim 2 \times N_{pix}^2 \times 8\) bytes

Optimization Strategies

For large-scale analyses:

  1. Hierarchical approach: Start with low resolution for initial estimates

  2. Block processing: Process multipole blocks separately to reduce memory

  3. Iterative methods: Use conjugate gradient for matrix operations

  4. Preconditioning: Improve convergence with appropriate preconditioners

Resource Scaling

Typical requirements for full-mission analyses:

Resolution

Active Pixels

Memory (GB)

Compute Time

Cores

nside=128

~50k

~2

~30 min

4-8

nside=256

~200k

~8

~3 hours

8-16

nside=512

~800k

~128

~1 day

16-32

nside=1024

~3.2M

~2000

~1 week

32-64

Configuration Options

Essential QML configuration parameters:

# Analysis type
do_cross: false  # Auto-correlation (true for cross-correlation)

# Input data
mapfiles1: ["map_T.fits", "map_Q.fits", "map_U.fits"]
covmatfile1: "noise_covariance.bin"
maskfile: "analysis_mask.fits"
input_convention: Dl   # "Cl" (default) or "Dl" for input spectra
output_convention: Dl  # "Cl" (default) or "Dl" for output spectra

# Multipole range
lmin: 2
lmax: 2000
nbands: 50  # Number of bandpowers

# Output files
outfile: "qml_power_spectra.dat"
outfilelhood: "likelihood_data.dat"

# Computation options
compute_fisher: true  # Include Fisher matrix computation
save_estimators: false  # Save estimator matrices (large files)

Advanced Features

Likelihood Analysis

The QML estimates come with full likelihood information for parameter fitting:

# Access likelihood data for cosmological parameter fitting
spectra = Spectra("config.yaml")
spectra.run()

if spectra.rank == 0:
    # Get bandpowers and covariance for likelihood
    bandpowers = spectra.get_power_spectra()
    covariance = spectra.get_covariance()

    # Use in cosmological parameter fitting pipeline

Systematic Error Studies

# Test robustness to different analysis choices
systematics = {
    'lmax_values': [1500, 2000, 2500, 3000],
    'mask_apodization': [0.0, 1.0, 2.0, 5.0],
    'multipole_binning': ['linear', 'logarithmic', 'hybrid']
}

for test_case in systematic_combinations:
    config = base_config.copy()
    config.update(test_case)

    spectra = Spectra(config)
    spectra.run()
    # Analyze systematic differences...

See Also