Upload 16 files

config.py

@@ -74,6 +74,7 @@ class HPCConfig:
     clamp_kurtosis_soft: float = 20.0
     clamp_sigma_soft: float = 6.0
     clamp_sigma_hard: float = 4.0
+    temper_imatrix: bool = True
 
     # ---- reporting & parallelism ----------------------------------------
     device: str = "auto"

cuda_backend.py

@@ -223,10 +223,31 @@ def _make_qp_cuda(values: "torch.Tensor", weights: "torch.Tensor", nmax: int = 1
     return d, codes.to(torch.uint8)
 
 
-def _build_seeds_cuda(X: "torch.Tensor", W: "torch.Tensor | None", clamp_val: float):
+def _block_clamp_cuda(X: "torch.Tensor", cfg: HPCConfig) -> "torch.Tensor":
+    """Per-block kurtosis-gated outlier clamp thresholds on GPU (shape (B, 1, 1))."""
+    b2 = X * X
+    mean_sq = torch.mean(b2, dim=1)
+    mean_4 = torch.mean(b2 * b2, dim=1)
+    kurt = torch.where(mean_sq > 1.0e-30, mean_4 / torch.clamp(mean_sq * mean_sq, min=1.0e-30), torch.full_like(mean_sq, 3.0))
+    sigma = torch.sqrt(torch.clamp(mean_sq, min=0.0))
+    mult = torch.where(
+        kurt <= cfg.clamp_kurtosis_off,
+        torch.full_like(kurt, 1.0e9),
+        torch.where(
+            kurt <= cfg.clamp_kurtosis_soft,
+            torch.full_like(kurt, float(cfg.clamp_sigma_soft)),
+            torch.full_like(kurt, float(cfg.clamp_sigma_hard)),
+        ),
+    )
+    clamp = torch.where(mean_sq > 1.0e-30, sigma * mult, torch.full_like(sigma, float("inf")))
+    return clamp[:, None, None]
+
+
+def _build_seeds_cuda(X: "torch.Tensor", W: "torch.Tensor | None", cfg: HPCConfig):
     B = X.shape[0]
     blocks = X.reshape(B, 16, 16)
-    clipped = torch.clamp(blocks, min=-clamp_val, max=clamp_val)
+    per_blk_clamp = _block_clamp_cuda(X, cfg)
+    clipped = torch.clamp(blocks, min=-per_blk_clamp, max=per_blk_clamp)
     if W is None:
         sigma_sq = torch.mean(blocks * blocks, dim=(1, 2), keepdim=True)
         weights = torch.sqrt(sigma_sq + clipped * clipped)
@@ -248,19 +269,19 @@ def _build_seeds_cuda(X: "torch.Tensor", W: "torch.Tensor | None", clamp_val: fl
         "dmin": dmin,
         "ls": ls,
         "lm": lm,
+        "clamp": per_blk_clamp.reshape(B, 1),
     }
 
 
-def _phase2_wls_cuda(X: "torch.Tensor", W: "torch.Tensor", seeds: dict, clamp_val: float):
-    clipped = torch.clamp(X, min=-clamp_val, max=clamp_val)
+def _phase2_wls_cuda(X: "torch.Tensor", W: "torch.Tensor", seeds: dict):
+    per_blk_clamp = seeds["clamp"]
+    clipped = torch.clamp(X, min=-per_blk_clamp, max=per_blk_clamp)
     base_d = seeds["d"].clone()
     base_m = seeds["dmin"].clone()
     ls = seeds["ls"].clone()
     lm = seeds["lm"].clone()
     scales = seeds["scales"]
     mins = seeds["mins"]
-    seed_d = seeds["d"]
-    seed_m = seeds["dmin"]
 
     for _ in range(5):
         codes = _codes_from_params_t(clipped, base_d, base_m, ls, lm).to(torch.float32)
@@ -276,8 +297,8 @@ def _phase2_wls_cuda(X: "torch.Tensor", W: "torch.Tensor", seeds: dict, clamp_va
         s_det = torch.where(valid, det, torch.ones_like(det))
         d_new = (Sbb * Sxa - Sab * Sxb) / s_det
         m_new = (Sab * Sxa - Saa * Sxb) / s_det
-        upd_d = valid & (d_new > 0.0) & (d_new < 4.0 * (seed_d + 1.0e-10))
-        upd_m = valid & (m_new > 0.0) & (m_new < 4.0 * (seed_m + 1.0e-10))
+        upd_d = valid & (d_new > 0.0) & torch.isfinite(d_new)
+        upd_m = valid & (m_new > 0.0) & torch.isfinite(m_new)
         base_d = torch.where(upd_d, _round_fp16_t(d_new), base_d)
         base_m = torch.where(upd_m, _round_fp16_t(m_new), base_m)
         sd = torch.where(base_d > 1.0e-15, base_d, torch.ones_like(base_d))
@@ -584,10 +605,10 @@ def _standalone_coupled_solve_cuda(
             d = _round_fp16_t(d)
             m = _round_fp16_t(m)
             codes = _codes_from_params_t(X, d, m, ls, lm)
-            obj, _ = _objective_t(X, W, d, m, ls, lm, codes, None, cfg)
-            take = obj < best_obj - 1.0e-12
+            _, sse_cur = _objective_t(X, W, d, m, ls, lm, codes, None, cfg)
+            take = sse_cur < best_obj - 1.0e-12
             if torch.any(take):
-                best_obj = torch.where(take, obj, best_obj)
+                best_obj = torch.where(take, sse_cur, best_obj)
                 best_d = torch.where(take, d, best_d)
                 best_m = torch.where(take, m, best_m)
                 best_ls = torch.where(take[:, None], ls, best_ls)
@@ -617,29 +638,6 @@ def _extended_dm_solve_and_ulp_cuda(
     Sxa = torch.sum(W * X * A, dim=1)
     Sxb = torch.sum(W * X * Bm, dim=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 = Bm[:, :half] + Bm[:, half:]
-            Xp = X[:, :half] + X[:, half:]
-            Saa = Saa + vw_w * torch.sum(Ap * Ap, dim=1)
-            Sab = Sab + vw_w * torch.sum(Ap * Bp, dim=1)
-            Sbb = Sbb + vw_w * torch.sum(Bp * Bp, dim=1)
-            Sxa = Sxa + vw_w * torch.sum(Xp * Ap, dim=1)
-            Sxb = Sxb + vw_w * torch.sum(Xp * Bp, dim=1)
-        if cfg.dc_lambda > 0.0:
-            dc_w = float(cfg.dc_lambda) / QK_K
-            rA = torch.sum(A, dim=1)
-            rB = torch.sum(Bm, dim=1)
-            rS = torch.sum(X, dim=1) + carry
-            Saa = Saa + dc_w * rA * rA
-            Sab = Sab + dc_w * rA * rB
-            Sbb = Sbb + dc_w * rB * rB
-            Sxa = Sxa + dc_w * rS * rA
-            Sxb = Sxb + dc_w * rS * rB
-
     det = Saa * Sbb - Sab * Sab
     valid = det > 1.0e-30
     s_det = torch.where(valid, det, torch.ones_like(det))
@@ -691,7 +689,7 @@ def _refine_trajectory_cuda(
     best_d, best_m = d0.clone(), m0.clone()
     best_ls, best_lm = ls0.clone(), lm0.clone()
     best_codes = codes0.clone()
-    best_obj, _ = _objective_t(X, W, best_d, best_m, best_ls, best_lm, best_codes, carry, cfg)
+    _, best_sse = _objective_t(X, W, best_d, best_m, best_ls, best_lm, best_codes, carry, cfg)
     d, m, ls, lm, codes = best_d.clone(), best_m.clone(), best_ls.clone(), best_lm.clone(), best_codes.clone()
 
     ls_deltas = torch.tensor([-3, -2, -1, 0, 1, 2, 3], dtype=torch.int64, device=device)
@@ -701,10 +699,10 @@ def _refine_trajectory_cuda(
 
     for _ in range(cfg.refine_passes):
         d_new, m_new = _extended_dm_solve_and_ulp_cuda(X, W, ls, lm, codes, d, m, carry, cfg)
-        obj_dm, _ = _objective_t(X, W, d_new, m_new, ls, lm, codes, carry, cfg)
-        imp_dm = obj_dm < best_obj - 1.0e-12
+        _, sse_dm = _objective_t(X, W, d_new, m_new, ls, lm, codes, carry, cfg)
+        imp_dm = sse_dm < best_sse - 1.0e-12
         if torch.any(imp_dm):
-            best_obj = torch.where(imp_dm, obj_dm, best_obj)
+            best_sse = torch.where(imp_dm, sse_dm, best_sse)
             best_d = torch.where(imp_dm, d_new, best_d)
             best_m = torch.where(imp_dm, m_new, best_m)
             best_ls = torch.where(imp_dm[:, None], ls, best_ls)
@@ -713,10 +711,10 @@ def _refine_trajectory_cuda(
         d, m = best_d.clone(), best_m.clone()
 
         codes_re = _codes_from_params_t(X, d, m, ls, lm)
-        obj_re, _ = _objective_t(X, W, d, m, ls, lm, codes_re, carry, cfg)
-        imp_re = obj_re < best_obj - 1.0e-12
+        _, sse_re = _objective_t(X, W, d, m, ls, lm, codes_re, carry, cfg)
+        imp_re = sse_re < best_sse - 1.0e-12
         if torch.any(imp_re):
-            best_obj = torch.where(imp_re, obj_re, best_obj)
+            best_sse = torch.where(imp_re, sse_re, best_sse)
             best_codes = torch.where(imp_re[:, None], codes_re, best_codes)
         codes = best_codes.clone()
 
@@ -766,10 +764,10 @@ def _refine_trajectory_cuda(
         lm_try = cand_lm[b_idx, s_idx, pick].to(torch.uint8)
         codes_try = _codes_from_params_t(X, d, m, ls_try, lm_try)
         d_try, m_try = _extended_dm_solve_and_ulp_cuda(X, W, ls_try, lm_try, codes_try, d, m, carry, cfg)
-        obj_try, _ = _objective_t(X, W, d_try, m_try, ls_try, lm_try, codes_try, carry, cfg)
-        imp_sub = obj_try < best_obj - 1.0e-12
+        _, sse_try = _objective_t(X, W, d_try, m_try, ls_try, lm_try, codes_try, carry, cfg)
+        imp_sub = sse_try < best_sse - 1.0e-12
         if torch.any(imp_sub):
-            best_obj = torch.where(imp_sub, obj_try, best_obj)
+            best_sse = torch.where(imp_sub, sse_try, best_sse)
             best_d = torch.where(imp_sub, d_try, best_d)
             best_m = torch.where(imp_sub, m_try, best_m)
             best_ls = torch.where(imp_sub[:, None], ls_try, best_ls)
@@ -777,7 +775,7 @@ def _refine_trajectory_cuda(
             best_codes = torch.where(imp_sub[:, None], codes_try, best_codes)
         d, m, ls, lm, codes = best_d.clone(), best_m.clone(), best_ls.clone(), best_lm.clone(), best_codes.clone()
 
-    return best_d, best_m, best_ls, best_lm, best_codes, best_obj
+    return best_d, best_m, best_ls, best_lm, best_codes, best_sse
 
 
 def _compute_dc_carry_fast(
@@ -1003,13 +1001,19 @@ def quantize_cuda(
             _check_stop()
             stop = min(start + batch_blocks, n_blocks)
             X_t = torch.from_numpy(X_np[start:stop]).to(dev)
-            W_t = (
-                torch.ones_like(X_t)
-                if W_np is None
-                else torch.from_numpy(W_np[start:stop]).to(dev)
-            )
-            seeds_t = _build_seeds_cuda(X_t, None if W_np is None else W_t, clamp_val)
-            bd_t, bm_t = _phase2_wls_cuda(X_t, W_t, seeds_t, clamp_val)
+            if W_np is None:
+                W_t = torch.ones_like(X_t)
+            else:
+                W_raw_t = torch.from_numpy(W_np[start:stop]).to(dev)
+                if cfg.temper_imatrix:
+                    sig2_t = torch.mean(X_t * X_t, dim=1, keepdim=True)
+                    W_t = torch.sqrt(torch.clamp(W_raw_t, min=1.0e-6)) * torch.sqrt(sig2_t + X_t * X_t)
+                    wm_t = torch.mean(W_t, dim=1, keepdim=True)
+                    W_t = torch.where(wm_t > 1.0e-30, W_t / torch.clamp(wm_t, min=1.0e-30), torch.ones_like(W_t))
+                else:
+                    W_t = W_raw_t
+            seeds_t = _build_seeds_cuda(X_t, None if W_np is None else W_t, cfg)
+            bd_t, bm_t = _phase2_wls_cuda(X_t, W_t, seeds_t)
 
             scales_seed[start:stop] = seeds_t["scales"].cpu().numpy()
             mins_seed[start:stop] = seeds_t["mins"].cpu().numpy()
@@ -1020,12 +1024,13 @@ def quantize_cuda(
                 X_t, W_t, seeds_t["scales"], seeds_t["mins"], bd_t, bm_t, cfg
             )
             sd_t, sm_t, sls_t, slm_t, sc_t, sobj_t = _standalone_coupled_solve_cuda(X_t, W_t, cfg)
+            _, sse_Solo = _objective_t(X_t, W_t, sd_t, sm_t, sls_t, slm_t, sc_t, None, cfg)
             solo_d[start:stop] = sd_t.cpu().numpy()
             solo_m[start:stop] = sm_t.cpu().numpy()
             solo_ls[start:stop] = sls_t.cpu().numpy()
             solo_lm[start:stop] = slm_t.cpu().numpy()
             solo_codes[start:stop] = sc_t.cpu().numpy()
-            solo_obj[start:stop] = sobj_t.cpu().numpy()
+            solo_obj[start:stop] = sse_Solo.cpu().numpy()
 
             elapsed = max(time.time() - t0, 1.0e-6)
             rate = stop / elapsed
@@ -1067,11 +1072,17 @@ def quantize_cuda(
             stop = min(start + batch_blocks, n_blocks)
             B_c = stop - start
             X_t = torch.from_numpy(X_np[start:stop]).to(dev)
-            W_t = (
-                torch.ones_like(X_t)
-                if W_np is None
-                else torch.from_numpy(W_np[start:stop]).to(dev)
-            )
+            if W_np is None:
+                W_t = torch.ones_like(X_t)
+            else:
+                W_raw_t = torch.from_numpy(W_np[start:stop]).to(dev)
+                if cfg.temper_imatrix:
+                    sig2_t = torch.mean(X_t * X_t, dim=1, keepdim=True)
+                    W_t = torch.sqrt(torch.clamp(W_raw_t, min=1.0e-6)) * torch.sqrt(sig2_t + X_t * X_t)
+                    wm_t = torch.mean(W_t, dim=1, keepdim=True)
+                    W_t = torch.where(wm_t > 1.0e-30, W_t / torch.clamp(wm_t, min=1.0e-30), torch.ones_like(W_t))
+                else:
+                    W_t = W_raw_t
             Xb = X_t.reshape(B_c, NS, 16)
             Wb = W_t.reshape(B_c, NS, 16)
             Sw = torch.sum(Wb, dim=2)

gguf.py

@@ -254,11 +254,17 @@ def _dequantize(raw: bytes, ti: TensorInfo) -> np.ndarray:
         return (d[:, None] * q).reshape(-1)
     if t == GGML_TYPE_Q2_K:
         blocks = np.frombuffer(raw, dtype=np.uint8).reshape(-1, 84)
+        n = blocks.shape[0]
         scales = blocks[:, :16]
         qs = blocks[:, 16:80]
         d = blocks[:, 80:82].copy().view("<f2").reshape(-1).astype(np.float32)
         dmin = blocks[:, 82:84].copy().view("<f2").reshape(-1).astype(np.float32)
-        codes = np.stack([unpack_codes(q) for q in qs])
+        codes = np.empty((n, QK_K), dtype=np.uint8)
+        for base in (0, 128):
+            q = qs[:, base // 4 : base // 4 + 32]
+            for s_lvl in range(4):
+                lo = base + 32 * s_lvl
+                codes[:, lo : lo + 32] = (q >> np.uint8(2 * s_lvl)) & np.uint8(3)
         return dequantize_params_batch(d, dmin, scales & 0x0F, scales >> 4, codes).reshape(-1)
     raise ValueError(f"unsupported source GGUF type {ti.type_name} for {ti.name}")
 

pipeline.py

@@ -52,7 +52,7 @@ from .q2k import (
     unpack_codes,
     weighted_sse,
 )
-from .seeds import build_seeds, seed_clamp_value
+from .seeds import block_clamp_values, build_seeds, seed_clamp_value
 from .sieve import block_chain_edges, sieve_measure_chain, subblock_chain_edges
 from .spectral import dc_nudge_codes, fold_weights, objective, refine_codes_hadamard
 from .viterbi import robust_temperature, viterbi_select
@@ -86,9 +86,10 @@ class QuantizedTensor:
         return np.concatenate([b.dequantize() for b in self.blocks])
 
 
-def _phase2_wls(X, W, seeds, clamp_val):
+def _phase2_wls(X, W, seeds, cfg: HPCConfig):
     n = X.shape[0]
-    clipped = np.clip(X, -clamp_val, clamp_val).astype(np.float32)
+    per_blk_clamp = block_clamp_values(X, cfg).reshape(n, 1)
+    clipped = np.clip(X, -per_blk_clamp, per_blk_clamp).astype(np.float32)
     base_d = seeds["d"].astype(np.float64).copy()
     base_m = seeds["dmin"].astype(np.float64).copy()
     ls = seeds["ls"].copy()
@@ -111,8 +112,8 @@ def _phase2_wls(X, W, seeds, clamp_val):
         safe_det = np.where(valid, det, 1.0)
         d_new = (Sbb * Sxa - Sab * Sxb) / safe_det
         m_new = (Sab * Sxa - Saa * Sxb) / safe_det
-        upd_d = valid & (d_new > 0.0) & (d_new < 4.0 * (seeds["d"] + 1.0e-10))
-        upd_m = valid & (m_new > 0.0) & (m_new < 4.0 * (seeds["dmin"] + 1.0e-10))
+        upd_d = valid & (d_new > 0.0) & np.isfinite(d_new)
+        upd_m = valid & (m_new > 0.0) & np.isfinite(m_new)
         base_d = np.where(upd_d, round_fp16(d_new.astype(np.float32)).astype(np.float64), base_d)
         base_m = np.where(upd_m, round_fp16(m_new.astype(np.float32)).astype(np.float64), base_m)
 
@@ -412,7 +413,7 @@ def _batch_extended_dm_solve_and_ulp(
     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."""
+    """Closed-form 2x2 Cramer solve + O(1) quadratic ULP polish on true weighted SSE."""
     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)
@@ -422,29 +423,6 @@ def _batch_extended_dm_solve_and_ulp(
     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)
@@ -578,7 +556,7 @@ def _refine_trajectory_algebraic(
     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
+        imp_dm = sse_dm < best_sse - 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)
@@ -591,7 +569,7 @@ def _refine_trajectory_algebraic(
 
         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
+        imp_re = sse_re < best_sse - 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)
@@ -645,7 +623,7 @@ def _refine_trajectory_algebraic(
         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
+        imp_sub = sse_try < best_sse - 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)
@@ -659,7 +637,7 @@ def _refine_trajectory_algebraic(
     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
+    imp_r = (sse_r <= best_sse * (1.0 + cfg.sse_budget)) & (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)
@@ -702,7 +680,7 @@ def _phase4_algebraic_batch(
     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(
+    d1, m1, ls1, lm1, c1, _, sse1 = _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
@@ -714,10 +692,10 @@ def _phase4_algebraic_batch(
         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(
+        d2, m2, ls2, lm2, c2, _, sse2 = _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
+        take2 = sse2 < sse1 - 1.0e-15
         fac_hits = int(np.sum(take2))
         d1 = np.where(take2, d2, d1)
         m1 = np.where(take2, m2, m1)
@@ -911,7 +889,17 @@ class Q2KQuantizer:
                 raise ValueError("importance must match the weight count")
             if np.any(imp < 0.0) or not np.all(np.isfinite(imp)):
                 raise ValueError("importance must be finite and non-negative")
-            W = imp.reshape(n_blocks, QK_K)
+            W_raw = imp.reshape(n_blocks, QK_K)
+            if cfg.temper_imatrix:
+                sig2 = np.mean(X.astype(np.float64) ** 2, axis=1, keepdims=True)
+                W = (
+                    np.sqrt(np.maximum(W_raw.astype(np.float64), 1.0e-6))
+                    * np.sqrt(sig2 + X.astype(np.float64) ** 2)
+                )
+                w_mean = np.mean(W, axis=1, keepdims=True)
+                W = np.where(w_mean > 1.0e-30, W / np.maximum(w_mean, 1.0e-30), 1.0).astype(np.float32)
+            else:
+                W = W_raw
         else:
             W = None
 
@@ -934,7 +922,7 @@ class Q2KQuantizer:
             f"Solving 5-iteration WLS across {n_blocks} superblocks",
             sigma=sigma,
         )
-        base_d, base_m = _phase2_wls(X, W, seeds, clamp_val)
+        base_d, base_m = _phase2_wls(X, W, seeds, cfg)
 
         candidate_errors = np.full((n_blocks, TOTAL_CANDIDATES), np.inf, dtype=np.float32)
         standalone: list[FactorResult] = []

seeds.py

@@ -156,16 +156,34 @@ def seed_weights(
     return np.sqrt(sigma_sq + np.minimum(block, clamp_val) ** 2)
 
 
+def block_clamp_values(blocks: np.ndarray, cfg: HPCConfig) -> np.ndarray:
+    """Per-block kurtosis-gated outlier clamp thresholds (shape (n_blocks, 1, 1))."""
+    b = np.asarray(blocks, dtype=np.float64).reshape(-1, QK_K_SEED := 256)
+    b2 = b * b
+    mean_sq = np.mean(b2, axis=1)
+    mean_fourth = np.mean(b2 * b2, axis=1)
+    kurt = np.where(mean_sq > 1.0e-30, mean_fourth / np.maximum(mean_sq * mean_sq, 1.0e-30), 3.0)
+    sigma = np.sqrt(np.maximum(mean_sq, 0.0))
+    mult = np.where(
+        kurt <= cfg.clamp_kurtosis_off,
+        1.0e9,
+        np.where(kurt <= cfg.clamp_kurtosis_soft, cfg.clamp_sigma_soft, cfg.clamp_sigma_hard),
+    )
+    clamp = np.where(mean_sq > 1.0e-30, sigma * mult, np.inf)
+    return clamp[:, None, None]
+
+
 def build_seeds(
     x: np.ndarray,
     importance: np.ndarray | None,
     cfg: HPCConfig,
-    clamp_val: float = float("inf"),
+    clamp_val: float | np.ndarray = float("inf"),
 ) -> dict[str, np.ndarray]:
     """Phase-1 seeds for all blocks: per-sub (scale,min) and shared (d,dmin)."""
     n_blocks = x.shape[0]
     blocks = np.asarray(x, dtype=np.float64).reshape(n_blocks, 16, 16)
-    clipped = np.clip(blocks, -clamp_val, clamp_val)
+    per_blk_clamp = block_clamp_values(blocks, cfg)
+    clipped = np.clip(blocks, -per_blk_clamp, per_blk_clamp)
     if importance is None:
         sigma_sq = np.mean(blocks * blocks, axis=(1, 2), keepdims=True)
         weights = np.sqrt(sigma_sq + clipped * clipped)