Fisher Matrix Analysis

The Fisher class provides comprehensive Fisher information matrix computation for cosmological parameter forecasting and optimal experiment design. The implementation features MPI parallelization, exact covariance matrix treatment, and integration with the broader CosmoForge analysis pipeline.

Mathematical Background

The Fisher information matrix quantifies the information content of observations about model parameters. For observations of Gaussian fields on the sphere, it is computed as:

\[F_{ij} = \frac{1}{2} \text{Tr}\left[ \mathbf{C}^{-1} \frac{\partial \mathbf{C}}{\partial \theta_i} \mathbf{C}^{-1} \frac{\partial \mathbf{C}}{\partial \theta_j} \right]\]

where: - \(\mathbf{C} = \mathbf{S} + \mathbf{N}\) is the total covariance matrix - \(\mathbf{S}\) is the signal covariance from theoretical power spectra - \(\mathbf{N}\) is the instrumental noise covariance - \(\theta_i\) are the cosmological parameters of interest

The inverse Fisher matrix provides the parameter covariance matrix under Gaussian posterior assumptions: \(\text{Cov}(\theta_i, \theta_j) = (F^{-1})_{ij}\).

Class Documentation

class qube.fisher.Fisher(params_file=None, compression=None, cache_derivatives=False, symmetry_mode=None, **kwargs)[source]

Bases: Core, MPISharedMemoryMixin

Fisher matrix computation for cosmological parameter estimation.

This class implements the Fisher information matrix calculation for forecasting cosmological parameter constraints from observations of spin-0 and spin-2 fields on the sphere. It inherits from cosmocore.Core and extends it with Fisher-specific functionality.

The Fisher matrix F_ij quantifies the information content of the data about parameters θ_i and θ_j, computed as:

F_ij = (1/2) * Tr[C^(-1) * ∂C/∂θ_i * C^(-1) * ∂C/∂θ_j]

The computation uses Core’s unified API which transparently handles harmonic and pixel computation bases when configured.

Parameters:

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

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

  • **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

Computed Fisher information matrix.

Type:

numpy.ndarray

n_ell

Number of multipole moments in the inference window (params.lmax - params.lmin + 1).

Type:

int

n_params

Total number of Fisher matrix parameters (nspectra * n_ell).

Type:

int

Examples

Basic Fisher matrix computation:

>>> from cosmoforge.qube import Fisher
>>> fisher = Fisher("config/fisher_analysis.yaml")
>>> fisher.run()
>>> F_matrix = fisher.get_fisher_matrix()

With harmonic basis:

>>> fisher = Fisher("config.yaml", compression={"method": "harmonic"})
>>> fisher.run()  # Uses harmonic basis transparently
__init__(params_file=None, compression=None, cache_derivatives=False, symmetry_mode=None, **kwargs)[source]

Initialize Fisher matrix computation class.

Parameters:

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

  • compression (dict, optional) –

    Computation basis configuration dictionary. Options:

    • method : str (“harmonic” or “pixel”)

    • epsilon : float (eigenvalue threshold for pixel basis)

    • basis : str (for pixel: “harmonic”, “noise_weighted”, etc.)

    • mode_fraction : float (alternative to epsilon)

  • cache_derivatives (bool, optional) – Whether to cache binned derivative matrices for Spectra to reuse. Default is False — Spectra recomputes via the pixel-direct fast path (single Legendre/Wigner pass per bin), which is amortized across MPI ranks. Setting True keeps the dense per-bin cache, which can cost n_params * n_pix^2 of memory (e.g. ~14 GB at QU_nside64 fsky=0.1) but skips the recompute. The harmonic-fast path always caches its small COO triplets regardless of this flag.

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

property lmax_signal: int

Signal-cov ceiling (ADR 0009).

Resolution order: explicit setter, then params.lmax_signal, then 4 * nside.

setup_signal_matrix()[source]

Compute signal covariance S = Sum_ell C_ell P_ell from fiducial spectra.

Returns:

Signal covariance matrix of shape (n_pix, n_pix).

Return type:

numpy.ndarray

prepare_covariance_matrices()[source]

Build total covariance C = N + S, then compute C^{-1}.

After this method, noise_cov1 (and noise_cov2 if cross-correlation) hold the inverted total covariance, not the original noise.

compute()[source]

