perf(mmle): use matrix-vector EAP projection

docs/changelog.d/mmle-fallback-eap-matmul.md

@@ -0,0 +1,6 @@
+# MMLE fallback EAP matrix-vector projection
+
+## Changed
+
+- Replaced the NumPy fallback MMLE EAP broadcast-and-reduce expression with an equivalent dense matrix-vector product, avoiding the additional posterior-shaped temporary array and permitting optimized numerical-library dispatch when available.
+- Added an independently reconstructed missing-data parity fixture, a source-level allocation-path regression, and APA 7th doctoring while preserving the Rust primary backend and all statistical contracts.

docs/doctoring/mmle-fallback-eap-matmul.md

@@ -0,0 +1,40 @@
+# Doctoring record: fallback MMLE EAP matrix-vector projection
+
+## Decision
+
+The NumPy reference/fallback MMLE implementation computes each respondent's expected-a-posteriori ability with the dense matrix-vector product
+
+\[
+\widehat{\theta}_p = \sum_q \pi_{pq} x_q,
+\]
+
+implemented as `posterior @ nodes`.
+
+This replaces an element-wise broadcast followed by an axis reduction. The statistical quantity, quadrature nodes, posterior weights, missing-data handling, estimator initialization, M-step, stopping rule, and returned transport remain unchanged.
+
+## Architectural boundary
+
+The repository's resolved production backend remains Rust. This change optimizes the explicitly retained NumPy reference/fallback path; it does not move production psychometric arithmetic from Rust into Python and does not introduce a second estimator contract.
+
+For a posterior matrix of shape `(n_persons, n_nodes)`, the former expression materialized an element-wise product with the same shape before reduction. NumPy documents `matmul` and the `@` operator as matrix-product operations and notes that optimized BLAS is used when possible. Exact performance depends on array shape, layout, linked numerical libraries, hardware, and runtime conditions, so this change makes no universal speedup claim.
+
+## Verification contract
+
+`tests/test_mmle_fallback_eap_matmul.py`:
+
+- reconstructs the one-iteration posterior independently for a realistic partially observed response matrix;
+- computes the former weighted-sum reference explicitly;
+- requires the public fallback result to agree at tight floating-point tolerance; and
+- pins the allocation-bounded matrix-vector source path so the broadcast temporary is not silently restored.
+
+The complete repository CI remains authoritative for Rust/PyO3 tests, fallback behavior, package acceptance, GPU no-skip evidence, fuzzing, security scans, and the production coverage/docstring gates.
+
+## Interpretation boundary
+
+Numerical parity of this projection does not establish parameter recovery, model fit, global optimality, construct validity, fairness, or operational readiness. Those claims require the repository's separate simulation, recovery, validation, and governance evidence.
+
+## References
+
+Bock, R. D., & Aitkin, M. (1981). Marginal maximum likelihood estimation of item parameters: Application of an EM algorithm. *Psychometrika, 46*(4), 443–459. https://doi.org/10.1007/BF02293801
+
+NumPy Developers. (2026). *numpy.matmul—NumPy v2.5 manual*. https://numpy.org/doc/stable/reference/generated/numpy.matmul.html

python/fast_mlsirm/estimators/mmle.py

@@ -184,7 +184,7 @@ def fit_mmle_2pl(
             break
 
     # ---- EAP ability for each person ----
-    theta = (posterior * nodes[None, :]).sum(axis=1)
+    theta = posterior @ nodes
 
     return {
         "a": a,

tests/test_mmle_fallback_eap_matmul.py

@@ -0,0 +1,79 @@
+"""Numerical and implementation contracts for fallback MMLE EAP projection."""
+
+from __future__ import annotations
+
+import inspect
+
+import numpy as np
+
+from fast_mlsirm.estimators.mmle import fit_mmle_2pl, gauss_hermite_nodes
+
+
+def test_one_iteration_eap_matches_independent_weighted_sum_reference() -> None:
+    """Dense matrix-vector projection must preserve the original EAP equation."""
+    y = np.array(
+        [
+            [1.0, 0.0, 1.0],
+            [0.0, 1.0, 0.0],
+            [1.0, 1.0, 0.0],
+            [0.0, 0.0, 1.0],
+        ],
+        dtype=np.float64,
+    )
+    observed = np.array(
+        [
+            [True, True, True],
+            [True, False, True],
+            [True, True, False],
+            [False, True, True],
+        ],
+        dtype=bool,
+    )
+    n_nodes = 9
+    seed = 7
+
+    y_filled = np.where(observed, y, 0.0)
+    obs_f = observed.astype(np.float64)
+    nodes, weights = gauss_hermite_nodes(n_nodes)
+    p_item = (y_filled * obs_f).sum(0) / np.clip(obs_f.sum(0), 1.0, None)
+    p_item = np.clip(p_item, 0.02, 0.98)
+    rng = np.random.default_rng(seed)
+    discrimination = np.ones(y.shape[1]) + 0.01 * rng.standard_normal(y.shape[1])
+    intercept = np.log(p_item / (1.0 - p_item))
+    logit = nodes[:, None] * discrimination[None, :] + intercept[None, :]
+    log_p1 = -np.logaddexp(0.0, -logit)
+    log_p0 = -np.logaddexp(0.0, logit)
+    log_joint = (
+        (y_filled * obs_f) @ log_p1.T
+        + ((1.0 - y_filled) * obs_f) @ log_p0.T
+        + np.log(weights)[None, :]
+    )
+    maximum = log_joint.max(axis=1, keepdims=True)
+    stabilized = np.exp(log_joint - maximum)
+    posterior = stabilized / stabilized.sum(axis=1, keepdims=True)
+    expected_theta = (posterior * nodes[None, :]).sum(axis=1)
+
+    result = fit_mmle_2pl(
+        y,
+        observed,
+        n_nodes=n_nodes,
+        max_iter=1,
+        seed=seed,
+    )
+
+    assert result["status"] == "max_iter_reached"
+    assert result["n_iter"] == 1
+    np.testing.assert_allclose(
+        np.asarray(result["theta"], dtype=np.float64),
+        expected_theta,
+        rtol=1e-13,
+        atol=1e-13,
+    )
+
+
+def test_eap_projection_retains_the_allocation_bounded_matmul_path() -> None:
+    """The fallback must not restore the posterior-by-node temporary array."""
+    source = inspect.getsource(fit_mmle_2pl)
+
+    assert "theta = posterior @ nodes" in source
+    assert "posterior * nodes[None, :]" not in source