From 2d49936e879a673554818054cfc51f46841795b7 Mon Sep 17 00:00:00 2001
From: Gregoire CATTAN <gregoire.cattan@ibm.com>
Date: Sat, 5 Sep 2026 23:11:35 +0200
Subject: [PATCH] Fix missing diagonal term in TransverseFieldIsing.from_qubo
 bias

The QUBO-to-Ising substitution x=(1+s)/2 requires the diagonal of the
symmetrized QUBO matrix to appear in the linear bias h_q, not just in
the row sums. Without it, h_q was systematically wrong and the solver
would converge confidently to incorrect optima. Verified the fix by
brute-force affine-consistency checks and by confirming the solver now
recovers the true optimum end-to-end.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
---
 p_kit/library/quantum.py |  2 +-
 tests/test_quantum.py    | 56 ++++++++++++++++++++++++++++++++++++++++
 2 files changed, 57 insertions(+), 1 deletion(-)

diff --git a/p_kit/library/quantum.py b/p_kit/library/quantum.py
index 35a7742..4558998 100644
--- a/p_kit/library/quantum.py
+++ b/p_kit/library/quantum.py
@@ -92,5 +92,5 @@ def from_qubo(cls, qubo, gamma, beta=1.0, n_replicas=10):
         qubo_sym = qubo + qubo.T - np.diag(np.diag(qubo))   # symmetrise, count diagonal once
         j_q = -qubo_sym / 4
         np.fill_diagonal(j_q, 0)
-        h_q = -np.sum(qubo_sym, axis=1) / 4
+        h_q = -(np.sum(qubo_sym, axis=1) + np.diag(qubo_sym)) / 4
         return cls(j_q, h_q, gamma, beta, n_replicas)
diff --git a/tests/test_quantum.py b/tests/test_quantum.py
index e42d57b..3576900 100644
--- a/tests/test_quantum.py
+++ b/tests/test_quantum.py
@@ -1,7 +1,10 @@
+import itertools
+
 import numpy as np
 import pytest
 from p_kit.library.quantum import TransverseFieldIsing
 from p_kit.solver import CaSuDaSolver
+from p_kit.solver.annealing import constant
 
 
 def _make(n=3, R=4, gamma=0.5, beta=1.0, seed=0):
@@ -80,3 +83,56 @@ def test_solve():
     solver = CaSuDaSolver(Nt=100, dt=0.1667, i0=0.8, seed=42)
     _, all_m, _ = solver.solve(c)
     assert all_m.shape[-1] == c.n_pbits
+
+
+def _brute_force_qubo_min(qubo):
+    n = qubo.shape[0]
+    best = None
+    for bits in itertools.product([0, 1], repeat=n):
+        x = np.array(bits, dtype=float)
+        val = x @ qubo @ x
+        if best is None or val < best:
+            best = val
+    return best
+
+
+@pytest.mark.parametrize("qubo", [
+    np.array([[1.0, -2.0], [0.0, 1.0]]),
+    np.random.default_rng(0).standard_normal((4, 4)),
+])
+def test_from_qubo_h_matches_energy(qubo):
+    """h_q/j_q must encode -x^T Q x up to an additive constant: the
+    "maximize" energy F(s) = h@s + 0.5*s@J@s (the same form used by
+    CaSuDaSolver/GibbsSolver) has to differ from -x^T Q x by a constant
+    across every spin configuration, not just at the optimum.
+    """
+    n = qubo.shape[0]
+    circuit = TransverseFieldIsing.from_qubo(qubo, gamma=0.0, n_replicas=1)
+    J, h = circuit.J, circuit.h
+
+    offsets = []
+    for bits in itertools.product([-1, 1], repeat=n):
+        s = np.array(bits, dtype=float)
+        x = (1 + s) / 2
+        F = h @ s + 0.5 * s @ J @ s
+        offsets.append(F + x @ qubo @ x)
+
+    np.testing.assert_allclose(offsets, offsets[0])
+
+
+def test_from_qubo_solver_finds_optimum():
+    qubo = np.array([
+        [1.0, -2.0, 0.0],
+        [0.0, 1.0, -2.0],
+        [0.0, 0.0, 1.0],
+    ])
+    true_min = _brute_force_qubo_min(qubo)
+
+    circuit = TransverseFieldIsing.from_qubo(qubo, gamma=0.0, n_replicas=1)
+    solver = CaSuDaSolver(Nt=2000, dt=0.1667, i0=2.0, seed=7)
+    _, all_m, _ = solver.solve(circuit, annealing_func=constant)
+
+    x_samples = (all_m[-500:] + 1) / 2
+    values = np.einsum("ij,jk,ik->i", x_samples, qubo, x_samples)
+    assert np.isclose(values.min(), true_min)
+    assert np.mean(np.isclose(values, true_min)) > 0.3