Compute Fisher matrix F_{ij} = (1/2) Tr[C⁻¹ dC^i C⁻¹ dC^j].

Supports single- and multi-spectrum, with or without computation basis. For harmonic basis, exploits sparse derivative structure for O(nnz_i × nnz_j) trace computation instead of O(n_modes²).

run()[source]

Execute the complete Fisher matrix analysis pipeline.

A computation basis can be enabled via the compression parameter.

Return type:

None

get_fisher_matrix()[source]

Retrieve the beam-smoothed Fisher information matrix.

Return type:

ndarray | None

Layout

Shape (nspectra * nbins, nspectra * nbins). Rows and columns run identically: each spectrum occupies a contiguous block of nbins columns; within a block, columns are ordered by bin index (bin 0 first, then bin 1, …). Spectrum order follows spectra_list, which is filled by the harmonic layer as:

  1. All auto-spectra in component order (component i first, then the kinds emitted by that component — e.g. for a spin-2 component: GG (EE), CC (BB), GC (EB)).

  2. All cross-spectra in lexicographic (i, j) order with i < j, each pair contributing the kinds emitted by the spin pair (e.g. spin-0 × spin-2 emits SG (TE), SC (TB)).

Example: for TQU (spins=[0, 2], labels=["T", "E", "B"]) under SymmetryMode.SYMMETRIC and nbins=10, the matrix is (60, 60) and the column blocks 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 get_bandpower_slices() to look up the slice for a given physical-label spectrum without having to compute these offsets by hand:

slices = fisher.get_bandpower_slices()
tt_block = fisher.get_fisher_matrix()[slices["TT"], slices["TT"]]
te_x_ee = fisher.get_fisher_matrix()[slices["TE"], slices["EE"]]
get_error_bars(*, as_dict=False)[source]

Compute parameter forecast errors from the Fisher matrix.

Parameters:

as_dict (bool, optional) – If False (default, back-compat), return a flat numpy array of shape (nspectra * nbins,) ordered as described in get_fisher_matrix(). If True, return a dict keyed by physical labels ("TT", "EE", "TE", …) with per-spectrum error arrays of shape (nbins,).

Return type:

ndarray | dict[str, ndarray] | None

Examples

Flat shape (default) — same as before:

err = fisher.get_error_bars()              # shape (60,) for TQU
ee_errors = err[10:20]                     # bins 0..9 of EE

Dict shape — label-keyed, no manual offsets:

err = fisher.get_error_bars(as_dict=True)
ee_errors = err["EE"]                      # shape (nbins,)
te_errors = err["TE"]
get_bandpower_slices()[source]

Return a mapping {label: slice} for the flat Fisher axis.

Each value is a slice selecting that spectrum’s bin block in get_fisher_matrix() (along either axis) or in get_error_bars() (the single axis). Labels are concatenated from params.labels per slot (e.g. "TT", "EE", "TE"); for multi-frequency setups the user’s per-slot labels make the keys unambiguous ("T100T143" etc.).

Returns None until compute() has run (spectra_list is populated as a side effect of the Fisher pipeline).

Return type:

dict[str, slice] | None

Examples

For TQU under SYMMETRIC mode with nbins=10 the returned mapping is:

{"TT": slice(0, 10),
 "EE": slice(10, 20),
 "BB": slice(20, 30),
 "EB": slice(30, 40),
 "TE": slice(40, 50),
 "TB": slice(50, 60)}

Typical usage:

F = fisher.get_fisher_matrix()
slices = fisher.get_bandpower_slices()
tt_diag = np.diag(F[slices["TT"], slices["TT"]])
te_x_ee = F[slices["TE"], slices["EE"]]  # cross-spectrum block

See also

get_fisher_block

Convenience wrapper that returns a single (label_i, label_j) block directly.

get_fisher_block(label_i, label_j=None)[source]

Extract one (label_i, label_j) block of the Fisher matrix.

Convenience wrapper over get_fisher_matrix() and get_bandpower_slices().

Parameters:

  • label_i (str) – Physical label for the row block (e.g. "TT", "EE").

  • label_j (str or None, optional) – Physical label for the column block. If None (default), returns the auto block (label_i, label_i).

Returns:

Block of shape (nbins, nbins). Returns None on worker ranks or before compute() has run.

Return type:

numpy.ndarray or None

Raises:

