From 73d69c953484e0f60e571f56aed471c1b34ea57e Mon Sep 17 00:00:00 2001
From: Gregoire CATTAN <gregoire.cattan@ibm.com>
Date: Mon, 7 Sep 2026 12:09:21 +0200
Subject: [PATCH] Fix factor-of-2 in PolyOptimizer quadratic coupling encoding

Both quadratic branches of _encode() wrote half the coupling weight they
should. The comments ("factor-of-2 symmetry accounted for") show the
reasoning: the author assumed a symmetric pair contributes 2*J[i,j] to the
energy and pre-divided by two.

It does not. Solvers read J through the 1/2 * s.J.s convention -- see the
energy readout in CaSuDaSolver.solve(), and note that `field = m @ J` is
exactly the gradient of that form -- under which the (i,j) and (j,i) entries
and the leading 1/2 cancel, so a symmetric pair contributes J[i,j]
undivided.

h was already correct, so the encoded landscape was not a rescaling of the
target polynomial but a distortion of it: the linear part kept its weight
while the quadratic part lost half of its own. Minimising -0.75x - 0.75y + xy
over {0,1}^2 flattens the true optima (1,0)/(0,1) at -0.75 into a three-way
tie with (1,1) at -0.5. Across 40 random polynomials, 37 encodings failed to
be an affine image of their polynomial; doubling J fixes all 40 exactly.

p-kit already had the right convention elsewhere: TransverseFieldIsing.from_qubo
uses -Q/4 for the pair coefficient under the identical (1+s)/2 mapping. The two
converters disagreed by exactly this factor of 2.

test_quadratic_J had pinned the buggy value, having been written from the
implementation rather than the mathematics, which is how this survived. It now
asserts the derived value. The new tests check that J/h reproduce the
polynomial up to an additive constant -- comparing only argmin is unsafe, since
a mis-scaled coupling produces ties that an argmin assertion passes on an
arbitrary tie-break.

DocplexOptimizer inherits the fix.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 p_kit/library/poly.py | 16 +++++++++---
 tests/test_poly.py    | 57 +++++++++++++++++++++++++++++++++++++++++--
 2 files changed, 67 insertions(+), 6 deletions(-)

diff --git a/p_kit/library/poly.py b/p_kit/library/poly.py
index 25732ee..7db6959 100644
--- a/p_kit/library/poly.py
+++ b/p_kit/library/poly.py
@@ -79,13 +79,18 @@ def _encode(self):
             elif mono[0] == mono[1]:
                 # Quadratic x_v^2:
                 #   h contribution:  coeff * C * 2^k  per bit k
-                #   J contribution (j<k):  coeff/4 * 2^j * 2^k  (factor-of-2 symmetry accounted for)
+                #   J contribution (j<k):  coeff/2 * 2^j * 2^k
+                # The (j,k) and (k,j) terms of the double sum both survive, so
+                # the s_j*s_k coefficient is 2 * coeff/4 * 2^j * 2^k. Solvers
+                # read J through the 1/2 * s.J.s convention (see CaSuDaSolver),
+                # under which a symmetric pair contributes exactly J[j,k] --
+                # not 2*J[j,k] -- so w is that coefficient, undivided.
                 v = mono[0]
                 for k in range(n):
                     self.h[self._idx(v, k)] += coeff * C * (2 ** k)
                 for j in range(n):
                     for k in range(j + 1, n):
-                        w = coeff / 4 * (2 ** j) * (2 ** k)
+                        w = coeff / 2 * (2 ** j) * (2 ** k)
                         ij, ik = self._idx(v, j), self._idx(v, k)
                         self.J[ij, ik] += w
                         self.J[ik, ij] += w
@@ -93,7 +98,10 @@ def _encode(self):
             else:
                 # Cross-term x_u * x_v (u != v):
                 #   h contribution:  coeff * C/2 * 2^k  per bit k, for both variables
-                #   J contribution:  coeff/8 * 2^j * 2^k  (factor-of-2 symmetry accounted for)
+                #   J contribution:  coeff/4 * 2^j * 2^k
+                # b_j*b_k expands to (1 + s_j + s_k + s_j*s_k)/4, so the pair
+                # coefficient is coeff/4 * 2^j * 2^k. See the x_v^2 branch for
+                # why that goes into J undivided.
                 u, v = mono
                 for k in range(n):
                     wk = 2 ** k
