Use Algebraic solver in phase 4

Slightly less efficient but whatever, it needs to complete in a day not week

pipeline.py

@@ -26,6 +26,8 @@ import numpy as np
 from .config import HPCConfig
 from .factored import (
     CAND_TO_QUHIT,
+    MULTS_D,
+    MULTS_M,
     N_CAND_M,
     TOTAL_CANDIDATES,
     FactorResult,
@@ -42,6 +44,7 @@ from .q2k import (
     codes_from_params,
     codes_from_params_batch,
     dequantize_params,
+    dequantize_params_batch,
     derive_subscales,
     make_block,
     pack_codes,
@@ -51,7 +54,7 @@ from .q2k import (
 )
 from .seeds import build_seeds, seed_clamp_value
 from .sieve import block_chain_edges, sieve_measure_chain, subblock_chain_edges
-from .spectral import dc_nudge_codes, objective, refine_codes_hadamard
+from .spectral import dc_nudge_codes, fold_weights, objective, refine_codes_hadamard
 from .viterbi import robust_temperature, viterbi_select
 
 
@@ -372,12 +375,364 @@ def _assemble_block(x, w, seeds, base_d, base_m, cidx, dc_carry, cfg):
     return d, m, ls, lm, codes, best[0], sse
 
 
-def _candidate_config(scales, mins, base_d, base_m, cidx):
-    d16, m16 = candidate_pair(base_d, base_m, cidx)
-    d = float(from_fp16_bits(d16))
-    m = float(from_fp16_bits(m16))
-    ls, lm = derive_subscales(scales, mins, d, m)
-    return d, m, ls, lm
+def _batch_objective(
+    X: np.ndarray,
+    W_f64: np.ndarray,
+    d: np.ndarray,
+    m: np.ndarray,
+    ls: np.ndarray,
+    lm: np.ndarray,
+    codes: np.ndarray,
+    carry: np.ndarray,
+    cfg: HPCConfig,
+) -> tuple[np.ndarray, np.ndarray]:
+    xh = dequantize_params_batch(d, m, ls, lm, codes).astype(np.float64)
+    e = X.astype(np.float64) - xh
+    sse = np.sum(W_f64 * e * e, axis=1)
+    if not cfg.use_spectral or (cfg.dc_lambda == 0.0 and cfg.vw_lambda == 0.0):
+        return sse, sse
+    fw = fold_weights(QK_K, cfg)
+    pen = np.zeros(X.shape[0], dtype=np.float64)
+    src = e
+    m_len = QK_K // 2
+    lvl = 0
+    while m_len >= 1:
+        dst = src[:, :m_len] + src[:, m_len : 2 * m_len]
+        if m_len == 1:
+            dst = dst + carry[:, None].astype(np.float64)
+        if fw[lvl] != 0.0:
+            pen += fw[lvl] * np.sum(dst * dst, axis=1)
+        src = dst
+        m_len //= 2
+        lvl += 1
+    return sse + pen, sse
+
+
+def _batch_extended_dm_solve_and_ulp(
+    X_f64: np.ndarray,
+    W_f64: np.ndarray,
+    ls: np.ndarray,
+    lm: np.ndarray,
+    codes: np.ndarray,
+    d_cur: np.ndarray,
+    m_cur: np.ndarray,
+    carry: np.ndarray,
+    cfg: HPCConfig,
+) -> tuple[np.ndarray, np.ndarray]:
+    """Closed-form 2x2 Cramer solve + O(1) quadratic ULP polish on the full extended objective."""
+    N = X_f64.shape[0]
+    A = np.repeat(ls.astype(np.float64), 16, axis=1) * codes.astype(np.float64)
+    B = np.repeat(lm.astype(np.float64), 16, axis=1)
+    Saa = np.sum(W_f64 * A * A, axis=1)
+    Sab = np.sum(W_f64 * A * B, axis=1)
+    Sbb = np.sum(W_f64 * B * B, axis=1)
+    Sxa = np.sum(W_f64 * X_f64 * A, axis=1)
+    Sxb = np.sum(W_f64 * X_f64 * B, axis=1)
+
+    if cfg.use_spectral:
+        if cfg.vw_lambda > 0.0:
+            vw_w = float(cfg.vw_lambda) / QK_K
+            half = QK_K // 2
+            Ap = A[:, :half] + A[:, half:]
+            Bp = B[:, :half] + B[:, half:]
+            Xp = X_f64[:, :half] + X_f64[:, half:]
+            Saa += vw_w * np.sum(Ap * Ap, axis=1)
+            Sab += vw_w * np.sum(Ap * Bp, axis=1)
+            Sbb += vw_w * np.sum(Bp * Bp, axis=1)
+            Sxa += vw_w * np.sum(Xp * Ap, axis=1)
+            Sxb += vw_w * np.sum(Xp * Bp, axis=1)
+        if cfg.dc_lambda > 0.0:
+            dc_w = float(cfg.dc_lambda) / QK_K
+            rA = np.sum(A, axis=1)
+            rB = np.sum(B, axis=1)
+            rS = np.sum(X_f64, axis=1) + carry.astype(np.float64)
+            Saa += dc_w * rA * rA
+            Sab += dc_w * rA * rB
+            Sbb += dc_w * rB * rB
+            Sxa += dc_w * rS * rA
+            Sxb += dc_w * rS * rB
+
+    det = Saa * Sbb - Sab * Sab
+    valid = det > 1.0e-30
+    s_det = np.where(valid, det, 1.0)
+    d_new = (Sbb * Sxa - Sab * Sxb) / s_det
+    m_new = (Sab * Sxa - Saa * Sxb) / s_det
+    d_opt = np.where(valid & (d_new > 0.0) & np.isfinite(d_new), d_new, d_cur)
+    m_opt = np.where(valid & (m_new >= 0.0) & np.isfinite(m_new), m_new, m_cur)
+
+    rad = max(int(cfg.ulp_radius), 1)
+    deltas = np.arange(-rad, rad + 1, dtype=np.int32)
+    d16_c = to_fp16_bits(d_opt.astype(np.float32)).astype(np.int32)
+    m16_c = to_fp16_bits(m_opt.astype(np.float32)).astype(np.int32)
+    d_bits = np.clip(d16_c[:, None] + deltas[None, :], 0, 0x7BFF).astype(np.uint16)
+    m_bits = np.clip(m16_c[:, None] + deltas[None, :], 0, 0x7BFF).astype(np.uint16)
+    dd = from_fp16_bits(d_bits).astype(np.float64)[:, :, None]
+    mm = from_fp16_bits(m_bits).astype(np.float64)[:, None, :]
+
+    q_form = (
+        Saa[:, None, None] * (dd * dd)
+        - 2.0 * Sab[:, None, None] * (dd * mm)
+        + Sbb[:, None, None] * (mm * mm)
+        - 2.0 * Sxa[:, None, None] * dd
+        + 2.0 * Sxb[:, None, None] * mm
+    ).reshape(N, -1)
+    best_idx = np.argmin(q_form, axis=1)
+    K = deltas.size
+    rows = np.arange(N)
+    return dd[rows, best_idx // K, 0].astype(np.float32), mm[rows, 0, best_idx % K].astype(np.float32)
+
+
+def _batch_algebraic_refine_codes(
+    X_f64: np.ndarray,
+    W_f64: np.ndarray,
+    d: np.ndarray,
+    m: np.ndarray,
+    ls: np.ndarray,
+    lm: np.ndarray,
+    codes: np.ndarray,
+    carry: np.ndarray,
+    cfg: HPCConfig,
+) -> np.ndarray:
+    """Closed-form 1D quadratic coordinate projection for Vesica + DC + weighted SSE."""
+    N = X_f64.shape[0]
+    out = codes.astype(np.int64).copy()
+    lam_vw = (float(cfg.vw_lambda) / QK_K) if cfg.use_spectral else 0.0
+    lam_dc = (float(cfg.dc_lambda) / QK_K) if cfg.use_spectral else 0.0
+    if lam_vw == 0.0 and lam_dc == 0.0:
+        return out.astype(np.uint8)
+
+    d_sub = np.repeat(d.astype(np.float64)[:, None] * ls.astype(np.float64), 16, axis=1)
+    m_sub = np.repeat(m.astype(np.float64)[:, None] * lm.astype(np.float64), 16, axis=1)
+    active = d_sub >= 1.0e-15
+    e = X_f64 - (d_sub * out - m_sub)
+    sse_live = np.sum(W_f64 * e * e, axis=1)
+    sse_cap = sse_live * (1.0 + cfg.sse_budget)
+    half = QK_K // 2
+    v = e[:, :half] + e[:, half:]
+    dc = np.sum(e, axis=1) + carry.astype(np.float64)
+    denom = (W_f64 + lam_vw + lam_dc) * (d_sub * d_sub)
+    safe_den = np.where(active & (denom > 1.0e-30), denom, 1.0)
+
+    v_full = np.concatenate([v, v], axis=1)
+    grad = d_sub * (W_f64 * e + lam_vw * v_full + lam_dc * dc[:, None])
+    dq = np.where(active, nearest_int((grad / safe_den).astype(np.float32)), 0)
+    q_new = np.clip(out + dq, 0, 3)
+
+    for b in range(N):
+        changed = np.where(q_new[b] != out[b])[0]
+        if changed.size == 0:
+            continue
+        xb, wb, dsb, msb = X_f64[b], W_f64[b], d_sub[b], m_sub[b]
+        eb, vb, dcb, sseb, capb = e[b], v[b], float(dc[b]), float(sse_live[b]), float(sse_cap[b])
+        e_cand = xb[changed] - (dsb[changed] * q_new[b, changed] - msb[changed])
+        de_cand = e_cand - eb[changed]
+        dsse_cand = wb[changed] * (e_cand * e_cand - eb[changed] * eb[changed])
+        pi_cand = np.where(changed < half, changed, changed - half)
+        v_old = vb[pi_cand]
+        v_new = v_old + de_cand
+        d_spec = lam_vw * (v_old * v_old - v_new * v_new) + lam_dc * (dcb * dcb - (dcb + de_cand) ** 2)
+        order = changed[np.argsort(-(d_spec - dsse_cand))]
+        for k in order:
+            q_try = int(q_new[b, k])
+            if q_try == int(out[b, k]):
+                continue
+            e_new = xb[k] - (dsb[k] * q_try - msb[k])
+            de = e_new - eb[k]
+            dsse = wb[k] * (e_new * e_new - eb[k] * eb[k])
+            if sseb + dsse > capb:
+                continue
+            p = k if k < half else k - half
+            v_np = vb[p] + de
+            dc_n = dcb + de
+            spec_gain = lam_vw * (vb[p] * vb[p] - v_np * v_np) + lam_dc * (dcb * dcb - dc_n * dc_n)
+            if spec_gain > dsse + 1.0e-18:
+                out[b, k] = q_try
+                eb[k] = e_new
+                vb[p] = v_np
+                dcb = dc_n
+                sseb += dsse
+    return out.astype(np.uint8)
+
+
+def _refine_trajectory_algebraic(
+    X: np.ndarray,
+    X_f64: np.ndarray,
+    W_f64: np.ndarray,
+    Xb: np.ndarray,
+    Wb: np.ndarray,
+    Sw: np.ndarray,
+    Sx: np.ndarray,
+    safe_Sw: np.ndarray,
+    d0: np.ndarray,
+    m0: np.ndarray,
+    ls0: np.ndarray,
+    lm0: np.ndarray,
+    codes0: np.ndarray,
+    carry: np.ndarray,
+    cfg: HPCConfig,
+):
+    N = X.shape[0]
+    best_d, best_m = d0.copy(), m0.copy()
+    best_ls, best_lm = ls0.copy(), lm0.copy()
+    best_codes = codes0.copy()
+    best_obj, best_sse = _batch_objective(X, W_f64, best_d, best_m, best_ls, best_lm, best_codes, carry, cfg)
+    d, m, ls, lm, codes = best_d.copy(), best_m.copy(), best_ls.copy(), best_lm.copy(), best_codes.copy()
+    ls_deltas = np.array([-3, -2, -1, 0, 1, 2, 3], dtype=np.int64)
+    lm_deltas = np.array([-1, 0, 1], dtype=np.int64)
+    n_idx = np.arange(N)[:, None]
+    s_idx = np.arange(NS)[None, :]
+
+    for _ in range(cfg.refine_passes):
+        d_new, m_new = _batch_extended_dm_solve_and_ulp(X_f64, W_f64, ls, lm, codes, d, m, carry, cfg)
+        obj_dm, sse_dm = _batch_objective(X, W_f64, d_new, m_new, ls, lm, codes, carry, cfg)
+        imp_dm = obj_dm < best_obj - 1.0e-15
+        if np.any(imp_dm):
+            best_obj = np.where(imp_dm, obj_dm, best_obj)
+            best_sse = np.where(imp_dm, sse_dm, best_sse)
+            best_d = np.where(imp_dm, d_new, best_d)
+            best_m = np.where(imp_dm, m_new, best_m)
+            best_ls = np.where(imp_dm[:, None], ls, best_ls)
+            best_lm = np.where(imp_dm[:, None], lm, best_lm)
+            best_codes = np.where(imp_dm[:, None], codes, best_codes)
+        d, m = best_d.copy(), best_m.copy()
+
+        codes_re = codes_from_params_batch(X, d, m, ls, lm)
+        obj_re, sse_re = _batch_objective(X, W_f64, d, m, ls, lm, codes_re, carry, cfg)
+        imp_re = obj_re < best_obj - 1.0e-15
+        if np.any(imp_re):
+            best_obj = np.where(imp_re, obj_re, best_obj)
+            best_sse = np.where(imp_re, sse_re, best_sse)
+            best_codes = np.where(imp_re[:, None], codes_re, best_codes)
+        codes = best_codes.copy()
+
+        q = codes.astype(np.float64).reshape(N, NS, 16)
+        Sqq = np.sum(Wb * q * q, axis=2)
+        Swq = np.sum(Wb * q, axis=2)
+        Sqx = np.sum(Wb * q * Xb, axis=2)
+        q_bar = Swq / safe_Sw
+        Vqq = np.maximum(Sqq - Swq * q_bar, 0.0)
+        Cqx = Sqx - Sx * q_bar
+        s_star = np.where(Vqq > 1.0e-30, np.maximum(Cqx / np.maximum(Vqq, 1.0e-30), 0.0), d[:, None] * ls)
+        sd = np.where(d[:, None] > 1.0e-15, d[:, None], 1.0)
+        sm = np.where(m[:, None] > 1.0e-15, m[:, None], 1.0)
+        ls_center = np.where(d[:, None] > 1.0e-15, np.clip(nearest_int((s_star / sd).astype(np.float32)), 0, 15), 0)
+
+        c_ls_base = np.clip(ls_center[:, :, None] + ls_deltas[None, None, :], 0, 15).astype(np.float64)
+        d_sub_base = d.astype(np.float64)[:, None, None] * c_ls_base
+        t_opt_base = np.maximum(0.0, (d_sub_base * Swq[:, :, None] - Sx[:, :, None]) / safe_Sw[:, :, None])
+        lm_center = np.where(
+            m[:, None, None] > 1.0e-15,
+            np.clip(nearest_int((t_opt_base / sm[:, :, None]).astype(np.float32)), 0, 15),
+            0,
+        )
+
+        cand_ls = np.repeat(c_ls_base, 3, axis=2)
+        cand_lm = np.clip(
+            np.repeat(lm_center, 3, axis=2) + np.tile(lm_deltas, 7)[None, None, :], 0, 15
+        ).astype(np.float64)
+
+        d_sub_c = d.astype(np.float64)[:, None, None] * cand_ls
+        m_sub_c = m.astype(np.float64)[:, None, None] * cand_lm
+        safe_dc = d_sub_c >= 1.0e-15
+        sd_c = np.where(safe_dc, d_sub_c, 1.0)
+        q_c = np.where(
+            safe_dc[:, :, :, None],
+            np.clip(
+                nearest_int(((Xb[:, :, None, :] + m_sub_c[:, :, :, None]) / sd_c[:, :, :, None]).astype(np.float32)),
+                0,
+                3,
+            ),
+            0,
+        ).astype(np.float64)
+        e_c = Xb[:, :, None, :] - (d_sub_c[:, :, :, None] * q_c - m_sub_c[:, :, :, None])
+        err_c = np.sum(Wb[:, :, None, :] * e_c * e_c, axis=3)
+        pick = np.argmin(err_c, axis=2)
+        ls_try = cand_ls[n_idx, s_idx, pick].astype(np.uint8)
+        lm_try = cand_lm[n_idx, s_idx, pick].astype(np.uint8)
+        codes_try = codes_from_params_batch(X, d, m, ls_try, lm_try)
+        d_try, m_try = _batch_extended_dm_solve_and_ulp(X_f64, W_f64, ls_try, lm_try, codes_try, d, m, carry, cfg)
+        obj_try, sse_try = _batch_objective(X, W_f64, d_try, m_try, ls_try, lm_try, codes_try, carry, cfg)
+        imp_sub = obj_try < best_obj - 1.0e-15
+        if np.any(imp_sub):
+            best_obj = np.where(imp_sub, obj_try, best_obj)
+            best_sse = np.where(imp_sub, sse_try, best_sse)
+            best_d = np.where(imp_sub, d_try, best_d)
+            best_m = np.where(imp_sub, m_try, best_m)
+            best_ls = np.where(imp_sub[:, None], ls_try, best_ls)
+            best_lm = np.where(imp_sub[:, None], lm_try, best_lm)
+            best_codes = np.where(imp_sub[:, None], codes_try, best_codes)
+        d, m, ls, lm, codes = best_d.copy(), best_m.copy(), best_ls.copy(), best_lm.copy(), best_codes.copy()
+
+    refined = _batch_algebraic_refine_codes(X_f64, W_f64, best_d, best_m, best_ls, best_lm, best_codes, carry, cfg)
+    d_r, m_r = _batch_extended_dm_solve_and_ulp(X_f64, W_f64, best_ls, best_lm, refined, best_d, best_m, carry, cfg)
+    obj_r, sse_r = _batch_objective(X, W_f64, d_r, m_r, best_ls, best_lm, refined, carry, cfg)
+    imp_r = obj_r < best_obj - 1.0e-15
+    if np.any(imp_r):
+        best_obj = np.where(imp_r, obj_r, best_obj)
+        best_sse = np.where(imp_r, sse_r, best_sse)
+        best_d = np.where(imp_r, d_r, best_d)
+        best_m = np.where(imp_r, m_r, best_m)
+        best_codes = np.where(imp_r[:, None], refined, best_codes)
+
+    return best_d, best_m, best_ls, best_lm, best_codes, best_obj, best_sse
+
+
+def _phase4_algebraic_batch(
+    X: np.ndarray,
+    W: np.ndarray | None,
+    scales_seed: np.ndarray,
+    mins_seed: np.ndarray,
+    base_d: np.ndarray,
+    base_m: np.ndarray,
+    best: np.ndarray,
+    carry: np.ndarray,
+    standalone: list[FactorResult],
+    cfg: HPCConfig,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, int, int]:
+    """Batched algebraic Phase-4 solver across N blocks (replaces brute-force 16x16 & ULP grids)."""
+    N = X.shape[0]
+    X_f64 = X.astype(np.float64)
+    W_f64 = np.ones_like(X_f64) if W is None else W.astype(np.float64)
+    Xb = X_f64.reshape(N, NS, 16)
+    Wb = W_f64.reshape(N, NS, 16)
+    Sw = np.sum(Wb, axis=2)
+    Sx = np.sum(Wb * Xb, axis=2)
+    safe_Sw = np.maximum(Sw, 1.0e-30)
+
+    di = best // N_CAND_M
+    mi = best % N_CAND_M
+    d_v = round_fp16(base_d * MULTS_D[di])
+    m_v = round_fp16(base_m * MULTS_M[mi])
+    sd_v = np.where(d_v[:, None] > 1.0e-15, d_v[:, None], 1.0)
+    sm_v = np.where(m_v[:, None] > 1.0e-15, m_v[:, None], 1.0)
+    ls_v = np.where(d_v[:, None] > 1.0e-15, np.clip(nearest_int(scales_seed / sd_v), 0, 15), 0).astype(np.uint8)
+    lm_v = np.where(m_v[:, None] > 1.0e-15, np.clip(nearest_int(mins_seed / sm_v), 0, 15), 0).astype(np.uint8)
+    codes_v = codes_from_params_batch(X, d_v, m_v, ls_v, lm_v)
+
+    d1, m1, ls1, lm1, c1, obj1, _ = _refine_trajectory_algebraic(
+        X, X_f64, W_f64, Xb, Wb, Sw, Sx, safe_Sw, d_v, m_v, ls_v, lm_v, codes_v, carry, cfg
+    )
+    cand_hits = 0
+    fac_hits = 0
+
+    if cfg.use_factor_floor:
+        d_f = np.array([f.d for f in standalone], dtype=np.float32)
+        m_f = np.array([f.dmin for f in standalone], dtype=np.float32)
+        ls_f = np.stack([f.ls for f in standalone], axis=0)
+        lm_f = np.stack([f.lm for f in standalone], axis=0)
+        codes_f = np.stack([f.codes for f in standalone], axis=0)
+        d2, m2, ls2, lm2, c2, obj2, _ = _refine_trajectory_algebraic(
+            X, X_f64, W_f64, Xb, Wb, Sw, Sx, safe_Sw, d_f, m_f, ls_f, lm_f, codes_f, carry, cfg
+        )
+        take2 = obj2 < obj1 - 1.0e-15
+        fac_hits = int(np.sum(take2))
+        d1 = np.where(take2, d2, d1)
+        m1 = np.where(take2, m2, m1)
+        ls1 = np.where(take2[:, None], ls2, ls1)
+        lm1 = np.where(take2[:, None], lm2, lm1)
+        c1 = np.where(take2[:, None], c2, c1)
+
+    return d1, m1, ls1, lm1, c1, cand_hits, fac_hits
 
 
 def _apply_dc_nudge(X, W, blocks, cfg, row_width):
@@ -617,62 +972,6 @@ class Q2KQuantizer:
                             eta=eta,
                             sigma=sigma,
                         )
-
-                _check_stop()
-                _emit(
-                    "Phase 3: 6-State Sieve & 36-State Viterbi DP",
-                    0,
-                    n_blocks,
-                    f"Running Sieve + 36-state Viterbi sequence optimization ({cfg.optimizer.upper()})",
-                    sigma=sigma,
-                )
-                best, temperature, _ = _phase3_select(candidate_errors, cfg)
-                carry = _compute_dc_carry(X, seeds, best, base_d, base_m, cfg, row_width)
-                _emit(
-                    "Phase 3.9: Viterbi Path & Rolling DC Carry",
-                    n_blocks,
-                    n_blocks,
-                    f"Boltzmann temp={temperature:.4e}, max |DC carry|={float(np.max(np.abs(carry))):.4e}",
-                    temperature=temperature,
-                    sigma=sigma,
-                )
-
-                t_asm = time.time()
-                last_emit = 0.0
-                p4_args = [
-                    (i, int(best[i]), float(carry[i]), standalone[i])
-                    for i in range(n_blocks)
-                ]
-                for i, (blk_obj, cand_hit, fac_hit, blk_err) in enumerate(
-                    pool.imap(_phase4_one, p4_args, chunksize=chunksize)
-                ):
-                    _check_stop()
-                    blocks.append(blk_obj)
-                    candidate_floor_used += cand_hit
-                    factor_floor_used += fac_hit
-                    running_sse += blk_err
-                    now = time.time()
-                    if (i + 1) % step_report == 0 or (i + 1) == n_blocks or (now - last_emit) >= 0.25:
-                        last_emit = now
-                        elapsed = max(now - t_asm, 1.0e-6)
-                        rate = (i + 1) / elapsed
-                        eta = (n_blocks - (i + 1)) / max(rate, 1.0e-9)
-                        cur_rmse = float(np.sqrt(running_sse / ((i + 1) * QK_K)))
-                        _emit(
-                            f"Phase 4: Sub-block Sieve + ULP + Hadamard Assembly ({n_workers} workers)",
-                            i + 1,
-                            n_blocks,
-                            f"blk {i+1}/{n_blocks} ({rate:.1f} blk/s, ETA {eta:.1f}s) "
-                            f"rmse={cur_rmse:.6e} rmse/sigma={cur_rmse/max(sigma, 1e-15):.4f} "
-                            f"floors(cand={candidate_floor_used}, fac={factor_floor_used})",
-                            rmse=cur_rmse,
-                            sigma=sigma,
-                            temperature=temperature,
-                            cand_floors=candidate_floor_used,
-                            fac_floors=factor_floor_used,
-                            rate=rate,
-                            eta=eta,
-                        )
                 pool.close()
                 pool.join()
             except BaseException:
@@ -703,61 +1002,72 @@ class Q2KQuantizer:
                         sigma=sigma,
                     )
 
+        _MP_STATE.clear()
+        _check_stop()
+        _emit(
+            "Phase 3: 6-State Sieve & 36-State Viterbi DP",
+            0,
+            n_blocks,
+            f"Running Sieve + 36-state Viterbi sequence optimization ({cfg.optimizer.upper()})",
+            sigma=sigma,
+        )
+        best, temperature, _ = _phase3_select(candidate_errors, cfg)
+        carry = _compute_dc_carry(X, seeds, best, base_d, base_m, cfg, row_width)
+        _emit(
+            "Phase 3.9: Viterbi Path & Rolling DC Carry",
+            n_blocks,
+            n_blocks,
+            f"Boltzmann temp={temperature:.4e}, max |DC carry|={float(np.max(np.abs(carry))):.4e}",
+            temperature=temperature,
+            sigma=sigma,
+        )
+
+        t_asm = time.time()
+        batch_sz = 256
+        for start in range(0, n_blocks, batch_sz):
             _check_stop()
-            _emit(
-                "Phase 3: 6-State Sieve & 36-State Viterbi DP",
-                0,
-                n_blocks,
-                f"Running Sieve + 36-state Viterbi sequence optimization ({cfg.optimizer.upper()})",
-                sigma=sigma,
+            stop = min(start + batch_sz, n_blocks)
+            W_chunk = None if W is None else W[start:stop]
+            d_b, m_b, ls_b, lm_b, codes_b, c_hits, f_hits = _phase4_algebraic_batch(
+                X[start:stop],
+                W_chunk,
+                seeds["scales"][start:stop],
+                seeds["mins"][start:stop],
+                base_d[start:stop],
+                base_m[start:stop],
+                best[start:stop],
+                carry[start:stop],
+                standalone[start:stop],
+                cfg,
             )
-            best, temperature, _ = _phase3_select(candidate_errors, cfg)
-            carry = _compute_dc_carry(X, seeds, best, base_d, base_m, cfg, row_width)
+            candidate_floor_used += c_hits
+            factor_floor_used += f_hits
+            for k in range(stop - start):
+                blk_obj = make_block(d_b[k], m_b[k], ls_b[k], lm_b[k], codes_b[k])
+                blocks.append(blk_obj)
+                running_sse += blk_obj.error(X[start + k])
+
+            now = time.time()
+            elapsed = max(now - t_asm, 1.0e-6)
+            rate = stop / elapsed
+            eta = (n_blocks - stop) / max(rate, 1.0e-9)
+            cur_rmse = float(np.sqrt(running_sse / (stop * QK_K)))
             _emit(
-                "Phase 3.9: Viterbi Path & Rolling DC Carry",
-                n_blocks,
+                "Phase 4: Algebraic Coupled Bilinear Assembly",
+                stop,
                 n_blocks,
-                f"Boltzmann temp={temperature:.4e}, max |DC carry|={float(np.max(np.abs(carry))):.4e}",
-                temperature=temperature,
+                f"blk {stop}/{n_blocks} ({rate:.1f} blk/s, ETA {eta:.1f}s) "
+                f"rmse={cur_rmse:.6e} rmse/sigma={cur_rmse/max(sigma, 1e-15):.4f} "
+                f"floors(cand={candidate_floor_used}, fac={factor_floor_used})",
+                rmse=cur_rmse,
                 sigma=sigma,
+                temperature=temperature,
+                cand_floors=candidate_floor_used,
+                fac_floors=factor_floor_used,
+                rate=rate,
+                eta=eta,
             )
 
-            t_asm = time.time()
-            last_emit = 0.0
-            for i in range(n_blocks):
-                _check_stop()
-                blk_obj, cand_hit, fac_hit, blk_err = _phase4_one(
-                    (i, int(best[i]), float(carry[i]), standalone[i])
-                )
-                blocks.append(blk_obj)
-                candidate_floor_used += cand_hit
-                factor_floor_used += fac_hit
-                running_sse += blk_err
-                now = time.time()
-                if (i + 1) % step_report == 0 or (i + 1) == n_blocks or (now - last_emit) >= 0.25:
-                    last_emit = now
-                    elapsed = max(now - t_asm, 1.0e-6)
-                    rate = (i + 1) / elapsed
-                    eta = (n_blocks - (i + 1)) / max(rate, 1.0e-9)
-                    cur_rmse = float(np.sqrt(running_sse / ((i + 1) * QK_K)))
-                    _emit(
-                        "Phase 4: Sub-block Sieve + ULP + Hadamard Assembly",
-                        i + 1,
-                        n_blocks,
-                        f"blk {i+1}/{n_blocks} ({rate:.1f} blk/s, ETA {eta:.1f}s) "
-                        f"rmse={cur_rmse:.6e} rmse/sigma={cur_rmse/max(sigma, 1e-15):.4f} "
-                        f"floors(cand={candidate_floor_used}, fac={factor_floor_used})",
-                        rmse=cur_rmse,
-                        sigma=sigma,
-                        temperature=temperature,
-                        cand_floors=candidate_floor_used,
-                        fac_floors=factor_floor_used,
-                        rate=rate,
-                        eta=eta,
-                    )
-
-        _MP_STATE.clear()
-
         if cfg.apply_dc_carry and (cfg.dc_lambda > 0.0 or cfg.dc_decay > 0.0):
             _check_stop()
             _emit(