diff --git a/CHANGELOG.md b/CHANGELOG.md
index d08845599..c1261ff91 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,7 @@
 
 ### Fixed
 
+- Public multigroup M2 fails closed without the compiled Rust core and delegates target/null projected quadratic forms to native `projected_m2`.
 - Public conditional-Rasch M2 fails closed without the compiled Rust core and delegates every result field to the native `m2_cmle_rasch_stat` entrypoint.
 - Reject ambiguous LLM-judge JSON with duplicate keys or unexpected top-level fields; require the exact mode-specific schema including advisory `accepted`.
 - Require compiled Rust ownership for public `s_x2` and `person_fit`, including prior-mean S-X² dispatch, with fail-closed errors when the core is missing.
diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs
index 03698d32e..352ba8ad5 100644
--- a/crates/fast-mlsirm-py/src/lib.rs
+++ b/crates/fast-mlsirm-py/src/lib.rs
@@ -16,7 +16,10 @@ use mlsirm_core::fitstats::{
     benjamini_hochberg as core_benjamini_hochberg, chi2_sf as core_chi2_sf,
     infit_outfit as core_infit_outfit, leniency_residuals as core_leniency_residuals,
     m2_cmle_rasch as core_m2_cmle_rasch, m2_rmsea2 as core_m2, person_fit as core_person_fit,
-    poly_local_dependence as core_poly_ld, poly_m2 as core_poly_m2, s_x2 as core_s_x2, SX2Config,
+    poly_local_dependence as core_poly_ld, poly_m2 as core_poly_m2,
+    projected_m2 as core_projected_m2,
+    projected_m2_workspace_elements as core_projected_m2_workspace_elements,
+    s_x2 as core_s_x2, SX2Config, PROJECTED_M2_MAX_WORKSPACE_ELEMENTS,
 };
 use mlsirm_core::linking::{
     irt_link as core_irt_link, link_fixed_item_parameters as core_link_fixed_item_parameters,
@@ -7034,6 +7037,69 @@ fn m2_stat(
 }
 
 
+
+/// Projected M2 quadratic form ownership entrypoint (dense residual / Delta / Xi).
+#[pyfunction]
+fn projected_m2(
+    residual: PyReadonlyArray1<'_, f64>,
+    delta: PyReadonlyArray2<'_, f64>,
+    xi: PyReadonlyArray2<'_, f64>,
+    n: f64,
+) -> PyResult<f64> {
+    let residual_arr = residual.as_array();
+    let delta_arr = delta.as_array();
+    let xi_arr = xi.as_array();
+    let s = residual_arr.len();
+    if delta_arr.shape()[0] != s {
+        return Err(PyValueError::new_err(
+            "delta rows must match residual length",
+        ));
+    }
+    if xi_arr.shape() != [s, s] {
+        return Err(PyValueError::new_err("xi must be residual_len x residual_len"));
+    }
+    let p = delta_arr.shape()[1];
+    let workspace_elements = core_projected_m2_workspace_elements(s, p)
+        .map_err(PyValueError::new_err)?;
+    if workspace_elements > PROJECTED_M2_MAX_WORKSPACE_ELEMENTS {
+        return Err(PyValueError::new_err(
+            "projected M2 workspace exceeds supported element budget",
+        ));
+    }
+    if !n.is_finite() {
+        return Err(PyValueError::new_err("n must be finite"));
+    }
+    if residual_arr.iter().any(|value| !value.is_finite()) {
+        return Err(PyValueError::new_err("residual values must be finite"));
+    }
+    if delta_arr.iter().any(|value| !value.is_finite()) {
+        return Err(PyValueError::new_err("delta values must be finite"));
+    }
+    if xi_arr.iter().any(|value| !value.is_finite()) {
+        return Err(PyValueError::new_err("xi values must be finite"));
+    }
+    let residual = residual.as_slice()?;
+    let delta_elements = s
+        .checked_mul(p)
+        .ok_or_else(|| PyValueError::new_err("projected M2 dimensions overflow"))?;
+    let xi_elements = s
+        .checked_mul(s)
+        .ok_or_else(|| PyValueError::new_err("projected M2 dimensions overflow"))?;
+    let mut delta_flat = vec![0.0_f64; delta_elements];
+    for row in 0..s {
+        for col in 0..p {
+            delta_flat[row * p + col] = delta_arr[[row, col]];
+        }
+    }
+    let mut xi_flat = vec![0.0_f64; xi_elements];
+    for row in 0..s {
+        for col in 0..s {
+            xi_flat[row * s + col] = xi_arr[[row, col]];
+        }
+    }
+    core_projected_m2(residual, &delta_flat, xi_flat, s, p, n).map_err(PyValueError::new_err)
+}
+
 /// Conditional-Rasch M2 (CMLE ownership path).
 #[pyfunction]
 #[pyo3(signature = (y, observed, n_persons, item_easiness))]
@@ -9169,6 +9235,7 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(leniency_residuals_stat, m)?)?;
     m.add_function(wrap_pyfunction!(m2_stat, m)?)?;
     m.add_function(wrap_pyfunction!(m2_cmle_rasch_stat, m)?)?;
