"""
Local identifiability diagnostics for the BLIM.
The BLIM is not identifiable in general: structurally different parameter
vectors (β, η, π) can induce the *same* distribution over response
patterns. Stefanutti, Heller, Anselmi and Robusto (2012) characterize
local identifiability through the rank of the Jacobian of the prediction
map θ ↦ (P(R))_R. Full column rank is a sufficient regular local criterion;
a deficient rank at an isolated singular point is not its converse.
The unrestricted parameter count is ``2|Q| + |K| - 1``. Forward- and
backward-gradedness of the structure (Spoto, Stefanutti & Vidotto 2013)
are simple structural conditions that already imply trade-off dimensions
(η_q ↔ π for forward-graded items, β_q ↔ π for backward-graded items).
These diagnostics are the counterpart of ``pks::jacobian()``,
``pks::is.forward.graded()`` and ``pks::is.backward.graded()`` in R.
Because the rank at a single point only *lower-bounds* the generic rank
of the model, :func:`check_identifiability` evaluates the Jacobian at
several interior points (a user-suppliable point plus deterministic
random draws) and reports the maximum rank; see its docstring for the
one-sided logic of the resulting verdict.
References:
Stefanutti, L., Heller, J., Anselmi, P., & Robusto, E. (2012).
Assessing the local identifiability of probabilistic knowledge
structures. Behavior Research Methods, 44(4), 1197-1211.
Spoto, A., Stefanutti, L., & Vidotto, G. (2013). Considerations about
the identification of forward- and backward-graded knowledge
structures. Journal of Mathematical Psychology, 57(5), 249-254.
Heller, J. (2017). Identifiability in probabilistic knowledge
structures. Journal of Mathematical Psychology, 77, 46-57.
"""
from __future__ import annotations
import warnings
from collections.abc import Mapping
from dataclasses import dataclass
import numpy as np
from numpy.typing import NDArray
from knowledgespaces._limits import check_domain_size
from knowledgespaces.estimation.constraints import BLIMConstraints, _compile_constraints
from knowledgespaces.estimation.prediction import _item_vector, _log_conditional
from knowledgespaces.structures.knowledge_structure import KnowledgeStructure
[docs]
def is_forward_graded(structure: KnowledgeStructure, item: str) -> bool:
"""Whether the structure is forward-graded in ``item``.
A knowledge structure 𝒦 is forward-graded in item q iff
``{K ∪ {q} : K ∈ 𝒦} ⊆ 𝒦`` (Spoto, Stefanutti & Vidotto 2013).
Forward-gradedness in q implies a trade-off between ``eta[q]`` and the
state probabilities, hence local non-identifiability of the BLIM.
Mirrors ``pks::is.forward.graded()``.
"""
if item not in structure.domain:
raise ValueError(f"Item {item!r} is not in the domain.")
states = set(structure.states)
return all(frozenset(state | {item}) in states for state in states)
[docs]
def is_backward_graded(structure: KnowledgeStructure, item: str) -> bool:
"""Whether the structure is backward-graded in ``item``.
A knowledge structure 𝒦 is backward-graded in item q iff
``{K \\ {q} : K ∈ 𝒦} ⊆ 𝒦`` (Spoto, Stefanutti & Vidotto 2013).
Backward-gradedness in q implies a trade-off between ``beta[q]`` and
the state probabilities, hence local non-identifiability of the BLIM.
Mirrors ``pks::is.backward.graded()``.
"""
if item not in structure.domain:
raise ValueError(f"Item {item!r} is not in the domain.")
states = set(structure.states)
return all(frozenset(state - {item}) in states for state in states)
def _as_vector(
value: float | Mapping[str, float] | np.ndarray,
items: list[str],
name: str,
) -> np.ndarray:
"""Expand scalar / mapping / array parameter into an item-aligned array."""
if isinstance(value, Mapping):
missing = [q for q in items if q not in value]
if missing:
raise ValueError(f"{name} is missing items {missing}.")
extra = set(value) - set(items)
if extra:
raise ValueError(f"{name} has unknown items {sorted(extra)}.")
vec = np.array([float(value[q]) for q in items], dtype=np.float64)
elif isinstance(value, np.ndarray):
vec = value.astype(np.float64)
if vec.shape != (len(items),):
raise ValueError(f"{name} array has shape {vec.shape}, expected ({len(items)},).")
else:
vec = np.full(len(items), float(value), dtype=np.float64)
# The finiteness check must come first: NaN passes every comparison
# below (NaN < x is always False) and would otherwise surface much
# later as a cryptic LinAlgError from the SVD.
if not np.all(np.isfinite(vec)) or np.any(vec <= 0) or np.any(vec >= 1):
raise ValueError(f"All {name} values must be finite and in (0, 1) for the Jacobian.")
return vec
[docs]
def blim_jacobian(
structure: KnowledgeStructure,
*,
beta: float | Mapping[str, float] | np.ndarray = 0.1,
eta: float | Mapping[str, float] | np.ndarray = 0.1,
pi: Mapping[frozenset[str], float] | np.ndarray | None = None,
constraints: BLIMConstraints | None = None,
max_items: int | None = None,
max_memory_bytes: int = 8_000_000_000,
) -> np.ndarray:
"""Analytic Jacobian of the BLIM prediction map at (β, η, π).
The prediction map sends the parameter vector
``θ = (β_1..β_Q, η_1..η_Q, π_2..π_|K|)`` (the probability of the first
canonical state is ``1 - Σ`` of the others) to the vector of response
pattern probabilities ``(P(R))_R`` over all ``2^Q`` patterns. The BLIM
has regular local identifiability at θ when this Jacobian has full
column rank ``npar = 2|Q| + |K| - 1`` (Stefanutti, Heller, Anselmi &
Robusto 2012). A deficient rank at a singular point alone does not
establish generic non-identifiability.
Parameters
----------
structure : KnowledgeStructure
The knowledge structure.
beta, eta : float, Mapping[str, float], or np.ndarray
Evaluation point, scalar, per-item mapping, or array aligned with
the sorted domain. Without constraints, values must lie in (0, 1).
With constraints, boundary values in [0, 1] are accepted for
evaluating the restricted polynomial map. Default 0.1 for both.
pi : Mapping, np.ndarray, or None
State probabilities at the evaluation point (canonical state
order: by size, then lexicographically). None (default) uses the
uniform prior.
constraints : BLIMConstraints or None
Restrict the evaluation point and Jacobian to the declared model.
Fixed values override inputs; free equality groups use their mean.
Columns are free beta groups, free eta groups (ordered by first
item in the sorted domain), then free pi coordinates. Fixed pi
removes all pi columns. Boundary derivatives use products excluding
the differentiated item, avoiding division by zero.
max_items : int or None
Override for the ``2^Q`` pattern-enumeration guard.
max_memory_bytes : int
Hard cap on the estimated peak allocation, which scales with
``2^Q × |K|`` (several concurrent float64 arrays of that shape).
A :class:`MemoryError` is raised before any large array is
allocated; a :class:`ResourceWarning` is emitted above 1 GB.
Default 8 GB.
Returns
-------
np.ndarray
Jacobian of shape ``(2^Q, npar)``, rows in lexicographic pattern
order, columns ordered ``(beta_q)_q, (eta_q)_q, (pi_k)_{k>=2}``.
Its rank equals the rank of the ``(2^Q - 1)``-row version used by
``pks::jacobian()``: every column sums to zero across patterns
(Σ_R P(R) = 1), so dropping one row does not change the rank.
Note that the rank at a single point only lower-bounds the
generic rank of the model (matrix rank is lower semicontinuous);
:func:`check_identifiability` therefore evaluates several points.
Raises
------
ValueError
On invalid parameters.
DomainTooLargeError
If ``|Q|`` exceeds the enumeration guard.
MemoryError
If the estimated peak allocation exceeds ``max_memory_bytes``.
"""
items = sorted(structure.domain)
states = sorted(structure.states, key=lambda s: (len(s), sorted(s)))
n_items = len(items)
n_states = len(states)
if max_items is not None:
check_domain_size(
n_items,
max_threshold=max_items,
context="blim_jacobian (enumerates all 2^Q response patterns)",
)
else:
check_domain_size(
n_items,
context="blim_jacobian (enumerates all 2^Q response patterns)",
)
compiled = None if constraints is None else _compile_constraints(constraints, items, states)
resolver = _as_vector if compiled is None else _item_vector
beta_vec = resolver(beta, items, "beta")
eta_vec = resolver(eta, items, "eta")
if compiled is not None:
compiled.beta.project(beta_vec, clip_free=False)
compiled.eta.project(eta_vec, clip_free=False)
if pi is None:
pi_vec: NDArray[np.float64] = np.full(n_states, 1.0 / n_states)
elif isinstance(pi, Mapping):
missing = [s for s in states if s not in pi]
if missing:
raise ValueError(f"pi is missing states {sorted(missing, key=sorted)}.")
if set(pi) - set(states):
raise ValueError("pi has unknown states.")
pi_vec = np.array([float(pi[s]) for s in states], dtype=np.float64)
else:
pi_vec = np.asarray(pi, dtype=np.float64)
if pi_vec.shape != (n_states,):
raise ValueError(f"pi array has shape {pi_vec.shape}, expected ({n_states},).")
if not np.all(np.isfinite(pi_vec)) or np.any(pi_vec < 0) or abs(pi_vec.sum() - 1.0) > 1e-8:
raise ValueError("pi must be finite, non-negative and sum to 1.")
if compiled is not None and compiled.pi is not None:
pi_vec = compiled.pi.copy()
# Memory preflight. Peak concurrent allocations, all rows = 2^Q:
# ~7 (patterns x states) float64 arrays (L, l_q with its two
# outer-product temporaries, ratio, expression temporaries and margin),
# the Jacobian J and its constrained projection, (patterns x npar) with
# npar = 2|Q| + |K| - 1,
# ~3 (patterns x items) arrays (the pattern matrix R materialized
# as int64 then float64, plus the 1 - R temporary).
# The (patterns x items) and (patterns x npar) terms dominate for
# sparse structures (small |K|), the (patterns x states) terms for
# powerset-like structures — the multipliers are calibrated as an
# conservative estimate on tracemalloc peaks across both shapes (see
# tests/test_identifiability.py), including Python 3.14 allocation
# overhead. This preflight is not an operating-system memory guarantee.
n_patterns = 2**n_items
npar_preflight = 2 * n_items + n_states - 1
estimated_bytes = n_patterns * (7 * n_states + 2 * npar_preflight + 3 * n_items) * 8
if estimated_bytes > max_memory_bytes:
raise MemoryError(
f"blim_jacobian would allocate ~{estimated_bytes / 1e9:.2f} GB "
f"({n_patterns} patterns x [7 x {n_states} states + "
f"2 x {npar_preflight} parameters + 3 x {n_items} items] x 8 bytes). "
f"Exceeds max_memory_bytes={max_memory_bytes / 1e9:.2f} GB. "
f"Use a sparser structure or pass max_memory_bytes=... to "
f"override."
)
if estimated_bytes > 1_000_000_000:
warnings.warn(
f"blim_jacobian will allocate ~{estimated_bytes / 1e9:.2f} GB "
f"({n_patterns} patterns x {n_states} states).",
category=ResourceWarning,
stacklevel=2,
)
# State membership S[k, q] and all 2^Q response patterns R[r, q]
item_idx = {q: i for i, q in enumerate(items)}
S = np.zeros((n_states, n_items), dtype=np.float64)
for k, state in enumerate(states):
for q in state:
S[k, item_idx[q]] = 1.0
R = ((np.arange(n_patterns)[:, None] >> np.arange(n_items - 1, -1, -1)) & 1).astype(np.float64)
# Response likelihood L[r, k] = P(R_r | K_k), computed in log space via
# matrix products so that no (patterns x states x items) 3D tensor is
# ever materialized. Note the remaining 2D arrays are still 2^Q x |K|
# — hence the memory preflight above.
# P(R_q = 1 | K) = S*(1-beta) + (1-S)*eta, strictly inside (0, 1) here.
p_correct = S * (1 - beta_vec) + (1 - S) * eta_vec # (n_states, n_items)
boundary = np.any(p_correct == 0) or np.any(p_correct == 1)
if boundary:
L = np.exp(_log_conditional(R, p_correct))
else:
log_pc = np.log(p_correct)
log_pic = np.log(1 - p_correct)
L = np.exp(R @ log_pc.T + (1 - R) @ log_pic.T)
npar = 2 * n_items + n_states - 1
J = np.zeros((n_patterns, npar))
# ∂P(R)/∂beta_q = Σ_K π_K [q ∈ K] (L / l_q) * (∂l_q/∂beta_q)
# with ∂l_q/∂beta_q = -1 if R_q = 1 else +1 (only for q ∈ K);
# ∂P(R)/∂eta_q analogous with +1 if R_q = 1 else -1 (only for q ∉ K).
for qi in range(n_items):
# l_q[r, k] = R[r,q]*p_correct[k,q] + (1-R[r,q])*(1-p_correct[k,q])
l_q = np.outer(R[:, qi], p_correct[:, qi]) + np.outer(1 - R[:, qi], 1 - p_correct[:, qi])
if boundary:
others = np.arange(n_items) != qi
ratio = np.exp(_log_conditional(R[:, others], p_correct[:, others]))
else:
ratio = L / l_q # interior point
sign = 1.0 - 2.0 * R[:, qi] # +1 if incorrect, -1 if correct
J[:, qi] = (ratio * (S[:, qi] * pi_vec)).sum(axis=1) * sign
J[:, n_items + qi] = (ratio * ((1 - S[:, qi]) * pi_vec)).sum(axis=1) * (-sign)
# ∂P(R)/∂pi_k = L(R | K_k) - L(R | K_1) for the free components k >= 2
# (pi_1 = 1 - Σ_{k>=2} pi_k).
J[:, 2 * n_items :] = L[:, 1:] - L[:, :1]
return J if compiled is None else compiled.restrict_jacobian(J, n_items)
[docs]
@dataclass(frozen=True)
class IdentifiabilityReport:
"""Result of a BLIM local-identifiability check.
Attributes
----------
rank : int
Numerical rank of the Jacobian of the prediction map, maximized
over the evaluated interior points (the user-supplied point plus
the random points drawn by :func:`check_identifiability`). Since
the rank at any point lower-bounds the generic rank, this is the
best available lower bound on the generic rank.
npar : int
Number of free parameters after fixed/equality constraints;
``2|Q| + |K| - 1`` for the unrestricted model.
locally_identifiable : bool
True iff ``rank == npar``. Full rank at *any* interior point
establishes generic local identifiability, up to the numerical
rank decision (see ``smallest_singular_value`` and
``rank_tolerance`` for its margin); generic statements are
silent on null sets such as the boundary, so for boundary
solutions re-evaluate at the fitted point. When False, no evaluated
point had full rank — strong evidence (not proof) that ``npar -
rank`` independent trade-off dimensions exist generically:
parameter estimates are not unique and only likelihood-based
quantities (log-likelihood, G², AIC, BIC, predicted pattern
probabilities) are comparable across implementations or fits.
rank_deficiency : int
``npar - rank``.
forward_graded : tuple[str, ...]
Items in which the structure is forward-graded (η ↔ π trade-off).
backward_graded : tuple[str, ...]
Items in which the structure is backward-graded (β ↔ π trade-off).
smallest_singular_value : float
The ``npar``-th singular value of the Jacobian at the point that
achieved the reported rank. When the rank is full, this is the
smallest counted singular value (compare with rank_tolerance);
when deficient, it is the smallest column-direction singular
value, not necessarily the first discarded one. It is zero if
there are fewer singular values than parameters.
rank_tolerance : float
The singular-value cutoff used for the rank decision at that
point (NumPy convention: ``max(J.shape) * eps * sigma_max``).
The numerical rank is the number of singular values above it.
"""
rank: int
npar: int
locally_identifiable: bool
rank_deficiency: int
forward_graded: tuple[str, ...]
backward_graded: tuple[str, ...]
smallest_singular_value: float = 0.0
rank_tolerance: float = 0.0
[docs]
def check_identifiability(
structure: KnowledgeStructure,
*,
beta: float | Mapping[str, float] | np.ndarray = 0.1,
eta: float | Mapping[str, float] | np.ndarray = 0.1,
pi: Mapping[frozenset[str], float] | np.ndarray | None = None,
constraints: BLIMConstraints | None = None,
n_points: int = 5,
seed: int = 0,
max_items: int | None = None,
max_memory_bytes: int = 8_000_000_000,
) -> IdentifiabilityReport:
"""Check local identifiability of the BLIM on a knowledge structure.
Computes the rank of the analytic Jacobian of the prediction map at
the given evaluation point (default: β = η = 0.1, uniform π) and, if
that rank is deficient, at up to ``n_points`` additional random
interior points (β_q, η_q ~ U(0.05, 0.3) independently per item,
π ~ Dirichlet(1, …, 1)) drawn from a fixed-seed generator so results
are deterministic. The reported rank is the *maximum* over the
evaluated points, and the model is declared locally identifiable iff
that rank equals ``npar = 2|Q| + |K| - 1`` (Stefanutti, Heller,
Anselmi & Robusto 2012).
Parameters
----------
structure : KnowledgeStructure
The knowledge structure.
beta, eta : float, Mapping[str, float], or np.ndarray
User-suppliable evaluation point, e.g. fitted estimates:
``check_identifiability(s, beta=est.beta_dict(), eta=est.eta_dict(),
pi=est.pi_dict(), n_points=0)``. Dictionaries preserve item
alignment even if estimation used an unsorted column order;
``n_points=0`` evaluates the fitted point alone. Default 0.1 for both.
pi : Mapping, np.ndarray, or None
State probabilities at the evaluation point. None (default) uses
the uniform prior.
n_points : int
Number of additional random interior points at which the Jacobian
is evaluated when the user-supplied point is rank deficient
(evaluation stops early once full rank is found). ``0`` restores
the single-point check. Default 5.
seed : int
Seed of the internal random generator that draws the extra
points; fixed by default so that repeated calls give identical
reports. Default 0.
constraints : BLIMConstraints or None
Restrict the evaluation point and Jacobian to the declared model.
Fixed values override inputs; free equality groups use their mean.
Columns are free beta groups, free eta groups (ordered by first
item in the sorted domain), then free pi coordinates. Fixed pi
removes all pi columns. Boundary derivatives use products excluding
the differentiated item, avoiding division by zero.
max_items : int or None
Override for the ``2^Q`` pattern-enumeration guard.
max_memory_bytes : int
Hard cap on the estimated peak allocation per Jacobian; see
:func:`blim_jacobian`. Default 8 GB.
Notes
-----
The rank of the Jacobian is constant — equal to the *generic rank* —
on an open dense subset of the parameter space, and can only *drop*
on the complementary null set (matrix rank is lower semicontinuous).
The mathematical criterion is one-sided: exact full rank at **any**
interior point establishes generic local identifiability. This
implementation uses a floating-point singular-value threshold, so
inspect the reported margin and tolerance. A deficient rank at
all evaluated points is strong evidence, not proof, of generic rank
deficiency. Evaluating only one point is not enough: symmetric
structures can place the symmetric default point β = η = 0.1 with
uniform π inside the exceptional null set (e.g. the knowledge space
{∅, {a, b}, {c, d}, Q} has generic rank 11 = npar but rank 7 at the
default point), which is why the random points are evaluated.
Examples
--------
>>> from knowledgespaces import space_from_prerequisites
>>> s = space_from_prerequisites(["a", "b"], [("a", "b")])
>>> report = check_identifiability(s)
>>> report.npar
6
"""
if n_points < 0:
raise ValueError(f"n_points must be non-negative, got {n_points}.")
J = blim_jacobian(
structure,
beta=beta,
eta=eta,
pi=pi,
constraints=constraints,
max_items=max_items,
max_memory_bytes=max_memory_bytes,
)
npar = J.shape[1]
def _rank_details(jac: np.ndarray) -> tuple[int, float, float]:
"""Numerical rank with its SVD details (NumPy matrix_rank convention)."""
if jac.shape[1] == 0:
return 0, 0.0, 0.0
s = np.linalg.svd(jac, compute_uv=False)
tol = float(s.max() * max(jac.shape) * np.finfo(jac.dtype).eps)
r = int(np.count_nonzero(s > tol))
sigma_npar = float(s[npar - 1]) if s.size >= npar else 0.0
return r, sigma_npar, tol
rank, sigma_npar, tol = _rank_details(J)
if rank < npar and n_points > 0:
# The user-supplied point may lie in the null set where the rank
# drops below the generic rank; probe random interior points and
# keep the maximum (a full-rank point settles the question).
rng = np.random.default_rng(seed)
n_items = len(structure.domain)
n_states = len(structure.states)
for _ in range(n_points):
beta_r = rng.uniform(0.05, 0.3, size=n_items)
eta_r = rng.uniform(0.05, 0.3, size=n_items)
pi_r = rng.dirichlet(np.ones(n_states))
J_r = blim_jacobian(
structure,
beta=beta_r,
eta=eta_r,
pi=pi_r,
constraints=constraints,
max_items=max_items,
max_memory_bytes=max_memory_bytes,
)
rank_r, sigma_r, tol_r = _rank_details(J_r)
if rank_r > rank:
rank, sigma_npar, tol = rank_r, sigma_r, tol_r
if rank == npar:
break
items = sorted(structure.domain)
forward = tuple(q for q in items if is_forward_graded(structure, q))
backward = tuple(q for q in items if is_backward_graded(structure, q))
return IdentifiabilityReport(
rank=rank,
npar=npar,
locally_identifiable=(rank == npar),
rank_deficiency=npar - rank,
forward_graded=forward,
backward_graded=backward,
smallest_singular_value=sigma_npar,
rank_tolerance=tol,
)