bqs 0.2.0

Better Quantum Software

bqs — Better Quantum Software

bqs is a Python framework for building scalable, parallelized Variational Quantum Algorithms (VQAs), with a particular focus on hardware-efficient implementations in the NISQ era.

The design and methodology implemented in this package are described in
Parallel Circuit Implementation of Variational Quantum Algorithms
(preprint: arXiv:2304.03037).

---

Installation

pip install bqs

---

Quick start

import numpy as np
from bqs import QAOA, COBYLA

# Define a symmetric QUBO matrix
Q = np.array([[0.0, 0.5, 0.5],
              [0.5, 0.0, 0.5],
              [0.5, 0.5, 0.0]])

alg = QAOA(Q, p=2)
alg.set_optimizer(COBYLA)
alg.optimize(num_samples_training=100, num_iterations=50, opt_arguments={"display": False})

samples = alg.circuit.sample(num_samples=200)   # (200, 3) array of ±1 spins
energy  = alg.objective_function.evaluate_samples(samples)

---

Background

Problem formulation

bqs works with optimization problems expressed as QUBOs (Quadratic Unconstrained Binary Optimization):

$$E(\mathbf{x}) = \mathbf{x}^T Q \mathbf{x} + \text{offset}, \qquad x_i \in {0, 1}$$

Internally, all algorithms convert the QUBO to its equivalent Ising Hamiltonian via the substitution $x_i = (1 - s_i)/2$, $s_i \in {-1, +1}$:

$$H(\mathbf{s}) = \mathbf{s}^T J \mathbf{s} + \mathbf{h}^T \mathbf{s} + \text{offset}$$

where $J = \frac{1}{2}\text{triu}(Q, k=1)$ and $h_i = \frac{1}{2}\sum_j Q_{ij}$.

Users can supply either a QUBO matrix $Q$ to the constructor, or the Ising triple $(J, h, \text{offset})$ directly via the from_ising classmethod.

QAOA

The Quantum Approximate Optimization Algorithm (QAOA) is a variational hybrid algorithm. It applies $p$ alternating layers of a cost unitary $e^{-i\gamma_l H}$ and a mixer unitary $e^{-i\beta_l \sum_i X_i}$ to an initial uniform superposition:

$$U(\boldsymbol{\gamma}, \boldsymbol{\beta}) = \prod_{l=1}^{p} e^{-i\beta_l H_\text{mix}} e^{-i\gamma_l H} \cdot H^{\otimes n}|0\rangle^{\otimes n}$$

The $2p$ real parameters $(\boldsymbol{\gamma}, \boldsymbol{\beta})$ are tuned by a classical optimizer to minimize the expected energy $\langle H \rangle$ estimated from circuit samples.

Problem decomposition (pQAOA)

For large problems, the full QUBO may require more qubits than available hardware supports. bqs implements a parallel circuit decomposition strategy: the global problem $P$ is split into $k$ smaller sub-problems ${SP_i}$, each solved on a smaller QAOA circuit. Samples from the sub-circuits are glued back together and evaluated against the global cost function.

This approach preserves solution quality while drastically reducing the per-circuit qubit count.

---

Algorithm classes

All algorithm classes live in bqs and share a common interface.

QAOA

Vanilla QAOA on the full problem.

QAOA(
    Q,                        # (n, n) symmetric QUBO matrix
    offset=0.0,               # constant energy offset
    p=1,                      # number of QAOA layers
    hyperparameters=None,     # (2*p,) array; random in [0, 2π) if None
    qubit_order=None,         # (n,) int array mapping variable i → physical qubit index
    quantum_objects=None,     # AbstractQuantumObjects backend; CirqQuantumObjects if None
)

Named constructor — bypass the internal QUBO→Ising conversion when the Ising form is already available:

