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
    quantum_objects=None,     # AbstractQuantumObjects backend; CirqQuantumObjects if None
    layout=None,              # {label: physical qubit}; see "Variables, labels and layout"
    variable_order=None,      # labels of Q's rows (any hashable values); default 0..n-1
)

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,
    quantum_objects=None,
    layout=None,
    variable_order=None,
)

pQAOA

Parallel QAOA: each sub-problem slice runs its own QAOA circuit with independent $2p$ hyperparameters. Every variable must belong to exactly one slice; partition says which (by default, slices take consecutive rows of Q).

pQAOA(
    Q,                        # (n_full, n_full) full-model QUBO
    slice_Qs=None,            # per-slice QUBOs; None → the diagonal blocks of Q (needs partition)
    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
    cartesian=False,          # True → Cartesian-product gluing; False → zip gluing
    quantum_objects=None,
    partition=None,           # label lists: partition[k][i] = label of slice k's row i;
                              #   None → consecutive rows of Q, in slice order
    variable_order=None,      # labels of Q's rows; required unless the partition is sorted ints
    layout=None,              # {label: physical qubit}; qubits may be reused across slices
)

Named constructor:

pQAOA.from_ising(
    J, h, offset=0.0,
    slice_isings=None,        # (J_k, h_k[, offset_k]) tuples; None → blocks of J, h
    p=1,
    hyperparameters=None,
    cartesian=False,
    quantum_objects=None,
    partition=None,
    variable_order=None,
    layout=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).

Choosing the slices: pass partition, with variables named however you like (see Variables, labels and layout). The slice models default to the diagonal blocks of Q, so often Q and partition are all you need:

pQAOA(Q, partition=[[0, 1, 3], [4, 5]])                  # 5 variables named 0, 1, 3, 4, 5
pQAOA(Q, partition=[[0, 3, 5], [1, 2, 4]],               # Q's rows are variables 0..5,
      variable_order=[0, 1, 2, 3, 4, 5])                 #   so say so
pQAOA(Q, partition=[["a", "d"], ["b", "c"]],
      variable_order=["a", "b", "c", "d"])               # string labels

Pass slice_Qs (or slice_isings) explicitly when your slice models differ from the plain blocks; row $i$ of slice_Qs[k] is variable partition[k][i]. The defaults are the block of the model you passed: with from_ising, the blocks of J and h, which differ from the Ising form of the QUBO blocks when slices are coupled.

pQAOASingleParameters

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

pQAOASingleParameters(
    Q, slice_Qs=None,
    offset=0.0, slice_offsets=None,
    p=1,
    hyperparameters=None,     # (2*p,) array; random if None
    cartesian=False,
    quantum_objects=None,
    partition=None,
    variable_order=None,
    layout=None,
)

Named constructor:

pQAOASingleParameters.from_ising(
    J, h, offset=0.0,
    slice_isings=None,
    p=1,
    hyperparameters=None,
    cartesian=False,
    quantum_objects=None,
    partition=None,
    variable_order=None,
    layout=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,                  # labels, (n_slice, num_slices):
                              #   var_map[i, k] = the variable that slice var i is in copy k
    offset=0.0,
    slice_offset=0.0,
    p=1,
    hyperparameters=None,
    quantum_objects=None,
    layout=None,              # {label: physical qubit}; var_map[i, :] must share one qubit
    variable_order=None,      # labels of Q's rows; required unless the copies are sorted ints
)

Named constructor:

SingleSliceQAOA.from_ising(
    J, h, offset=0.0,
    slice_J, slice_h, slice_offset=0.0,
    var_map,
    p=1,
    hyperparameters=None,
    quantum_objects=None,
    layout=None,
    variable_order=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

Variables, labels and layout

Variables are named with labels of your choice — integers, strings, tuples, anything hashable. Two inputs tie the labels to the matrices:

• variable_order[i] is the label of row $i$ of Q (or of J, h) — and of column $i$ of the returned samples, which always come back in Q's row order.
• partition[k][i] is the label of row $i$ of slice $k$ (and of slice_Qs[k]). For SingleSliceQAOA, var_map[i, k] plays this role.

variable_order can be omitted when it is unambiguous: with no partition (labels are 0..n-1), or when the partition lists integer labels in increasing order — [[0, 1, 3], [4, 5]] means Q's rows are 0, 1, 3, 4, 5. Otherwise it is required and a ValueError says so; for a natural-order Q split as [[0, 3, 5], [1, 2, 4]], pass variable_order=[0, 1, 2, 3, 4, 5]. The partition must name every variable of variable_order exactly once.

alg = pQAOA(Q, partition=[["a", "d"], ["b", "c"]], variable_order=["a", "b", "c", "d"])
x = alg.circuit.sample(num_samples=100, bitstring=True)    # x[:, i] is variable "abcd"[i]

from bqs.utils import samples_to_dicts
samples_to_dicts(x, alg.circuit.variable_order)[0]          # {'a': 1, 'b': 0, 'c': 0, 'd': 1}

samples_to_dicts is a convenience for reading results by name; it builds one dict per sample, so keep the array for large batches.

layout says which physical qubit runs each variable: a dict {label: qubit}. Values are backend qubit objects such as cirq.GridQubit(3, 4), or integers passed to the backend's qubit(). It has no effect on the samples, only on which qubits the circuits use. Every variable must be mapped and the qubits must be distinct within one circuit; different slices may reuse the same qubits. Without a layout, integer labels run on quantum_objects.qubit(label) and other labels on quantum_objects.qubit(row position).

from cirq import GridQubit
good = [GridQubit(3, 4), GridQubit(3, 5)]
pQAOA(Q, partition=[["a", "d"], ["b", "c"]], variable_order=["a", "b", "c", "d"],
      layout={"a": good[0], "d": good[1],      # slice 0
              "b": good[0], "c": good[1]})     # slice 1 reuses the same qubits

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, bitstring=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], cartesian=True)))

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

# pQAOA: non-contiguous slices {0, 2} and {1, 3, 4, 5}; slice models default to the blocks of Q
print("pQAOA (partition) energy:",
      run(pQAOA(Q_full, partition=[[0, 2], [1, 3, 4, 5]], variable_order=range(N))))

# 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;
        column j is the j-th key of qubit_name_to_object (dict iteration order).
        Keys are not necessarily 0..n-1: pQAOA slices use their full-model
        indices (e.g. {0, 3, 5}), so iterate over the dict instead of range(n).
        """
        ...

    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, variable_order=None, quantum_objects=None):
        super().__init__(Q=Q, my_extra_param=my_extra_param, offset=offset, p=p,
                          hyperparameters=hyperparameters, variable_order=variable_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):
        self.hyperparameters = self._validate_hp_2d(self.hyperparameters, self._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

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, variable_order=None, quantum_objects=None):
        super().__init__(Q=Q, warm_start=np.asarray(warm_start, dtype=float),
                          offset=offset, p=p, hyperparameters=hyperparameters,
                          variable_order=variable_order, quantum_objects=quantum_objects)

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

    # _set_objective_function is inherited from BaseQAOA: it builds a QAOACostFunction from
    # circuit._J / _h / _ising_offset (samples come back in Q's row order, like J and h).

    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 and get_circuit. BaseQAOA already provides _set_objective_function, scoring samples with the circuit's J, h (samples are in Q's row order).
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

# Read samples by variable label
from bqs.utils import samples_to_dicts
samples_to_dicts(samples, alg.circuit.variable_order)   # [{label: value, ...}, ...]

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