KeyError – If either label is not in the current get_bandpower_slices() mapping. The error message lists the legal labels.

Examples

>>> tt_block = fisher.get_fisher_block("TT")
>>> te_x_ee = fisher.get_fisher_block("TE", "EE")
>>> bb_eb = fisher.get_fisher_block("BB", "EB")  # likely zero for std cosmo
get_window_matrix()[source]

Retrieve the window matrix for QML power spectrum estimation.

The window matrix W relates the expected QML estimates to the beam-smoothed theory spectrum: <q> = W @ C_theory. It encodes the mode coupling induced by partial sky coverage, beam, and pixel window effects.

Layout

Shape (n_params, n_params) with the same flat layout as get_fisher_matrix() — each spectrum occupies a contiguous block of nbins rows/columns, in spectra_list order. When using this matrix in an inference loop, the theory vector C_theory must follow the same flat ordering; see get_bandpower_slices() or pass a label-keyed dict to the convolve_theory_func returned by Spectra.get_power_spectra() (mode="convolved").

returns:

Window matrix of shape (n_params, n_params) where n_params = n_spectra * n_bins. Returns None for worker processes or if computation hasn’t completed.

rtype:

numpy.ndarray or None

Notes

The window matrix is the beam-smoothed Fisher matrix:

W_{bb'} = (1/2) Tr[C⁻¹ dC^b C⁻¹ dC^{b'}]

where dC^b = Sum_ell w_{b,ell} b²_ell dC^ell includes the binning weights and beam smoothing factors.

Used by the “convolved” normalization mode, where instead of deconvolving the window, the theory is convolved for comparison.

get_bandpower_window_function(per_ell_fisher=None)[source]

Compute the bandpower window function for binned QML inference.

The window matrix W maps a per-ℓ theory spectrum to the expected binned bandpower estimates from the deconvolved QML output:

<C_hat_b> = W @ C_ell_theory

Use this in a Gaussian likelihood for parameter inference:

mu_b(θ) = W @ C_ell(θ) -2 ln L = (d - mu)ᵀ F_b (d - mu)

where d = Spectra.get_power_spectra(mode="deconvolved") and F_b = Fisher.get_fisher_matrix().

This differs from get_window_matrix(), which returns the bin-to-bin Fisher used in “convolved” mode and assumes a constant spectrum within each bin. The bandpower window function is the correct quantity for inference when bins are wide enough that the theory varies appreciably within them.

Parameters:

per_ell_fisher (np.ndarray, optional) – Pre-computed per-ℓ Fisher matrix of shape (n_ell, n_ell) (single spectrum). If None, the per-ℓ Fisher is computed on demand by re-running the Fisher pipeline with delta_ell=1 using the current configuration. The result is cached for subsequent calls.

Returns:

W – Bandpower window matrix. Rows correspond to bin indices, columns to multipoles ℓ=2..lmax. Returns None on worker processes (rank != 0).

Return type:

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

Notes

The window function encodes:

  • Mode coupling from sky cuts.

  • Beam smoothing (b²_ℓ).

  • Pixel window functions.

  • The Fisher-weighting that QML naturally applies within each bin.

The theory C_ℓ should be the unbeamed physical spectrum (standard CAMB/CLASS output), since beam² is already absorbed into the derivatives used to build the window.

For best statistical accuracy near lmax, use the buffer approach: estimate to lmax_buffer = lmax_science + margin (a few bin widths or ~1/fsky multipoles), then drop the buffer rows from W, d and F_b before forming the likelihood. This absorbs mode coupling from ℓ > lmax_science into discarded buffer bins.

Scope: currently supports single-spectrum analysis only — calling it on a multi-spectrum Fisher (nspectra > 1) raises NotImplementedError. For multi-spectrum likelihoods today, use get_window_matrix() together with the convolve_theory_func returned by Spectra.get_power_spectra() (mode="convolved"); that path accepts label-keyed dict input and applies the multi- spectrum window matrix. Multi-spectrum extension of this per-ℓ window function is planned.

Examples

Standard inference workflow:

>>> fisher = Fisher("config.yaml", compression={"method": "harmonic"})
>>> fisher.set_binning(bins)
>>> fisher.run()
>>> W = fisher.get_bandpower_window_function()
>>> F_b = fisher.get_fisher_matrix()
>>> # In the MCMC loop:
>>> mu = W @ cl_theory_per_ell        # expected bandpower
>>> r = data_bandpower - mu           # residual
>>> neg2lnL = r @ F_b @ r             # quadratic form