QAOA.from_ising(
    J,                        # (n, n) upper-triangular Ising coupling matrix
    h,                        # (n,) Ising local-field vector
    offset=0.0,
    p=1,
    hyperparameters=None,
    qubit_order=None,
    quantum_objects=None,
)

---

pQAOA

Parallel QAOA: each sub-problem slice runs its own QAOA circuit with independent $2p$ hyperparameters. The full problem must be decomposable as a concatenation of the slice variables.

pQAOA(
    Q,                        # (n_full, n_full) full-model QUBO
    slice_Qs,                 # list of per-slice QUBO matrices
    offset=0.0,
    slice_offsets=None,       # per-slice offsets; all 0 if None
    p=1,
    hyperparameters=None,     # (n_slices, 2*p) array; random if None
    qubit_orders=None,        # list of per-slice qubit order arrays
    cartesian=True,           # True → Cartesian-product gluing; False → zip gluing
    quantum_objects=None,
)

Named constructor:

pQAOA.from_ising(
    J, h, offset=0.0,
    slice_isings,             # list of (J_k, h_k) or (J_k, h_k, offset_k) tuples
    p=1,
    hyperparameters=None,
    qubit_orders=None,
    cartesian=True,
    quantum_objects=None,
)

Sample gluing:
With cartesian=True, if each slice circuit produces $m$ samples, the glued output has $m^k$ rows — every combination of slice outcomes. With cartesian=False, slices are zipped row-by-row, producing $m$ rows (each slice must produce the same number of samples).

---

pQAOASingleParameters

Same as pQAOA but all slices share a single set of $2p$ hyperparameters.

pQAOASingleParameters(
    Q, slice_Qs,
    offset=0.0, slice_offsets=None,
    p=1,
    hyperparameters=None,     # (2*p,) array; random if None
    qubit_orders=None,
    cartesian=True,
    quantum_objects=None,
)

Named constructor:

pQAOASingleParameters.from_ising(
    J, h, offset=0.0,
    slice_isings,
    p=1,
    hyperparameters=None,
    qubit_orders=None,
    cartesian=True,
    quantum_objects=None,
)

---

SingleSliceQAOA

A single shared circuit is run repeatedly for each copy of the slice. A variable map then remaps each copy's outputs back onto the full-model columns.

This is most useful when the global problem has a repetitive structure (e.g., identical sub-blocks), so one circuit suffices for all slices.

SingleSliceQAOA(
    Q,                        # (n_full, n_full) full-model QUBO
    slice_Q,                  # (n_slice, n_slice) QUBO for the repeated circuit
    var_map,                  # int array (n_slice, num_slices):
                              #   var_map[i, k] = full-model column for slice var i in copy k
    offset=0.0,
    slice_offset=0.0,
    p=1,
    hyperparameters=None,
    qubit_order=None,
    quantum_objects=None,
)

Named constructor:

SingleSliceQAOA.from_ising(
    J, h, offset=0.0,
    slice_J, slice_h, slice_offset=0.0,
    var_map,
    p=1,
    hyperparameters=None,
    qubit_order=None,
    quantum_objects=None,
)

Example var_map — 6-variable problem with 2-variable slice repeated 3 times:

var_map = np.array([[0, 2, 4],   # slice var 0 → full-model columns 0, 2, 4
                    [1, 3, 5]])  # slice var 1 → full-model columns 1, 3, 5

---

Running the optimization

All algorithm classes share the same three-step interface:

# 1. Instantiate
alg = QAOA(Q, p=2)

# 2. Attach an optimizer
alg.set_optimizer(COBYLA)

# 3. Optimize
alg.optimize(
    num_samples_training=100,     # samples per cost-function evaluation
    num_iterations=200,           # max optimizer iterations
    opt_arguments={"display": False},
    # circuit_options={"parallel": True},  # run slice circuits in parallel (pQAOA)
)

After optimization the best hyperparameters are stored in alg.hyperparameters and alg.circuit.hyperparameters.

Sampling results

