bayesmith.diagnose.identifiability
Source-derived signatures and contracts. See the reference index for coverage and conventions.
Module contract
Cross-block identifiability: the rank test a per-block guard cannot perform.
Every convergence guard in :mod:`bayesmith.exact` is computed **from one
block**. That is not an oversight, it is arithmetic: a residual ``|M x - b|``
and a condition bound ``kappa(A^T N^-1 A + S^-1)`` are both properties of the
operator for the block being solved, so neither can see a degeneracy whose
two halves live in *different* blocks. The linearity check cannot see it
either, and it is right not to: each conditional of a bilinear model
genuinely is affine.
The failure that follows is silent and large, and it is the measured case
this module was built against (ported from
``rheplicant.inference.identifiability``, re-measured on this package's own
graph vocabulary): an alternating solve over ``gain x T_ant`` with a free
antenna temperature per (time, frequency) cell reports a clean per-block
kappa and residual while sitting thousands of kelvin from the truth. Nothing
in the sweep is wrong; the *partition* is, and no per-block number is
entitled to say so.
What can say so is the rank of the Jacobian of the prediction with respect
to **all** the selected latents at once. Re-measured on the graph port of
that fixture::
free-per-cell T_ant, tone ON (5000 K) n_par=72 rank=64 nullity=8
free-per-cell T_ant, tone OFF n_par=72 rank=64 nullity=8
(3,3)-basis T_coeff, tone ON (5000 K) n_par=17 rank=17 nullity=0
(3,3)-basis T_coeff, tone OFF n_par=17 rank=16 nullity=1
Read that as: a known calibration tone buys **exactly nothing** against a
free-per-cell antenna temperature -- the free cell at the tone's channel
absorbs the gain sample by sample, so the nullity stays at ``n_time`` either
way -- and **everything** against a frequency-smooth one, where a delta at
one channel is not in the span of three smooth basis functions and cannot be
reabsorbed.
Three things about the method are not decoration.
**The Jacobian's columns are normalised.** A latent whose natural scale is
1e3 and one whose scale is 1e-3 produce columns differing by 1e6 in norm, and
a rank verdict taken on those reports the choice of units rather than the
identifiability of the model. Column normalisation measures each parameter in
units of its own effect on the prediction, which is the only scale-free
question there is to ask.
**It refuses single precision rather than forcing double.** rheplicant runs
this diagnostic inside a process-global x64 switch of its own; this package's
rule is that ``src/`` never touches ``jax.config``, so the caller opens
``with jax.enable_x64(True):`` -- graph construction included -- and a
float32 Jacobian is refused by name (see
:func:`~bayesmith.diagnose.local.refuse_single_precision`). The refusal is
load-bearing, not ceremony, and the boundary was re-measured under THIS
regime rather than ported: in float64 the motivating model's null direction
sits at **7.479e-17** of the largest singular value; computed in float32 the
same direction surfaces at **3.117e-8** -- above the default tolerance -- and
the degenerate model is reported as fully identified.
**The result is named.** An anonymous index into a flattened vector tells a
user they have a problem and nothing about which;
:meth:`IdentifiabilityReport.direction` hands back ``{"gain": ..., "t_ant":
...}``, shaped like the latents, so "the degenerate direction is *this*
combination" is something you can read and act on.
One thing about it is a limit rather than a feature. **Cost:** a dense
Jacobian and a dense SVD, ``O(n_data * n_par * min(...))`` time and
``n_data * n_par`` float64 words of memory. This is a design-time diagnostic
for tens to a few thousand parameters -- the size a Gibbs partition is
*chosen* at -- not something to run inside a sweep over a 10^6-coefficient
block. For one block's conditioning without forming anything, use
:func:`bayesmith.exact.solve.condition_bound`.DEFAULT_RANK_RTOL
Defined in bayesmith.diagnose.identifiability, line 154.
DEFAULT_RANK_RTOL: float = 1e-08IdentifiabilityReport
Defined in bayesmith.diagnose.identifiability, line 158.
class IdentifiabilityReport
names: tuple[str, ...]
shapes: tuple[tuple[int, ...], ...]
spans: tuple[tuple[int, int], ...]
n_par: int
n_data: int
rank: int
nullity: int
singular_values: np.ndarray
null_space: np.ndarray
jacobian: np.ndarray
column_norms: np.ndarray
rtol: float
threshold: floatWhat the joint Jacobian's rank says about a set of latents.
Deliberately a plain frozen dataclass rather than an ``eqx.Module``, for
the same reason :class:`~bayesmith.exact.block.LinearBlock` is: this is a
derived linear-algebra verdict, not a differentiable model. It holds
**numpy** arrays, not JAX ones -- a float64 JAX array that escapes the
caller's x64 context truncates, with a warning, the moment a
default-precision caller touches it, throwing away exactly the precision
the diagnostic exists to obtain. And ``rank``/``nullity`` are Python
ints, which is why this function cannot be jitted: a rank is a decision,
and a traced decision is one you cannot branch on.
Attributes:
names: the latents analysed, in the order the caller asked for.
shapes: their shapes, in the same order.
spans: ``(start, stop)`` of each latent within the flat parameter
vector, in the same order.
n_par: total number of real parameters.
n_data: size of the flattened prediction over every observed node.
``nullity`` can never be below ``n_par - n_data``: more
parameters than data points is a null space by counting alone.
rank: number of singular values of the COLUMN-NORMALISED Jacobian
strictly above :attr:`threshold`.
nullity: ``n_par - rank`` -- the dimension of the space of parameter
perturbations the prediction is blind to.
singular_values: ``(n_par,)`` descending. When ``n_data < n_par`` the
SVD returns only ``n_data`` values and the rest are exact zeros;
they are included rather than dropped, so ``rank`` is always the
count of entries above the threshold and never needs a caveat.
null_space: ``(nullity, n_par)`` orthonormal rows, in the
column-normalised coordinates the rank verdict is taken in. Use
:meth:`direction` for raw latent coordinates.
jacobian: ``(n_data, n_par)`` column-normalised, as analysed. Rows
are the observed nodes in sorted name order, each flattened --
:func:`~bayesmith.exact.fisher.dense_operator`'s own layout.
column_norms: the SAFE norms the normalisation divided by -- the real
norm where a column is live, 1.0 where it is exactly zero.
rtol: the relative tolerance used.
threshold: ``rtol * singular_values[0]``, the absolute cutoff.IdentifiabilityReport.weakest_identified
weakest_identified(self) -> float``s[rank-1] / s[0]`` -- how well the worst identified direction is seen. The headline number: how much less the data says about the direction it constrains least than about the one it constrains most. ``0.0`` when nothing at all is identified, which is the only case where the ratio has no meaning.
IdentifiabilityReport.direction
direction(self, index: int) -> dict[str, np.ndarray]One null direction in RAW latent coordinates, split by name. Add a small multiple of this to the latents and the prediction does not move, to first order -- that is the whole content of the report, and it is the form a caller acts in. Note that it is NOT the raw SVD row: the SVD is taken of the column-normalised Jacobian, so a null vector there has to be divided by the column norms again to become a perturbation of the parameters themselves. Returned with unit 2-norm over the flat vector; the scale is arbitrary, only the direction means anything. For per-latent weights comparable across quantities in different units, use :meth:`participation` instead -- in raw kelvin and dimensionless gain, a null direction's two halves are not comparable numbers.
IdentifiabilityReport.participation
participation(self, index: int) -> dict[str, float]Fraction of a null direction carried by each latent, summing to 1. Measured in the COLUMN-NORMALISED coordinates, not raw ones: a 3000 K antenna temperature and a gain near 1 cannot be compared in their own units, and a raw-unit share would report which quantity is numerically larger rather than which one the degeneracy involves. In normalised coordinates the bilinear ``gain x T_ant`` degeneracy comes out at 0.50/0.50, which is the true statement about it.
identifiability
Defined in bayesmith.diagnose.identifiability, line 306.
identifiability(graph: Graph, *, names: Sequence[str] | str | None=None, at: dict[str, jax.Array] | None=None, rtol: float=DEFAULT_RANK_RTOL) -> IdentifiabilityReportRank of the joint Jacobian: what the data cannot tell apart.
The diagnostic that sees ACROSS blocks. See the module docstring for why
no per-block guard can, and for the measured case that motivates it.
Args:
graph: the model.
names: which latents to differentiate with respect to -- a sequence,
or a bare string for one. ``None`` (the default) means all of
them, in declaration order. A subset asks the **conditional**
question a Gibbs block faces -- "is this block identified, with
the others held fixed?" -- and the answer is routinely *yes* for
every block of a partition whose joint model is degenerate. That
is the whole reason this function takes the joint by default.
at: values for the latents NOT selected (and evaluation values for
those that are). Identifiability is a LOCAL property of a
nonlinear model, so a sweep has to ask it where the sampler
currently is; defaults to each latent's prior centre, which is
right exactly once.
rtol: singular values at or below ``rtol * s_max`` are called null.
See :data:`DEFAULT_RANK_RTOL` for how the default is chosen and
when to override it.
Returns:
An :class:`IdentifiabilityReport`. The three numbers to read first
are ``n_par``, ``nullity`` and, when the nullity is non-zero,
``report.participation(0)`` -- which names the latents the degenerate
direction mixes::
report = identifiability(graph)
if report.nullity:
print(report.nullity, "blind directions;", report.participation(0))
Raises:
GraphError: if ``names`` is empty, repeats a latent, or names one the
graph does not declare; if ``at`` names one, or a value does not
broadcast; if a selected latent is complex or non-floating; if
the graph has no observed node or an observed node's distribution
carries no ``loc``; or if the Jacobian comes back float32 -- run
the call (graph construction included) inside
``with jax.enable_x64(True):``.