From fd5720833a73dac222d0663eb753b9dc3bb32bfd Mon Sep 17 00:00:00 2001 From: Gregoire CATTAN Date: Sat, 23 May 2026 08:44:48 +0200 Subject: [PATCH] add paoa --- README.md | 25 ++++++++++ examples/docplex_poly.py | 7 ++- examples/quantum_tfim.py | 69 ++++++++++++++++++++++++++++ p_kit/library/__init__.py | 5 +- p_kit/library/quantum.py | 96 +++++++++++++++++++++++++++++++++++++++ tests/test_quantum.py | 82 +++++++++++++++++++++++++++++++++ 6 files changed, 281 insertions(+), 3 deletions(-) create mode 100644 examples/quantum_tfim.py create mode 100644 p_kit/library/quantum.py create mode 100644 tests/test_quantum.py diff --git a/README.md b/README.md index a9282ba..62c8a22 100644 --- a/README.md +++ b/README.md @@ -123,4 +123,29 @@ class LogicCircuit: self.and_gate.output.connect(self.or_gate.input1) ``` +### Higher-level Library + +`p_kit.library` provides ready-made circuits for common problem types. + +| Class | Purpose | Extra p-bits mean… | +|---|---|---| +| `TSP` | Travelling Salesman Problem | one p-bit per (city, order) pair | +| `PolyOptimizer` | Classical polynomial minimisation (QUBO) | binary bits for integer precision | +| `TransverseFieldIsing` | Quantum-inspired annealing (TFIM emulation) | Trotter replicas in imaginary time | + +Both `PolyOptimizer` and `TransverseFieldIsing` can target the same QUBO/Ising problems, but they correspond to different parts of QAOA: + +- **`PolyOptimizer`** encodes only the QAOA *cost Hamiltonian* H_C = -Σ J_ij σ^z_i σ^z_j. It uses pure Boltzmann sampling to find the minimum — equivalent to QAOA with zero mixer strength (the classical baseline). + +- **`TransverseFieldIsing`** emulates the full QAOA *dynamics*: the cost Hamiltonian (intra-replica couplings) plus the mixer Hamiltonian H_B = -Γ Σ σ^x_i (inter-replica Trotter couplings). The `gamma` parameter is the QAOA mixer strength and `n_replicas` corresponds to the number of QAOA layers. This is what p-kit calls **PAOA** — a probabilistic emulation of QAOA for stoquastic Hamiltonians, running at room temperature on classical hardware (see [Camsari et al., 2019](https://arxiv.org/abs/1809.04028), Fig. 9). + +To use `TransverseFieldIsing` with a QUBO problem: + +```python +from p_kit.library.quantum import TransverseFieldIsing + +# Q is your (n, n) QUBO matrix +circuit = TransverseFieldIsing.from_qubo(Q, gamma=0.5, beta=2.0, n_replicas=20) +``` + For more information, visit our [documentation](https://github.com/IBM/p-kit/wiki). diff --git a/examples/docplex_poly.py b/examples/docplex_poly.py index cf66738..cfcb17f 100644 --- a/examples/docplex_poly.py +++ b/examples/docplex_poly.py @@ -1,5 +1,8 @@ -from docplex.mp.model import Model -from p_kit.library.docplex_optimizer import DocplexOptimizer +try: + from docplex.mp.model import Model + from p_kit.library.docplex_optimizer import DocplexOptimizer +except ImportError: + raise SystemExit(0) # skip gracefully if docplex is not installed from p_kit.solver.csd_solver import CaSuDaSolver from p_kit.solver.annealing import constant, execute from p_kit.visualization.histplot import histplot diff --git a/examples/quantum_tfim.py b/examples/quantum_tfim.py new file mode 100644 index 0000000..edba23d --- /dev/null +++ b/examples/quantum_tfim.py @@ -0,0 +1,69 @@ +# Transverse-Field Ising Model (TFIM) emulated with p-bits +# via the Suzuki-Trotter path integral decomposition. +# +# Ref: Camsari et al., arXiv:1809.04028, Fig. 9 +# +# A stoquastic n-qubit TFIM is mapped to a p-circuit with +# n * n_replicas p-bits. Intra-replica couplings reproduce +# the Ising Hamiltonian; inter-replica couplings encode the +# transverse field via the Trotter coupling K. + +import matplotlib.pyplot as plt +import numpy as np + +from p_kit.library.quantum import TransverseFieldIsing +from p_kit.solver.annealing import linear, execute +from p_kit.solver.csd_solver import CaSuDaSolver + +# --- Problem definition -------------------------------------------------- +# 3-qubit ferromagnetic TFIM: J_ij = 1 (spins want to align) +n = 3 +J_q = np.array([ + [0., 1., 1.], + [1., 0., 1.], + [1., 1., 0.], +]) +h_q = np.zeros(n) + +# Transverse field strength (quantum fluctuations). +# gamma=0 → classical Ising (all replicas lock together) +# gamma>>0 → strong quantum fluctuations +gamma = 0.5 +beta = 2.0 # inverse temperature +n_replicas = 20 + +# --- Build p-circuit from Trotter decomposition -------------------------- +circuit = TransverseFieldIsing(J_q, h_q, gamma=gamma, beta=beta, n_replicas=n_replicas) +print(f"Qubits: {n} | Replicas: {n_replicas} | Total p-bits: {circuit.n_pbits}") +print(f"Trotter coupling K = {-0.5 * np.log(np.tanh(beta * gamma / n_replicas)):.4f}") + +# --- Solve with annealing ------------------------------------------------ +solver = CaSuDaSolver(Nt=int(5e4), dt=0.1, i0=0, seed=42) +n_shots = 100 +samples = execute(solver, circuit, linear, n_shots=n_shots, n_last_samples=50, n_jobs=-1) + +# Recover qubit magnetizations by averaging over replicas +# samples shape: (n_shots * n_last_samples, n * n_replicas) +qubit_m = samples.reshape(-1, n_replicas, n) # (samples, replicas, qubits) +qubit_m_avg = qubit_m.mean(axis=1) # (samples, qubits): per-qubit mean over replicas + +# --- Plot ---------------------------------------------------------------- +fig, axes = plt.subplots(1, 2, figsize=(12, 4)) +fig.suptitle(f"TFIM emulated via Suzuki-Trotter (n={n}, R={n_replicas}, γ={gamma}, β={beta})") + +# Left: histogram of qubit-0 averaged magnetization +axes[0].hist(qubit_m_avg[:, 0], bins=30, density=True, color="steelblue", edgecolor="white") +axes[0].set_xlabel("⟨m₀⟩ (averaged over replicas)") +axes[0].set_ylabel("Probability density") +axes[0].set_title("Qubit-0 magnetization distribution") + +# Right: inter-qubit correlation +corr = np.mean(qubit_m_avg[:, :, None] * qubit_m_avg[:, None, :], axis=0) +im = axes[1].imshow(corr, vmin=-1, vmax=1, cmap="RdBu_r") +axes[1].set_title("Qubit-qubit correlations ⟨mᵢ mⱼ⟩") +axes[1].set_xticks(range(n)) +axes[1].set_yticks(range(n)) +plt.colorbar(im, ax=axes[1]) + +plt.tight_layout() +plt.show() diff --git a/p_kit/library/__init__.py b/p_kit/library/__init__.py index 15e8e39..770d02d 100644 --- a/p_kit/library/__init__.py +++ b/p_kit/library/__init__.py @@ -1,3 +1,6 @@ from . import tsp +from . import poly +from . import docplex_optimizer +from . import quantum -__all__ = ["tsp"] +__all__ = ["tsp", "poly", "docplex_optimizer", "quantum"] diff --git a/p_kit/library/quantum.py b/p_kit/library/quantum.py new file mode 100644 index 0000000..35a7742 --- /dev/null +++ b/p_kit/library/quantum.py @@ -0,0 +1,96 @@ +import numpy as np +from p_kit.psl.p_circuit import PCircuit + + +class TransverseFieldIsing(PCircuit): + """Maps a stoquastic transverse-field Ising model to a p-bit circuit + via the Suzuki-Trotter (path integral) decomposition. + + An n-qubit quantum system becomes an n * n_replicas p-bit circuit. + P-bit (i, tau) maps to index tau * n + i (replica-major ordering). + + This is the probabilistic analog of QAOA (PAOA): intra-replica couplings + encode the QAOA cost Hamiltonian H_C = -sum_ij J_ij sz_i sz_j, while + inter-replica Trotter couplings encode the QAOA mixer H_B = -gamma sum_i sx_i. + ``gamma`` is the mixer strength and ``n_replicas`` the number of QAOA layers. + Use ``PolyOptimizer`` instead for the classical baseline (gamma=0 equivalent). + + Parameters + ---------- + j_q : np.ndarray, shape (n, n) + Symmetric Ising coupling matrix of the quantum system (QAOA cost H_C). + h_q : np.ndarray, shape (n,) + Longitudinal field (bias) of the quantum system. + gamma : float + Transverse field / QAOA mixer strength. Use 0 for a purely classical Ising model. + beta : float, default 1.0 + Inverse temperature (1/T). + n_replicas : int, default 10 + Number of Trotter replicas (QAOA layers). Higher values give a better approximation. + + Notes + ----- + Ref: Camsari et al., arXiv:1809.04028, Fig. 9 + + .. versionadded:: 0.0.1 + """ + + def __init__(self, j_q, h_q, gamma, beta=1.0, n_replicas=10): + n = len(h_q) + R = n_replicas + super().__init__(n * R) + + J = np.zeros((n * R, n * R)) + h_vec = np.zeros(n * R) + + # Intra-replica: spatial couplings and bias scaled by beta/R + scale = beta / R + for tau in range(R): + s = tau * n + J[s:s + n, s:s + n] += scale * j_q + h_vec[s:s + n] += scale * h_q + + # Inter-replica: Trotter coupling K = -(1/2) ln(tanh(beta*gamma/R)) + # K > 0 ferromagnetically couples same qubit across adjacent replicas. + # gamma=0 means classical limit (no inter-replica coupling). + if gamma > 0: + K = -0.5 * np.log(np.tanh(beta * gamma / R)) + seen = set() + for tau in range(R): + tau_next = (tau + 1) % R + pair = (min(tau, tau_next), max(tau, tau_next)) + if pair in seen: + continue + seen.add(pair) + idx = np.arange(n) + J[tau * n + idx, tau_next * n + idx] += K + J[tau_next * n + idx, tau * n + idx] += K + + np.fill_diagonal(J, 0) + self.J = J + self.h = h_vec + + @classmethod + def from_qubo(cls, qubo, gamma, beta=1.0, n_replicas=10): + """Build a TransverseFieldIsing circuit from a QUBO matrix. + + Minimises x^T Q x with x ∈ {0, 1} by mapping to the Ising + spin variables s ∈ {-1, +1} via x_i = (1 + s_i) / 2. + + Parameters + ---------- + qubo : np.ndarray, shape (n, n) + QUBO matrix (upper or full triangle). + gamma : float + Transverse field strength. + beta : float, default 1.0 + Inverse temperature (1/T). + n_replicas : int, default 10 + Number of Trotter replicas. + """ + qubo = np.asarray(qubo, dtype=float) + 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 + return cls(j_q, h_q, gamma, beta, n_replicas) diff --git a/tests/test_quantum.py b/tests/test_quantum.py new file mode 100644 index 0000000..e42d57b --- /dev/null +++ b/tests/test_quantum.py @@ -0,0 +1,82 @@ +import numpy as np +import pytest +from p_kit.library.quantum import TransverseFieldIsing +from p_kit.solver import CaSuDaSolver + + +def _make(n=3, R=4, gamma=0.5, beta=1.0, seed=0): + rng = np.random.default_rng(seed) + J_q = rng.uniform(-1, 1, (n, n)) + J_q = (J_q + J_q.T) / 2 + np.fill_diagonal(J_q, 0) + h_q = rng.uniform(-1, 1, n) + return TransverseFieldIsing(J_q, h_q, gamma=gamma, beta=beta, n_replicas=R), n, R + + +def test_shape(): + c, n, R = _make() + assert c.n_pbits == n * R + assert c.J.shape == (n * R, n * R) + assert c.h.shape == (n * R,) + + +def test_symmetry(): + c, _, _ = _make() + np.testing.assert_array_almost_equal(c.J, c.J.T) + + +def test_zero_diagonal(): + c, n, R = _make() + np.testing.assert_array_equal(np.diag(c.J), np.zeros(n * R)) + + +def test_intra_replica_scaling(): + n, R, beta = 2, 3, 2.0 + J_q = np.array([[0.0, 1.0], [1.0, 0.0]]) + h_q = np.array([0.5, -0.5]) + c = TransverseFieldIsing(J_q, h_q, gamma=0.0, beta=beta, n_replicas=R) + scale = beta / R + for tau in range(R): + s = tau * n + np.testing.assert_array_almost_equal(c.J[s:s + n, s:s + n], scale * J_q) + np.testing.assert_array_almost_equal(c.h[s:s + n], scale * h_q) + + +def test_no_inter_replica_when_gamma_zero(): + n, R = 2, 3 + J_q = np.array([[0.0, 1.0], [1.0, 0.0]]) + h_q = np.zeros(n) + c = TransverseFieldIsing(J_q, h_q, gamma=0.0, n_replicas=R) + for tau in range(R): + tau_next = (tau + 1) % R + assert c.J[tau * n, tau_next * n] == 0.0 + + +def test_inter_replica_coupling_value(): + n, R = 2, 4 + beta, gamma = 1.0, 0.5 + J_q = np.zeros((n, n)) + h_q = np.zeros(n) + c = TransverseFieldIsing(J_q, h_q, gamma=gamma, beta=beta, n_replicas=R) + K = -0.5 * np.log(np.tanh(beta * gamma / R)) + # qubit 0: replica 0 connected to replica 1 at index [0, n] + assert c.J[0, n] == pytest.approx(K) + assert c.J[n, 0] == pytest.approx(K) + + +def test_r2_no_double_counting(): + """For R=2 the ring has one unique temporal bond per qubit — not two.""" + n, R = 2, 2 + beta, gamma = 1.0, 0.5 + J_q = np.zeros((n, n)) + h_q = np.zeros(n) + c = TransverseFieldIsing(J_q, h_q, gamma=gamma, beta=beta, n_replicas=R) + K = -0.5 * np.log(np.tanh(beta * gamma / R)) + assert c.J[0, n] == pytest.approx(K) # not 2*K + + +def test_solve(): + c, _, _ = _make(n=2, R=4) + 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