See also

get_window_matrix

Bin-to-bin window for convolved mode.

get_fisher_matrix

Inverse covariance for the bandpower likelihood.

Key Methods

Initialization and Setup

Fisher.__init__(params_file=None, compression=None, cache_derivatives=False, symmetry_mode=None, **kwargs)[source]

Initialize Fisher matrix computation class.

Parameters:

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

  • compression (dict, optional) –

    Computation basis configuration dictionary. Options:

    • method : str (“harmonic” or “pixel”)

    • epsilon : float (eigenvalue threshold for pixel basis)

    • basis : str (for pixel: “harmonic”, “noise_weighted”, etc.)

    • mode_fraction : float (alternative to epsilon)

  • cache_derivatives (bool, optional) – Whether to cache binned derivative matrices for Spectra to reuse. Default is False — Spectra recomputes via the pixel-direct fast path (single Legendre/Wigner pass per bin), which is amortized across MPI ranks. Setting True keeps the dense per-bin cache, which can cost n_params * n_pix^2 of memory (e.g. ~14 GB at QU_nside64 fsky=0.1) but skips the recompute. The harmonic-fast path always caches its small COO triplets regardless of this flag.

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

Fisher.setup_signal_matrix()[source]

Compute signal covariance S = Sum_ell C_ell P_ell from fiducial spectra.

Returns:

Signal covariance matrix of shape (n_pix, n_pix).

Return type:

numpy.ndarray

Fisher.prepare_covariance_matrices()[source]

Build total covariance C = N + S, then compute C^{-1}.

After this method, noise_cov1 (and noise_cov2 if cross-correlation) hold the inverted total covariance, not the original noise.

Core Computation

Fisher.compute()[source]

Compute Fisher matrix F_{ij} = (1/2) Tr[C⁻¹ dC^i C⁻¹ dC^j].

Supports single- and multi-spectrum, with or without computation basis. For harmonic basis, exploits sparse derivative structure for O(nnz_i × nnz_j) trace computation instead of O(n_modes²).

Fisher.run()[source]

Execute the complete Fisher matrix analysis pipeline.

A computation basis can be enabled via the compression parameter.

Return type:

None

Results Retrieval

Fisher.get_fisher_matrix()[source]

Retrieve the beam-smoothed Fisher information matrix.

Return type:

ndarray | None

Layout

Shape (nspectra * nbins, nspectra * nbins). Rows and columns run identically: each spectrum occupies a contiguous block of nbins columns; within a block, columns are ordered by bin index (bin 0 first, then bin 1, …). Spectrum order follows spectra_list, which is filled by the harmonic layer as:

  1. All auto-spectra in component order (component i first, then the kinds emitted by that component — e.g. for a spin-2 component: GG (EE), CC (BB), GC (EB)).

  2. All cross-spectra in lexicographic (i, j) order with i < j, each pair contributing the kinds emitted by the spin pair (e.g. spin-0 × spin-2 emits SG (TE), SC (TB)).

Example: for TQU (spins=[0, 2], labels=["T", "E", "B"]) under SymmetryMode.SYMMETRIC and nbins=10, the matrix is (60, 60) and the column blocks 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 get_bandpower_slices() to look up the slice for a given physical-label spectrum without having to compute these offsets by hand:

slices = fisher.get_bandpower_slices()
tt_block = fisher.get_fisher_matrix()[slices["TT"], slices["TT"]]
te_x_ee = fisher.get_fisher_matrix()[slices["TE"], slices["EE"]]
Fisher.get_error_bars(*, as_dict=False)[source]

Compute parameter forecast errors from the Fisher matrix.

Parameters:

as_dict (bool, optional) – If False (default, back-compat), return a flat numpy array of shape (nspectra * nbins,) ordered as described in get_fisher_matrix(). If True, return a dict keyed by physical labels ("TT", "EE", "TE", …) with per-spectrum error arrays of shape (nbins,).

Return type:

ndarray | dict[str, ndarray] | None

Examples

Flat shape (default) — same as before:

err = fisher.get_error_bars()              # shape (60,) for TQU
ee_errors = err[10:20]                     # bins 0..9 of EE

