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 implemented 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:

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)

Building new algorithm classes

Beyond swapping the backend or the optimizer, you can add an entirely new algorithm variant to the framework — e.g. your own flavor of QAOA — by subclassing BaseQAOA (or BaseVQA directly, if you're not doing QAOA at all). This is the most powerful extension point, but also the least "plug-and-play" one, so read this section fully before starting.

The three interfaces involved

Every algorithm in bqs is assembled from three pieces:

A new algorithm variant usually means writing a BaseQAOA subclass paired with a matching VariationalCircuit subclass (and, only if your scoring differs from plain Ising energy, a CostFunction subclass). optimize() itself is already implemented on BaseQAOA and rarely needs to be overridden.

How construction works

BaseVQA.__init__ takes the constructor arguments as keyword attributes and sets them on self before calling the hooks that need them:

class BaseVQA(ABC):
    def __init__(self, **attributes):
        for name, value in attributes.items():
            setattr(self, name, value)
        self.optimizer = None
        self._validate_hyperparameters()
        self._set_circuit()
        self._set_objective_function()

So a subclass never assigns anything to self before calling super().__init__() — it just forwards everything as keywords, and by the time _validate_hyperparameters/_set_circuit/_set_objective_function run, every attribute is already there:

class MyQAOAVariant(BaseQAOA):
    def __init__(self, Q, my_extra_param, offset=0.0, p=1,
                 hyperparameters=None, qubit_order=None, quantum_objects=None):
        super().__init__(Q=Q, my_extra_param=my_extra_param, offset=offset, p=p,
                          hyperparameters=hyperparameters, qubit_order=qubit_order,
                          quantum_objects=quantum_objects)

Adding a from_ising constructor

If your variant should also accept (J, h, offset) directly, don't branch on it inside __init__. Build the attribute dict for the Ising case and hand it to cls._from_attributes(...), which allocates the instance and runs the same BaseVQA.__init__ path — no cls.__new__ boilerplate to write yourself:

@classmethod
def from_ising(cls, J, h, offset=0.0, p=1, my_extra_param=None, **kwargs):
    return cls._from_attributes(
        Q=None,
        _J_in=np.asarray(J, dtype=float),   # stash Ising inputs under a private name
        _h_in=np.asarray(h, dtype=float),
        offset=offset, p=p, my_extra_param=my_extra_param, **kwargs,
    )

Your _set_circuit() then checks if self.Q is not None and falls back to self._J_in / self._h_in otherwise — that's the branch point, not the constructor.

Non-flat hyperparameters

BaseQAOA assumes a flat (2*p,) hyperparameter array by default, via the _validate_hyperparameters hook:

class BaseQAOA(BaseVQA):
    def _validate_hyperparameters(self):
        self.hyperparameters = self._validate_hp_1d(self.hyperparameters, self.p)

If your variant's hyperparameters don't fit that shape — e.g. pQAOA's per-slice (n_slices, 2*p) array — override _validate_hyperparameters to normalize/validate whatever shape you need; everything else about construction stays the same:

class pQAOA(BaseQAOA):
    def _validate_hyperparameters(self):
        n_slices = len(self.slice_Qs) if self.Q is not None else len(self._slice_isings)
        self.hyperparameters = self._validate_hp_2d(self.hyperparameters, n_slices, self.p)

Worked example

A small variant of QAOA that adds a caller-supplied bias to the local fields before sampling:

import numpy as np
from bqs.algorithms.qaoa import BaseQAOA
from bqs.utils import QAOACircuit, QAOACostFunction

class WarmStartQAOA(BaseQAOA):
    """QAOA that nudges the Ising local fields by a fixed 'warm_start' bias
    before the circuit is built — a simple way to bias sampling toward a
    known good solution."""

    def __init__(self, Q, warm_start, offset=0.0, p=1,
                 hyperparameters=None, qubit_order=None, quantum_objects=None):
        super().__init__(Q=Q, warm_start=np.asarray(warm_start, dtype=float),
                          offset=offset, p=p, hyperparameters=hyperparameters,
                          qubit_order=qubit_order, quantum_objects=quantum_objects)

    def _set_circuit(self):
        self.circuit = QAOACircuit(
            Q=self.Q, offset=self.offset, p=self.p,
            hyperparameters=self.hyperparameters, qubit_order=self.qubit_order,
            quantum_objects=self.quantum_objects,
        )
        self.circuit._h = self.circuit._h + self.warm_start   # apply the bias

    def _set_objective_function(self):
        self.objective_function = QAOACostFunction(
            self.circuit._J, self.circuit._h, self.circuit._ising_offset
        )

    def get_circuit(self):
        return self.circuit.qasm


alg = WarmStartQAOA(Q, warm_start=np.array([0.1, -0.1, 0.0]), p=2)

Reusing QAOACircuit and QAOACostFunction as-is (rather than writing new circuit/cost classes) is the cheapest way to prototype a variant — write your own VariationalCircuit subclass only once you actually need different gates or a different sampling scheme.

Writing a custom cost function

If your variant's scoring isn't plain Ising energy, subclass CostFunction:

from bqs.utils.vqa_utils import CostFunction

class MyCostFunction(CostFunction):
    def evaluate_samples(self, solutions):
        """
        solutions: list of dicts, {variable_name: measured_value}, or an
        equivalent batch of samples.
        Returns a single float — the mean energy/cost over the batch.
        """
        ...

Summary checklist for a new BaseQAOA subclass

1. Pass every subclass-specific value as a keyword to super().__init__(...) — there's no attribute-ordering rule to follow.
2. Decide whether BaseQAOA's 1-D hyperparameter validation fits your data. If not (e.g. 2-D per-slice hyperparameters), override _validate_hyperparameters.
3. Implement _set_circuit, _set_objective_function, and get_circuit.
4. Reuse existing VariationalCircuit/CostFunction classes where possible; only write new ones when gates or scoring genuinely differ.
5. If adding from_ising, build the Ising-case attribute dict and pass it to cls._from_attributes(...).

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 — see Building new algorithm classes if you're adding a branch to this tree yourself.

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