torchsky.map.operators.nuisance¶

Bases: _NuisancePrecision

Woodbury precision for externally calibrated Gaussian random effects.

A calibration independent of the read/noise term supplies coefficient mean mean and covariance Sigma = L L.T via covariance_factor=L of shape (n_modes, calibration_rank). Rank-deficient or zero covariance is supported directly, without inverting Sigma. With F=B L, apply W0 - W0 F solve(I + F.T W0 F, F.T W0). Calibration here declares a random-effect uncertainty model, not an arbitrary penalty or sky prior.

Call center(data) before passing data to normal-equation assembly; apply and quadratic always consume an already-centered residual. logdet_covariance_correction is logdet(I + F.T W0 F); add it to the base observed covariance log determinant when constructing a Gaussian likelihood. The base normalization remains the caller's responsibility. For parameter-dependent covariance the quadratic alone is insufficient.

The callback, mask, dtype/device, factor lifetime and observed-support contracts are the same as for ProfiledNuisancePrecision. In its factor construction/application costs, replace k by the number of columns of L; forming B L additionally costs O(N k r) for r calibration-factor columns. No jitter is added. Shared calibration across exposures requires one joint basis/precision action, rather than independent per-exposure elimination.

_cholesky: Tensor | None = None

_factor = _flat_basis @ covariance_factor

_mean_image = reshape(data_shape)

_weighted_factor = _weighted_columns(_factor)

logdet_covariance_correction = 2.0 * sum()

treatment = 'calibrated_random_effect'

apply(residual: Tensor) -> Tensor

center(data: Tensor) -> Tensor

Bases: _NuisancePrecision

Precision after profiling fixed unknown linear coefficients.

For y = A x + B beta + epsilon, this applies W0 - W0 B (B.T W0 B)^dagger B.T W0. basis is B with shape (*data_shape, n_modes). The base callback must be a fixed linear, symmetric positive-semidefinite precision with matching shape/dtype/device. Small nuisance-space checks cannot establish that contract outside its span.

valid zeros input and output support. For correlated noise the callback must already represent the precision of the observed measurements: masking a full-data precision does not marginalize missing observations.

rank_rtol defaults to n_modes * finfo(dtype).eps; eigenvalues above max(rank_atol, rank_rtol * max_abs_eigenvalue) are retained. This declared numerical-rank policy never regularizes or lifts zero eigenvalues. Rank depends on basis scaling: discarding a physically relevant nuisance mode can leave apparent precision, so choose tolerances that resolve the declared nuisance span. nuisance_rank reports the retained rank; discarded modes are not certified identifiable. An unmeasured sky/background gauge stays unmeasured. Negative Gram eigenvalues exceeding matrix-product roundoff are rejected, independently of the chosen rank tolerance.

Full-rank blocks use Cholesky solves. Singular blocks use the small Hermitian Moore-Penrose inverse. Autograd is supported on a fixed-rank stratum; rank changes are not differentiable. Factors are built once, so rebuild after changing basis, mask, or callback parameters, and for a fresh autograd graph. Storage is O(N k + k^2), construction uses k base actions and O(N k^2) contractions, and each application uses one base action and O(N k + k^2). CPU/CUDA support float32/float64; MPS requires float32. The declared spectral_device is CPU for MPS and the data device otherwise. Only the k-by-k Gram matrix and scalar roundoff bound go to that device for eigenvalue/rank validation and singular Hermitian inversion. Its pseudoinverse returns to the data device through differentiable Torch transfers. Pixel arrays, base actions, contractions, Cholesky factors/solves and output stay on the data device. This explicit small-matrix policy does not enable global MPS CPU fallback or replace PSD/rank checks with a shifted positive-definiteness test.

_cholesky: Tensor | None = None

_pseudoinverse: Tensor | None = None

_weighted_basis = _weighted_columns(_flat_basis)

nuisance_rank = 0

rank_atol = float(rank_atol)

rank_rtol = (
    _n_modes * eps
    if rank_rtol is None
    else float(rank_rtol)
)

treatment = 'profiled_unknown'

apply(residual: Tensor) -> Tensor