Dict shape — label-keyed, no manual offsets:

err = fisher.get_error_bars(as_dict=True)
ee_errors = err["EE"]                      # shape (nbins,)
te_errors = err["TE"]
Fisher.get_window_matrix()[source]

Retrieve the window matrix for QML power spectrum estimation.

The window matrix W relates the expected QML estimates to the beam-smoothed theory spectrum: <q> = W @ C_theory. It encodes the mode coupling induced by partial sky coverage, beam, and pixel window effects.

Layout

Shape (n_params, n_params) with the same flat layout as get_fisher_matrix() — each spectrum occupies a contiguous block of nbins rows/columns, in spectra_list order. When using this matrix in an inference loop, the theory vector C_theory must follow the same flat ordering; see get_bandpower_slices() or pass a label-keyed dict to the convolve_theory_func returned by Spectra.get_power_spectra() (mode="convolved").

returns:

Window matrix of shape (n_params, n_params) where n_params = n_spectra * n_bins. Returns None for worker processes or if computation hasn’t completed.

rtype:

numpy.ndarray or None

Notes

The window matrix is the beam-smoothed Fisher matrix:

W_{bb'} = (1/2) Tr[C⁻¹ dC^b C⁻¹ dC^{b'}]

where dC^b = Sum_ell w_{b,ell} b²_ell dC^ell includes the binning weights and beam smoothing factors.

Used by the “convolved” normalization mode, where instead of deconvolving the window, the theory is convolved for comparison.

Usage Examples

Basic Fisher Matrix Computation

from cosmoforge.qube import Fisher

# Initialize with configuration file
fisher = Fisher("config/fisher_analysis.yaml")

# Run complete computation pipeline
fisher.run()

# Get results (master process only)
if fisher.rank == 0:
    F_matrix = fisher.get_fisher_matrix()
    param_errors = fisher.get_error_bars()

    print(f"Fisher matrix shape: {F_matrix.shape}")
    print(f"Parameter errors: {param_errors}")

MPI Parallel Execution

# Run with 8 MPI processes
mpirun -n 8 python -c "
from cosmoforge.qube import Fisher
fisher = Fisher('config.yaml')
fisher.run()
"

Parameter Forecasting

# Compute parameter constraints for different survey configurations
surveys = ['current', 'near_future', 'ultimate']

for survey in surveys:
    fisher = Fisher(f"config/{survey}_config.yaml")
    fisher.run()

    if fisher.rank == 0:
        errors = fisher.get_error_bars()
        # Analyze parameter constraints...

Computational Notes

Performance Characteristics

  • Scaling: \(O(N_{ell}^2 \times N_{pix}^3)\) where \(N_{ell}\) is the number of multipole parameters and \(N_{pix}\) is the number of active pixels

  • Memory: Dominated by covariance matrix storage \(O(N_{pix}^2)\)

  • Parallelization: Fisher matrix elements distributed across MPI processes

Optimization Strategies

  • Use appropriate multipole range (lmin, lmax) for target precision

  • Leverage sky cuts to reduce active pixel count

  • Distribute computation across multiple nodes for large analyses

  • Monitor memory usage for high-resolution pixelizations

Typical Resource Requirements

For temperature + polarization analysis:

HEALPix nside

Active Pixels

Memory (GB)

Compute Time

Recommended Cores

256

~200k

~4

~1 hour

4-8

512

~800k

~16

~8 hours

8-16

1024

~3.2M

~250

~2 days

16-32

2048

~12.8M

~4000

~2 weeks

32-64

Configuration Parameters

Key YAML configuration parameters for Fisher analysis:

# HEALPix parameters
nside: 512
lmin: 2
lmax: 3000

# Field configuration
nfields: 3  # T, Q, U
physical_labels: ["T", "Q", "U"]

# Analysis mode
do_cross: false  # Auto-correlation analysis

# Input files
covmatfile1: "noise_cov_primary.bin"
clfile: "fiducial_spectra.txt"
beamfile: "beam_transfer.fits"

# Output files
outfilefisher: "fisher_matrix.dat"
outinvcovmatfile1: "inv_cov_primary.bin"

See Also

  • qube.spectra.Spectra : QML power spectrum estimation using Fisher results

  • cosmocore : Underlying computational infrastructure

  • Fisher Matrix Analysis Pipeline : Main execution script for Fisher analysis