pespace.inference.fisher

Fisher matrix computation for the Whittle likelihood.

In the high-SNR regime the likelihood can be expanded as a series in 1 / rho; keeping the leading order, the posterior distribution with uniform priors is approximated by a multivariate Gaussian distribution with means given by the fiducial parameters and covariance matrix given by the inverse of the Fisher information matrix (Vallisneri 2008),

\[F_{ij} = \left\langle \frac{\partial \tilde{s}}{\partial \theta_i} \middle| \frac{\partial \tilde{s}}{\partial \theta_j} \right\rangle,\]

where the noise-weighted inner product is

\[\langle a | b \rangle = 4 \mathrm{Re} \sum_k \frac{a^*(f_k) b(f_k)}{S_n(f_k) T}.\]

Functions

fisher_matrix(likelihood, fiducial_params[, ...])	Compute the Fisher information matrix at one or more fiducial parameter sets.
sample_from_fisher(fisher, num[, rkey_seed])	Draw random samples from the Gaussian approximation N(0, F^{-1}).
wrap_angle_params(samples)	Wrap the angle parameters of a sample set into finite canonical ranges.

pespace.inference.fisher.fisher_matrix(likelihood, fiducial_params, batch_size=None)

Compute the Fisher information matrix at one or more fiducial parameter sets.

The derivatives of the model response with respect to all parameters of the likelihood are computed with forward-mode automatic differentiation, and the Fisher matrix is assembled from the noise-weighted inner products of the derivatives, summed over all detectors and requested channels,

F_ij = 4 * df * Re sum_{c, f} conj(ds_i(f)) * ds_j(f) / psd_c(f).

A single parameter set returns a single matrix; a batch of parameter sets is evaluated with vmap (optionally split into chunks by batch_size), returning one matrix per set.

Parameters:

• likelihood (WhittleLikelihoodJax) – Likelihood instance defining the waveform model, the detectors with their noise power spectral densities, and the requested TDI channels. The observed data of the detectors is not used: the Fisher matrix only depends on the model and the noise power spectral densities.

• fiducial_params (dict) – Fiducial parameters at which the Fisher matrix is evaluated, e.g. the injected parameters of the simulated signal. Must contain all parameters of likelihood.likelihood_param_type (the waveform parameters followed by the response parameters); additional entries are ignored. A single parameter set is given by scalar values; a batch of parameter sets is given by 1-D arrays whose (common) leading axis enumerates the sets, with scalar entries broadcast over the batch.

• batch_size (int or None, optional) – Number of parameter sets evaluated at once within vmap. None (default) evaluates the whole batch with a single vmap call; an integer evaluates the batch in chunks of that size, bounding the memory footprint of large batches. Ignored for a single parameter set.

Returns:

Fisher matrix of shape (N_params, N_params) for a single parameter set, or (N_batch, N_params, N_params) for a batch of N_batch parameter sets. N_params is the number of parameters of likelihood.likelihood_param_type, and the rows and columns follow the order of its fields.

Return type:

jnp.ndarray

Notes

• The posterior distribution is approximated by a multivariate Gaussian with mean fiducial_params and covariance matrix jnp.linalg.inv(fisher_matrix(...)).

• When the observed data are noiseless and equal to the response of fiducial_params, the Fisher matrix equals minus the Hessian of log_likelihood() evaluated at fiducial_params.

pespace.inference.fisher.sample_from_fisher(fisher, num, rkey_seed=0)

Draw random samples from the Gaussian approximation N(0, F^{-1}).

The high-SNR posterior approximation associated with a Fisher matrix F is a multivariate Gaussian whose mean is given by the fiducial parameters and whose covariance matrix is F^{-1} (see fisher_matrix()). This function draws zero-mean samples from that Gaussian; the fiducial parameter values must be added by the caller.

The samples are obtained without forming or inverting the covariance matrix. The Fisher matrix is first rescaled to unit diagonal with the dimensionless matrix F' = F / (sqrt(diag(F)) outer sqrt(diag(F))), which removes the disparity of the parameter units (for a representative mixed-unit gravitational-wave case the condition number drops from ~1e19 to ~1e8, leaving a comfortable margin within float64). The scaled matrix is factorized with a Cholesky decomposition, F' = L L^T, and each sample is computed by solving the triangular system L^T x' = z for standard normal z, so that Cov(x') = L^{-T} L^{-1} = F'^{-1}; the samples are then rescaled back to the original units. The solve is backward stable and avoids both the explicit matrix inversion and the loss of symmetry and positive definiteness that make a naive inv(F)-then-sample route unreliable at such condition numbers.

Parameters:

• fisher (array_like, shape (N, N)) – Fisher information matrix, symmetric positive definite, with rows and columns ordered as the parameters of the likelihood.

• num (int) – Number of samples to draw.

• rkey_seed (int, optional) – Seed used to draw the standard normal variates. Defaults to 0, so repeated calls return identical samples.

Returns:

Zero-mean samples with covariance matrix fisher^{-1}, in the same parameter order and units as fisher.

Return type:

jnp.ndarray, shape (num, N)

Notes

• Double precision should be enabled (jax_enable_x64); in float32 the diagonal rescaling cannot compensate the extreme condition number of the mixed-unit Fisher matrix.

• The function is compatible with jax.jit; jax.vmap over the Fisher matrices (with a common seed) draws samples for several parameter sets at once.

pespace.inference.fisher.wrap_angle_params(samples)

Wrap the angle parameters of a sample set into finite canonical ranges.

The following entries of samples are wrapped if present; all other entries are left unchanged:

• ecliptic_longitude: wrapped modulo 2*pi into [0, 2*pi);

• ecliptic_latitude: folded by reflection about the boundaries -pi/2 and pi/2 into [-pi/2, pi/2];

• inclination: folded by reflection about the boundaries 0 and pi into [0, pi];

• polarization: wrapped modulo pi into [0, pi) (the detector response depends on the polarization angle through its double angle, so pi is its fundamental period);

• reference_phase: wrapped modulo 2*pi into [0, 2*pi).

These are the finite ranges conventionally used for the angle parameters of the space-borne detector response; they match the sampling ranges adopted in the example scripts.

Parameters:

samples (mapping) – Mapping from parameter names to sample values, e.g. a dict of arrays or a pandas.DataFrame; values may be scalars or arrays of any array library (Python, NumPy or JAX).

Returns:

A copy of samples with the angle parameters restricted to their finite ranges; neither the input container nor its arrays are modified. The copy is shallow: the wrapped entries are replaced by newly created arrays, while the remaining entries are shared with the input (a deep copy is not needed, since no array is modified in place).

Return type:

mapping

Notes

• The wrapping is exact for the periodic angles (ecliptic_longitude, polarization and reference_phase) and is a reflection fold for the bounded angles (ecliptic_latitude and inclination).

• In this likelihood the reference_phase is periodic with period pi rather than 2*pi (the strain depends on 2 * reference_phase), so the two wells at phi and phi + pi are exactly symmetric; the wrapping above keeps the conventional [0, 2*pi) range, and the wells can be identified with reference_phase % pi when needed.

• The wrapping is compatible with jax.jit; the reflection folds use jnp.abs, so those results are JAX arrays even for non-JAX inputs.