samples = alg.circuit.sample(num_samples=500)
# → np.ndarray (500, n) of ±1 spin values

samples_binary = alg.circuit.sample(num_samples=500, original_basis=True)
# → np.ndarray (500, n) of {0, 1} binary values

energy = alg.objective_function.evaluate_samples(samples)
# → float: mean Ising energy over the sample batch

---

Complete example

import numpy as np
from bqs import QAOA, pQAOA, pQAOASingleParameters, SingleSliceQAOA, COBYLA
from bqs.utils.bqm import qubo_to_ising

# 6-variable chain: (0,0)-(0,1), (0,0)-(1,0), (1,0)-(1,1), (1,1)-(2,0), (2,0)-(2,1)
N = 6
Q_full = np.zeros((N, N))
for i, j in [(0,1),(0,2),(2,3),(3,4),(4,5)]:
    Q_full[i,j] = Q_full[j,i] = 0.5

# Sub-problem slices
Q_slice1 = np.zeros((2, 2)); Q_slice1[0,1] = Q_slice1[1,0] = 0.5
Q_slice2 = np.zeros((4, 4))
for i, j in [(0,1),(1,2),(2,3)]:
    Q_slice2[i,j] = Q_slice2[j,i] = 0.5

def run(alg, num_samples=50, num_iter=20):
    alg.set_optimizer(COBYLA)
    alg.optimize(num_samples_training=num_samples, num_iterations=num_iter,
                 opt_arguments={"display": False})
    return alg.objective_function.evaluate_samples(
        alg.circuit.sample(num_samples=num_samples))

# Vanilla QAOA
print("QAOA energy:", run(QAOA(Q_full, p=2)))

# pQAOA: per-slice hyperparameters, Cartesian-product gluing
print("pQAOA energy:", run(pQAOA(Q_full, [Q_slice1, Q_slice2])))

# pQAOA: zip gluing
print("pQAOA (zip) energy:", run(pQAOA(Q_full, [Q_slice1, Q_slice2], cartesian=False)))

# Shared hyperparameters across slices
print("pQAOASingleParameters energy:", run(pQAOASingleParameters(Q_full, [Q_slice1, Q_slice2])))

# Single shared circuit for all slices
Q_ss = np.zeros((2,2)); Q_ss[0,1] = Q_ss[1,0] = 0.5
var_map = np.array([[0,2,4],[1,3,5]])
print("SingleSliceQAOA energy:", run(SingleSliceQAOA(Q_full, Q_ss, var_map)))

# Using Ising inputs directly (skips the internal QUBO→Ising conversion)
J, h, off = qubo_to_ising(Q_full)
print("QAOA (from_ising):", run(QAOA.from_ising(J, h, off, p=2)))

---

Quantum backend

The quantum_objects parameter accepts any subclass of AbstractQuantumObjects. This allows you to plug in a different quantum simulator or hardware SDK without changing any algorithm code.

The default backend is CirqQuantumObjects, which uses Google Cirq's statevector simulator.

Using the default backend explicitly

from bqs import QAOA, CirqQuantumObjects

alg = QAOA(Q, quantum_objects=CirqQuantumObjects())

Writing a custom backend

Subclass AbstractQuantumObjects and implement the required interface:

from bqs import AbstractQuantumObjects
import numpy as np