+    m.add_function(wrap_pyfunction!(projected_m2, m)?)?;
     m.add_function(wrap_pyfunction!(poly_m2, m)?)?;
     m.add_function(wrap_pyfunction!(poly_local_dependence, m)?)?;
     m.add_function(wrap_pyfunction!(poly_dif, m)?)?;
diff --git a/crates/mlsirm-core/src/fitstats.rs b/crates/mlsirm-core/src/fitstats.rs
index c387705e7..d93eb0b3c 100644
--- a/crates/mlsirm-core/src/fitstats.rs
+++ b/crates/mlsirm-core/src/fitstats.rs
@@ -1916,8 +1916,48 @@ fn chol_solve(l: &[f64], n: usize, b: &[f64]) -> Vec<f64> {
     x
 }
 
+/// Maximum supported dense workspace for projected M2, in `f64` elements.
+pub const PROJECTED_M2_MAX_WORKSPACE_ELEMENTS: usize = 4_000_000;
+
+/// Return the peak dense element budget required by projected M2.
+pub fn projected_m2_workspace_elements(
+    n_moments: usize,
+    n_parameters: usize,
+) -> Result<usize, String> {
+    let xi_elements = n_moments
+        .checked_mul(n_moments)
+        .ok_or_else(|| "projected M2 dimensions overflow".to_string())?;
+    let delta_elements = n_moments
+        .checked_mul(n_parameters)
+        .ok_or_else(|| "projected M2 dimensions overflow".to_string())?;
+    let parameter_elements = n_parameters
+        .checked_mul(n_parameters)
+        .ok_or_else(|| "projected M2 dimensions overflow".to_string())?;
+    let two_moments = n_moments
+        .checked_mul(2)
+        .ok_or_else(|| "projected M2 dimensions overflow".to_string())?;
+    let two_parameters = n_parameters
+        .checked_mul(2)
+        .ok_or_else(|| "projected M2 dimensions overflow".to_string())?;
+
+    [
+        xi_elements,
+        delta_elements,
+        delta_elements,
+        parameter_elements,
+        two_moments,
+        two_parameters,
+    ]
+    .into_iter()
+    .try_fold(0_usize, |total, elements| {
+        total
+            .checked_add(elements)
+            .ok_or_else(|| "projected M2 dimensions overflow".to_string())
+    })
+}
+
 /// Evaluate the M2 projected quadratic form without forming either inverse.