@@ -101,7 +109,7 @@ def _encode(self):
                     self.h[self._idx(v, k)] += coeff * C / 2 * wk
                 for j in range(n):
                     for k in range(n):
-                        w = coeff / 8 * (2 ** j) * (2 ** k)
+                        w = coeff / 4 * (2 ** j) * (2 ** k)
                         ij, ik = self._idx(u, j), self._idx(v, k)
                         self.J[ij, ik] += w
                         self.J[ik, ij] += w
diff --git a/tests/test_poly.py b/tests/test_poly.py
index 55f3169..110f636 100644
--- a/tests/test_poly.py
+++ b/tests/test_poly.py
@@ -44,9 +44,62 @@ def test_quadratic_h():
 
 
 def test_quadratic_J():
-    # 3*x^2 with n_bits=2: J[0,1] = 3/4 * 2^0 * 2^1 = 1.5
+    # 3*x^2 with n_bits=2: J[0,1] = 3/2 * 2^0 * 2^1 = 3.0
     circuit = PolyOptimizer({('x', 'x'): 3}, ['x'], n_bits=2, minimize=False)
-    assert np.isclose(circuit.J[0, 1], 1.5)
+    assert np.isclose(circuit.J[0, 1], 3.0)
+
+
+def _landscape(circuit, coeffs, variables, n_bits):
+    """Yield (f(x), F(s)) over every p-bit state.
+
+    F is the quantity the solvers actually see: h.s + 1/2 s.J.s, matching
+    CaSuDaSolver's energy readout and the gradient its update rule ascends.
+    """
+    import itertools
+    h = np.asarray(circuit.h).reshape(-1)
+    for bits in itertools.product([-1, 1], repeat=len(variables) * n_bits):
+        s = np.array(bits, dtype=float)
+        F = s @ h + 0.5 * (s @ circuit.J @ s)
+        vals = circuit.decode(s)
+        f = 0.0
+        for mono, c in coeffs.items():
+            if len(mono) == 1:
+                f += c * vals[mono[0]]
+            elif len(mono) == 2:
+                f += c * vals[mono[0]] * vals[mono[1]]
+        yield f, F
+
+
+@pytest.mark.parametrize("coeffs,variables,n_bits", [
+    ({('x',): -0.75, ('y',): -0.75, ('x', 'y'): 1}, ['x', 'y'], 1),
+    ({('x', 'x'): 1, ('x',): -3}, ['x'], 2),
+    ({('x', 'x'): 2, ('y', 'y'): -1, ('x', 'y'): 3, ('x',): 1}, ['x', 'y'], 2),
+])
+def test_encoding_is_faithful(coeffs, variables, n_bits):
+    """J/h must reproduce the polynomial up to an additive constant.
+
+    Comparing only argmin is not enough: a mis-scaled coupling flattens the
+    landscape into ties that an argmin check passes by arbitrary tie-break.
+    """
+    circuit = PolyOptimizer(coeffs, variables, n_bits=n_bits, minimize=False)
+    resid = np.array([F - f for f, F in _landscape(circuit, coeffs, variables, n_bits)])
+    assert np.allclose(resid, resid[0])
+
+
+def test_cross_term_ranks_optimum_strictly():
+    # min of -0.75x - 0.75y + xy over {0,1}^2 is -0.75 at (1,0) and (0,1);
+    # (1,1) scores -0.5 and must stay strictly worse. A halved xy coupling
+    # ties all three.
+    coeffs = {('x',): -0.75, ('y',): -0.75, ('x', 'y'): 1}
+    circuit = PolyOptimizer(coeffs, ['x', 'y'], n_bits=1, minimize=True)
+    # minimize=True negates, and the solvers ascend F, so the optimum is argmax F
+    best = max(F for _, F in _landscape(circuit, coeffs, ['x', 'y'], 1))
+    ranked = sorted(
+        (F, f) for f, F in _landscape(circuit, coeffs, ['x', 'y'], 1)
+    )
+    assert np.isclose(ranked[-1][1], -0.75)
+    assert not np.isclose(ranked[-1][0], ranked[-3][0])
+    assert np.isclose(best, ranked[-1][0])
 
 
 def test_minimize_negates():