class MyBackend(AbstractQuantumObjects):

    def qubit(self, index, **kwargs):
        """Return the backend's qubit object for physical qubit `index`."""
        ...

    def circuit(self, n_qubits):
        """Create and return an empty circuit for `n_qubits` qubits."""
        ...

    def append(self, circuit, operation):
        """Append `operation` to `circuit` in-place."""
        ...

    def sample(self, circuit, qubit_name_to_object, num_samples, **kwargs) -> np.ndarray:
        """
        Run `circuit` for `num_samples` shots.

        Returns
        -------
        np.ndarray of shape (num_samples, n) with values in {-1, +1}.
        qubit_name_to_object maps integer variable index i → qubit object.
        """
        ...

    def H(self, qubit):    ...   # Hadamard gate
    def X(self, qubit):    ...   # Pauli-X gate
    def Y(self, qubit):    ...   # Pauli-Y gate
    def Z(self, qubit):    ...   # Pauli-Z gate
    def CNOT(self, q1, q2):  ... # CNOT gate
    def rx(self, qubit, angle): ...   # Rx(angle) rotation
    def ry(self, qubit, angle): ...   # Ry(angle) rotation
    def rz(self, qubit, angle): ...   # Rz(angle) rotation

alg = QAOA(Q, quantum_objects=MyBackend())

---

Optimizers

Three optimizers are included out of the box:

Class	Description
COBYLA	Gradient-free; wraps scipy.optimize.minimize with method 'COBYLA'. Good default choice.
SPSA	Simultaneous Perturbation Stochastic Approximation. Useful for noisy landscapes.
MonteCarlo	Random hyperparameter search; keeps the best set found. Useful as a baseline.

All are passed as classes (not instances) to set_optimizer:

alg.set_optimizer(COBYLA)
alg.set_optimizer(SPSA)
alg.set_optimizer(MonteCarlo)

Writing a custom optimizer

Subclass Optimizer and implement _run:

from bqs.utils.vqa_utils import CostFunction, VariationalCircuit
from bqs.utils.optimizers import Optimizer
from scipy.optimize import minimize, OptimizeResult

class NelderMead(Optimizer):

    def __init__(self, cost_function: CostFunction, circuit: VariationalCircuit):
        super().__init__(cost_function, circuit)

    def _run(self, objective_function, num_iterations: int, display: bool = False, **kwargs):
        result = minimize(
            objective_function,
            self.hyperparameters,
            method="Nelder-Mead",
            options={"maxiter": num_iterations, "disp": display},
            **kwargs,
        )
        return result   # must be a scipy OptimizeResult (has .x and .fun)

alg.set_optimizer(NelderMead)

---

BQM utilities

from bqs.utils.bqm import qubo_to_ising, ising_to_qubo, qubo_energy, ising_energy

# QUBO → Ising
J, h, offset = qubo_to_ising(Q)
# J: (n, n) upper-triangular; h: (n,) local fields; offset: float

# Ising → QUBO
Q2, new_offset = ising_to_qubo(J, h, offset)

# Evaluate a single sample
x = np.array([0, 1, 1, 0, 1, 0])          # binary {0,1}
s = np.array([-1, 1, 1, -1, 1, -1])       # spin {-1,+1}
e_qubo  = qubo_energy(Q, 0.0, x)
e_ising = ising_energy(J, h, offset, s)   # equals e_qubo

---

Architecture overview

BaseVQA (abstract)
├── BaseQAOA (abstract)
│   ├── QAOA
│   ├── pQAOA
│   │   └── pQAOASingleParameters
│   └── SingleSliceQAOA

VariationalCircuit (abstract)
├── QAOACircuit
├── pQAOACircuit
│   └── pQAOASingleParametersCircuit
└── SingleSliceQAOACircuit  (inherits QAOACircuit)

CostFunction (abstract)
└── QAOACostFunction

Optimizer (abstract)
├── COBYLA
├── SPSA
└── MonteCarlo

AbstractQuantumObjects (abstract)
└── CirqQuantumObjects

Each BaseQAOA subclass wires together one VariationalCircuit, one QAOACostFunction, and an Optimizer. Users interact only with the BaseQAOA layer; the circuit and cost-function objects are set up automatically on construction.

---

Reference

Michele Cattelan and Sheir Yarkoni,
Parallel Circuit Implementation of Variational Quantum Algorithms,
npj Quantum Information, 2025.
https://www.nature.com/articles/s41534-025-00982-6