-fn projected_m2(
+pub fn projected_m2(
     e: &[f64],
     delta: &[f64],
     mut xi: Vec<f64>,
@@ -1927,6 +1967,39 @@ fn projected_m2(
 ) -> Result<f64, String> {
     let s = n_moments;
     let p = n_parameters;
+    let delta_elements = s
+        .checked_mul(p)
+        .ok_or_else(|| "projected M2 dimensions overflow".to_string())?;
+    let xi_elements = s
+        .checked_mul(s)
+        .ok_or_else(|| "projected M2 dimensions overflow".to_string())?;
+    let workspace_elements = projected_m2_workspace_elements(s, p)?;
+    if workspace_elements > PROJECTED_M2_MAX_WORKSPACE_ELEMENTS {
+        return Err(
+            "projected M2 workspace exceeds supported element budget".to_string(),
+        );
+    }
+    if e.len() != s {
+        return Err("residual length must match n_moments".to_string());
+    }
+    if delta.len() != delta_elements {
+        return Err("delta length must equal n_moments * n_parameters".to_string());
+    }
+    if xi.len() != xi_elements {
+        return Err("xi length must equal n_moments * n_moments".to_string());
+    }
+    if !n.is_finite() {
+        return Err("n must be finite".to_string());
+    }
+    if e.iter().any(|value| !value.is_finite()) {
+        return Err("residual values must be finite".to_string());
+    }
+    if delta.iter().any(|value| !value.is_finite()) {
+        return Err("delta values must be finite".to_string());
+    }
+    if xi.iter().any(|value| !value.is_finite()) {
+        return Err("xi values must be finite".to_string());
+    }
     cholesky_lower(&mut xi, s)?;
     let u = chol_solve(&xi, s, e); // Xi^-1 e
     let mut w = vec![0.0_f64; s * p]; // Xi^-1 Delta
diff --git a/crates/mlsirm-core/tests/projected_m2_validation.rs b/crates/mlsirm-core/tests/projected_m2_validation.rs
new file mode 100644
index 000000000..62d4b2c3d
--- /dev/null
+++ b/crates/mlsirm-core/tests/projected_m2_validation.rs
@@ -0,0 +1,65 @@
+//! Validation contracts for the public projected-M2 Rust boundary.
+
+use mlsirm_core::fitstats::projected_m2;
+
+#[test]
+fn projected_m2_rejects_short_residual_without_panicking() {
+    let result = projected_m2(
+        &[0.25],
+        &[1.0, 0.5],
+        vec![1.0, 0.0, 0.0, 1.0],
+        2,
+        1,
+        100.0,
+    );
+
+    assert!(result.is_err(), "short residual must return Err");
+}
+
+#[test]
+fn projected_m2_rejects_short_delta_without_panicking() {
+    let result = projected_m2(
+        &[0.25, -0.25],
+        &[1.0],
+        vec![1.0, 0.0, 0.0, 1.0],
+        2,
+        1,
+        100.0,
+    );
+
+    assert!(result.is_err(), "short delta must return Err");
+}
+
+#[test]
+fn projected_m2_rejects_short_covariance_without_panicking() {
+    let result = projected_m2(
+        &[0.25, -0.25],
+        &[1.0, 0.5],
+        vec![1.0, 0.0, 0.0],
+        2,
+        1,
+        100.0,
+    );
+
+    assert!(result.is_err(), "short covariance must return Err");
+}
+
+#[test]
+fn projected_m2_rejects_dimension_overflow_without_allocating() {
+    let result = projected_m2(&[], &[], Vec::new(), usize::MAX, 2, 1.0);
+
+    assert!(
+        matches!(result, Err(message) if message.contains("overflow")),
+        "overflowing dimensions must fail before indexing or allocation"
+    );
+}
+
+#[test]
+fn projected_m2_rejects_nonfinite_sample_size_before_arithmetic() {
+    let result = projected_m2(&[0.25], &[1.0], vec![1.0], 1, 1, f64::INFINITY);
+
+    assert!(
+        matches!(result, Err(message) if message.contains("n must be finite")),
+        "non-finite n must fail before numerical execution"
+    );
+}
diff --git a/docs/changelog.d/818-multigroup-m2-rust-projection.md b/docs/changelog.d/818-multigroup-m2-rust-projection.md
new file mode 100644
index 000000000..ee33d8366
--- /dev/null
+++ b/docs/changelog.d/818-multigroup-m2-rust-projection.md
@@ -0,0 +1,5 @@
+# Multigroup M2 Rust projection ownership
+
+## Fixed
+
+- Public `m2_multigroup` fails closed without the compiled Rust core and delegates target/null projected M2 quadratic forms to `projected_m2`.
diff --git a/python/fast_mlsirm/fitstats.py b/python/fast_mlsirm/fitstats.py
index 09079e1dc..6c0e17e08 100644
--- a/python/fast_mlsirm/fitstats.py
+++ b/python/fast_mlsirm/fitstats.py
@@ -2737,9 +2737,24 @@ def m2_multigroup(
         null_delta[rows, cols] = root_n * null_d
         null_xi_blocks.append(null_xi)
 
-    m2_value = _projected_m2_numpy(residual, delta, _block_diag(xi_blocks), 1.0)
-    null_m2 = _projected_m2_numpy(
-        null_residual, null_delta, _block_diag(null_xi_blocks), 1.0
+    core = _core_module()
+    if core is None or not hasattr(core, "projected_m2"):
+        raise RuntimeError("fit statistics require the compiled Rust core")
+    m2_value = float(
+        core.projected_m2(
+            np.ascontiguousarray(residual, dtype=np.float64),
+            np.ascontiguousarray(delta, dtype=np.float64),
+            np.ascontiguousarray(_block_diag(xi_blocks), dtype=np.float64),
+            1.0,
+        )
+    )
+    null_m2 = float(
+        core.projected_m2(
+            np.ascontiguousarray(null_residual, dtype=np.float64),
+            np.ascontiguousarray(null_delta, dtype=np.float64),
+            np.ascontiguousarray(_block_diag(null_xi_blocks), dtype=np.float64),
+            1.0,
+        )
     )
     df = float(n_groups * s - p)
     null_df = float(n_groups * s - n_groups * y0.shape[1])
diff --git a/tests/test_fitstats_multigroup_rust_ownership.py b/tests/test_fitstats_multigroup_rust_ownership.py
new file mode 100644
index 000000000..712410306
--- /dev/null
+++ b/tests/test_fitstats_multigroup_rust_ownership.py
@@ -0,0 +1,144 @@
+"""Fail-first Rust ownership contract for multigroup M2 (issue #627)."""
+
+from types import SimpleNamespace
+
+import numpy as np
+import pytest
+
+import fast_mlsirm.fitstats as fitstats
+
+
+def _multigroup_case():
+    """Return a small connected two-group MIRT case for ownership tests."""
+    rng = np.random.default_rng(627)
+    n_persons = 120
+    n_items = 6
+    group_id = np.repeat(np.arange(2, dtype=np.int64), n_persons // 2)
+    factor_id = np.zeros(n_items, dtype=np.int64)
+    group_mean = np.array([[0.0], [0.35]], dtype=np.float64)
+    group_sd = np.ones((2, 1), dtype=np.float64)
+    theta = rng.standard_normal(n_persons) + group_mean[group_id, 0]
+    intercept = np.linspace(-0.9, 0.9, n_items)
+    probability = 1.0 / (1.0 + np.exp(-(theta[:, None] + intercept[None, :])))
+    responses = (rng.random((n_persons, n_items)) < probability).astype(np.float64)
+    params = SimpleNamespace(
+        alpha=np.zeros(n_items),
+        b=intercept,
+        zeta=np.zeros((n_items, 1)),
+        tau=-30.0,
+    )
+    return responses, factor_id, params, group_id, group_mean, group_sd
+
+
+def _reject_reference_projection(*_args, **_kwargs):
+    """Fail if a public ownership test reaches the NumPy reference kernel."""
+    raise AssertionError("public multigroup M2 entered NumPy projection arithmetic")
+
+
+def _run_multigroup(case):
+    """Call the public multigroup M2 path with bounded quadrature."""
+    responses, factor_id, params, group_id, group_mean, group_sd = case
+    return fitstats.m2_multigroup(
+        responses,
+        factor_id,
+        params,
+        "MIRT",
+        group_id,
+        group_mean,
+        group_sd,
+        q_theta=7,
+        q_xi=7,
+    )
+
+
+def test_multigroup_m2_rejects_missing_core_before_reference_projection(monkeypatch):
+    """Public multigroup M2 must fail closed when the compiled core is absent."""
+    monkeypatch.setattr(fitstats, "_core_module", lambda: None)
+    monkeypatch.setattr(fitstats, "_projected_m2_numpy", _reject_reference_projection)
+
+    with pytest.raises(RuntimeError, match="fit statistics require the compiled Rust core"):
+        _run_multigroup(_multigroup_case())
+
+
+def test_multigroup_m2_requires_rust_projection_entrypoint(monkeypatch):
+    """A compiled core without the projection entrypoint is still non-passing."""
+    core_without_projection = SimpleNamespace(chi2_sf=lambda _x, _df: 0.5)
+    monkeypatch.setattr(fitstats, "_core_module", lambda: core_without_projection)
+    monkeypatch.setattr(fitstats, "_projected_m2_numpy", _reject_reference_projection)
+
+    with pytest.raises(RuntimeError, match="fit statistics require the compiled Rust core"):
+        _run_multigroup(_multigroup_case())
+
+
+def test_multigroup_m2_delegates_target_and_null_projection_to_rust(monkeypatch):
+    """Both public multigroup projection calls must cross the Rust boundary."""
+    calls = []
+
+    class ProjectionCore:
+        @staticmethod
+        def projected_m2(residual, delta, xi, n):
+            residual = np.asarray(residual)
+            delta = np.asarray(delta)
+            xi = np.asarray(xi)
+            calls.append((residual.shape, delta.shape, xi.shape, float(n)))
+            # Keep target/null finite and distinct while exercising downstream indices.
+            return 8.0 if len(calls) == 1 else 12.0
+
+        @staticmethod
+        def chi2_sf(_x, _df):
+            return 0.5
+
+    monkeypatch.setattr(fitstats, "_core_module", lambda: ProjectionCore())
+    monkeypatch.setattr(fitstats, "_projected_m2_numpy", _reject_reference_projection)
+
+    result = _run_multigroup(_multigroup_case())
+
+    assert np.isfinite(result.m2)
+    assert result.m2 == 8.0
+    assert result.null_m2 == 12.0
+    assert len(calls) == 2
+    for residual_shape, delta_shape, xi_shape, n in calls:
+        assert residual_shape == (42,)
+        assert delta_shape[0] == residual_shape[0]
+        assert xi_shape == (42, 42)
+        assert n == 1.0
+
+
+@pytest.mark.parametrize("target", ["residual", "delta", "xi", "n"])
+def test_projected_m2_native_rejects_nonfinite_inputs(target):
+    """The PyO3 boundary must reject non-finite values before Rust arithmetic."""
+    core = fitstats._core_module()
+    if core is None:
+        pytest.skip("compiled Rust core is unavailable in this test environment")
+
+    residual = np.array([0.25, -0.25], dtype=np.float64)
+    delta = np.array([[1.0], [0.5]], dtype=np.float64)
+    xi = np.eye(2, dtype=np.float64)
+    n = 100.0
+    if target == "residual":
+        residual[0] = np.nan
+    elif target == "delta":
+        delta[0, 0] = np.inf
+    elif target == "xi":
+        xi[0, 0] = np.nan
+    else:
+        n = np.inf
+
+    with pytest.raises(ValueError, match="finite"):
+        core.projected_m2(residual, delta, xi, n)
+
+
+def test_projected_m2_native_rejects_oversized_broadcast_before_copy():
+    """Logical broadcast shapes must hit the resource guard before materialization."""
+    core = fitstats._core_module()
+    if core is None:
+        pytest.skip("compiled Rust core is unavailable in this test environment")
+
+    size = 2049
+    residual = np.zeros(size, dtype=np.float64)
+    delta = np.broadcast_to(np.array([[1.0]], dtype=np.float64), (size, 1))
+    xi = np.broadcast_to(np.array([[1.0]], dtype=np.float64), (size, size))
+
+    assert not xi.flags.owndata
+    with pytest.raises(ValueError, match="workspace.*supported element budget"):
+        core.projected_m2(residual, delta, xi, 1.0)