QUBO / Ising API

The quprep.qubo module provides QUBO and Ising problem formulation, QAOA circuit generation, and D-Wave export.

Classical reference solvers (solve_brute, solve_sa) are available in quprep.qubo.solver for benchmarking quantum results against classical baselines. They are not part of the public quprep.qubo namespace.

---

Core types

quprep.qubo.converter.QUBOResult(Q, offset=0.0, variable_map=None, n_original=None)

Result of a QUBO conversion.

Attributes:

Name	Type	Description
Q	(ndarray, shape(n, n))	Upper-triangular QUBO matrix. Q[i,i] are linear (bias) terms; Q[i,j] for i<j are quadratic couplings. Suitable for D-Wave samplers, QAOA, and OpenQAOA.
offset	float	Constant offset term (does not affect the optimal solution).
variable_map	dict	Mapping from variable names to matrix indices. Empty by default.
n_original	int	Number of original (non-slack) variables. Equals Q.shape[0] unless inequality constraints introduced slack variables.

Functions

evaluate(x)

Evaluate the QUBO objective for a given binary assignment.

Computes \(x^T Q x + ext{offset}\).

Parameters:

Name	Type	Description	Default
x	(array - like, shape(n))	Binary assignment vector with values in {0, 1}.	required

Returns:

Type	Description
float	Objective value including the constant offset.

Examples:

>>> import numpy as np
>>> from quprep.qubo.problems.maxcut import max_cut
>>> adj = np.array([[0,1,1],[1,0,1],[1,1,0]], dtype=float)
>>> q = max_cut(adj)
>>> q.evaluate(np.array([0, 1, 1]))  # cut between node 0 and {1,2}
-2.0

from_dict(d) classmethod

Deserialize from a dict produced by to_dict().

Parameters:

Name	Type	Description	Default
d	dict	Dict with keys Q, offset, variable_map, n_original.	required

Returns:

Type	Description
QUBOResult	Reconstructed QUBO problem.

to_dict()

Serialize to a plain Python dict (JSON-compatible).

Returns:

Type	Description
dict	Keys: Q (nested list), offset, variable_map, n_original. Pass directly to json.dumps() to serialize.

to_dwave()

Export QUBO as a D-Wave Ocean SDK-compatible linear/quadratic dict.

Returns a dict mapping (i, j) tuples to coefficients, where:

• i == j encodes linear (bias) terms (from the diagonal of Q)
• i < j encodes quadratic (coupling) terms (from upper triangle of Q)

Zero entries are omitted. The result can be passed directly to dimod.BinaryQuadraticModel.from_qubo() or dwave.system.DWaveSampler.

Returns:

Type	Description
dict	{(i, j): float} with i <= j.

Examples:

>>> import numpy as np
>>> from quprep.qubo.problems.maxcut import max_cut
>>> adj = np.array([[0,1],[1,0]], dtype=float)
>>> max_cut(adj).to_dwave()
{(0, 0): -1.0, (1, 1): -1.0, (0, 1): 2.0}

to_ising()

Convert QUBO to Ising (h, J) form.

Returns:

Type	Description
IsingResult	Ising model equivalent to this QUBO problem.

quprep.qubo.ising.IsingResult(h, J, offset=0.0) dataclass

Ising model representation.

Attributes:

Name	Type	Description
h	(ndarray, shape(n))	Linear (bias) coefficients.
J	(ndarray, shape(n, n))	Quadratic (coupling) coefficients, upper-triangular (J[i,j] for i < j).
offset	float	Constant energy offset.

Functions

to_qubo()

Convert Ising (h, J) back to QUBO form.

Uses the inverse transformation \(x_i = (s_i + 1) / 2\):

\[ Q_{ii} = 2h_i - 2\sum_{j \neq i} J^{\text{sym}}_{ij}, \quad Q_{ij} = 4J_{ij} \; (i<j), \quad \text{offset} = \text{offset}_0 + \textstyle\sum J_{\text{upper}} - \sum_i h_i \]

Returns:

Type	Description
QUBOResult	QUBO form of the Ising model.

---

Conversion

quprep.qubo.converter.to_qubo(cost_matrix, constraints=None, penalty=10.0)

Convert a cost matrix and optional linear equality constraints to QUBO.

The QUBO objective is: minimize \(x^T Q x\), \(x \in \{0, 1\}^n\)

The input cost_matrix can be any square real matrix. It is converted to upper-triangular QUBO form:

\[ Q_{ii} = M_{ii} \; \text{(linear/bias)}, \quad Q_{ij} = M_{ij} + M_{ji} \; (i < j) \; \text{(quadratic coupling)} \]

Constraint penalties are added on top via the Lagrangian approach.

Parameters:

Name	Type	Description	Default
cost_matrix	(ndarray, shape(n, n))	Quadratic cost matrix. Diagonal entries encode linear terms.	required
constraints	list of dicts	Equality and inequality constraints. Each dict must have: "A" : np.ndarray, shape (m, n) or (n,) "b" : np.ndarray or scalar "type" : "eq" (default) or "ineq" (Ax <= b via slack variables) "penalty" : float (optional, falls back to global penalty) Inequality constraints expand the variable count by adding binary slack variables. result.n_original records the original count.	None
penalty	float	Default Lagrange multiplier. Heuristic: 10x the largest |cost| entry.	10.0

Returns:

Type	Description
QUBOResult	result.n_original = n (original variables). If inequality constraints are present, result.Q.shape[0] > n.

Source code in quprep/qubo/converter.py

def to_qubo(
    cost_matrix: np.ndarray,
    constraints: list[dict] | None = None,
    penalty: float = 10.0,
) -> QUBOResult:
    r"""
    Convert a cost matrix and optional linear equality constraints to QUBO.

    The QUBO objective is: minimize $x^T Q x$, $x \in \{0, 1\}^n$

    The input `cost_matrix` can be any square real matrix. It is converted to
    upper-triangular QUBO form:

    $$
    Q_{ii} = M_{ii} \; \text{(linear/bias)}, \quad
    Q_{ij} = M_{ij} + M_{ji} \; (i < j) \; \text{(quadratic coupling)}
    $$

    Constraint penalties are added on top via the Lagrangian approach.

    Parameters
    ----------
    cost_matrix : np.ndarray, shape (n, n)
        Quadratic cost matrix. Diagonal entries encode linear terms.
    constraints : list of dicts, optional
        Equality and inequality constraints. Each dict must have:
            "A"       : np.ndarray, shape (m, n) or (n,)
            "b"       : np.ndarray or scalar
            "type"    : "eq" (default) or "ineq" (Ax <= b via slack variables)
            "penalty" : float (optional, falls back to global penalty)
        Inequality constraints expand the variable count by adding binary
        slack variables. ``result.n_original`` records the original count.
    penalty : float
        Default Lagrange multiplier. Heuristic: 10x the largest |cost| entry.

    Returns
    -------
    QUBOResult
        ``result.n_original`` = n (original variables).
        If inequality constraints are present, ``result.Q.shape[0] > n``.
    """
    M = np.asarray(cost_matrix, dtype=float)
    if M.ndim != 2 or M.shape[0] != M.shape[1]:
        raise ValueError(f"cost_matrix must be square 2-D, got shape {M.shape}.")
    n = M.shape[0]

    # Build upper-triangular QUBO from cost matrix
    Q = np.zeros((n, n))
    np.fill_diagonal(Q, np.diag(M))
    upper_mask = np.triu(np.ones((n, n), dtype=bool), k=1)
    Q[upper_mask] = M[upper_mask] + M.T[upper_mask]

    total_offset = 0.0
    total_slack = 0

    if constraints:
        from quprep.qubo.constraints import equality_penalty, inequality_penalty

        for c in constraints:
            A = np.asarray(c["A"], dtype=float)
            b = np.asarray(c["b"], dtype=float)
            lam = float(c.get("penalty", penalty))
            ctype = c.get("type", "eq")

            if ctype == "ineq":
                # Expand Q to accommodate slack variables
                Q_pen, pen_offset, n_slack = inequality_penalty(A, b, lam)
                # Pad current Q to match expanded size
                new_size = Q.shape[0] + n_slack
                Q_expanded = np.zeros((new_size, new_size))
                Q_expanded[: Q.shape[0], : Q.shape[0]] = Q
                Q = Q_expanded + Q_pen
                total_offset += pen_offset
                total_slack += n_slack
            else:
                # Equality constraint — pad A to current Q size if slack was added
                cur_n = Q.shape[0]
                if A.ndim == 1:
                    A = A.reshape(1, -1)
                if A.shape[1] < cur_n:
                    pad = np.zeros((A.shape[0], cur_n - A.shape[1]))
                    A = np.concatenate([A, pad], axis=1)
                Q_pen, pen_offset = equality_penalty(A, b, lam)
                Q += Q_pen
                total_offset += pen_offset

    return QUBOResult(Q=Q, offset=total_offset, n_original=n)

quprep.qubo.ising.qubo_to_ising(qubo)

Convert a QUBOResult to Ising (h, J) form.

Uses the substitution \(x_i = (s_i + 1) / 2\):

\[ J_{ij} = \frac{Q_{ij}}{4} \; (i<j), \quad h_i = \frac{Q_{ii}}{2} + \sum_{j \neq i} \frac{Q^{\text{sym}}_{ij}}{4}, \quad \text{offset} = \text{offset}_0 + \frac{\sum_i Q_{ii}}{2} + \frac{\sum_{i<j} Q_{ij}}{4} \]

Parameters:

Name	Type	Description	Default
qubo	QUBOResult	QUBO problem in upper-triangular form.	required

Returns:

Type	Description
IsingResult	Ising model with bias vector h and coupling matrix J.

Source code in quprep/qubo/ising.py

def qubo_to_ising(qubo) -> IsingResult:
    r"""
    Convert a QUBOResult to Ising (h, J) form.

    Uses the substitution $x_i = (s_i + 1) / 2$:

    $$
    J_{ij} = \frac{Q_{ij}}{4} \; (i<j), \quad
    h_i = \frac{Q_{ii}}{2} + \sum_{j \neq i} \frac{Q^{\text{sym}}_{ij}}{4}, \quad
    \text{offset} = \text{offset}_0 + \frac{\sum_i Q_{ii}}{2} + \frac{\sum_{i<j} Q_{ij}}{4}
    $$

    Parameters
    ----------
    qubo : QUBOResult
        QUBO problem in upper-triangular form.

    Returns
    -------
    IsingResult
        Ising model with bias vector h and coupling matrix J.
    """
    Q = qubo.Q

    # Symmetric off-diagonal view — diagonal must be zero so row sums are purely off-diagonal
    Q_sym = Q + Q.T - 2.0 * np.diag(np.diag(Q))

    # Coupling matrix (upper-triangular)
    J = np.triu(Q, k=1) / 4.0

    # Bias vector
    off_diag_row_sums = Q_sym.sum(axis=1)  # Q_sym[i,i] = 0 by construction
    h = np.diag(Q) / 2.0 + off_diag_row_sums / 4.0

    # Constant offset
    offset = (
        qubo.offset
        + np.sum(np.diag(Q)) / 2.0
        + np.sum(np.triu(Q, k=1)) / 4.0
    )

    return IsingResult(h=h, J=J, offset=offset)

quprep.qubo.ising.ising_to_qubo(ising)

Convert an Ising model back to QUBO form.

Applies the inverse substitution \(x_i = (s_i + 1) / 2\):

\[ Q_{ii} = 2h_i - 2\sum_{j \neq i} J^{\text{sym}}_{ij}, \quad Q_{ij} = 4J_{ij} \; (i<j), \quad \text{offset} = \text{offset}_0 + \sum_{\text{upper}} J - \sum_i h_i \]

Parameters:

Name	Type	Description	Default
ising	IsingResult	Ising model in (h, J) form.	required

Returns:

Type	Description
QUBOResult	QUBO form of the Ising model.

Examples:

>>> import numpy as np
>>> from quprep.qubo.problems.maxcut import max_cut
>>> from quprep.qubo.ising import ising_to_qubo, qubo_to_ising
>>> adj = np.array([[0,1,1],[1,0,1],[1,1,0]], dtype=float)
>>> q = max_cut(adj)
>>> q2 = ising_to_qubo(qubo_to_ising(q))
>>> np.allclose(q.Q, q2.Q)
True

Source code in quprep/qubo/ising.py

def ising_to_qubo(ising: IsingResult) -> QUBOResult:
    r"""
    Convert an Ising model back to QUBO form.

    Applies the inverse substitution $x_i = (s_i + 1) / 2$:

    $$
    Q_{ii} = 2h_i - 2\sum_{j \neq i} J^{\text{sym}}_{ij}, \quad
    Q_{ij} = 4J_{ij} \; (i<j), \quad
    \text{offset} = \text{offset}_0 + \sum_{\text{upper}} J - \sum_i h_i
    $$

    Parameters
    ----------
    ising : IsingResult
        Ising model in (h, J) form.

    Returns
    -------
    QUBOResult
        QUBO form of the Ising model.

    Examples
    --------
    >>> import numpy as np
    >>> from quprep.qubo.problems.maxcut import max_cut
    >>> from quprep.qubo.ising import ising_to_qubo, qubo_to_ising
    >>> adj = np.array([[0,1,1],[1,0,1],[1,1,0]], dtype=float)
    >>> q = max_cut(adj)
    >>> q2 = ising_to_qubo(qubo_to_ising(q))
    >>> np.allclose(q.Q, q2.Q)
    True
    """
    return ising.to_qubo()

---

Problem library

quprep.qubo.problems.max_cut(adjacency)

Build the QUBO for the Max-Cut problem.

Parameters:

Name	Type	Description	Default
adjacency	(ndarray, shape(n, n))	Weighted adjacency matrix of the graph. Symmetric; self-loops (diagonal) are ignored. For unweighted graphs, use a 0/1 matrix.	required

Returns:

Type	Description
QUBOResult	Variable i=0..n-1 is 1 if node i is in partition S, 0 otherwise. Minimizing the QUBO objective maximises the cut weight.

Source code in quprep/qubo/problems/maxcut.py

def max_cut(adjacency: np.ndarray) -> QUBOResult:
    """
    Build the QUBO for the Max-Cut problem.

    Parameters
    ----------
    adjacency : np.ndarray, shape (n, n)
        Weighted adjacency matrix of the graph.
        Symmetric; self-loops (diagonal) are ignored.
        For unweighted graphs, use a 0/1 matrix.

    Returns
    -------
    QUBOResult
        Variable i=0..n-1 is 1 if node i is in partition S, 0 otherwise.
        Minimizing the QUBO objective maximises the cut weight.
    """
    adj = np.asarray(adjacency, dtype=float)
    n = adj.shape[0]
    # Symmetrize in case input is directed
    W = (adj + adj.T) / 2.0
    np.fill_diagonal(W, 0.0)

    Q = np.zeros((n, n))
    for i in range(n):
        for j in range(i + 1, n):
            w = W[i, j]
            if w != 0.0:
                Q[i, j] += 2.0 * w
                Q[i, i] -= w
                Q[j, j] -= w

    var_map = {f"x{i}": i for i in range(n)}
    return QUBOResult(Q=Q, offset=0.0, variable_map=var_map)

quprep.qubo.problems.knapsack

0/1 Knapsack QUBO formulation.

Knapsack: given n items with values v_i and weights w_i, select a subset to maximize total value without exceeding capacity W.

QUBO formulation (minimization):

\[\min -\sum_i v_i x_i + \lambda\!\left(\sum_i w_i x_i - W\right)^2\]

The capacity constraint is enforced via a quadratic penalty \(\lambda\). The penalty coefficient should be at least \(\max(v_i)\) to ensure feasibility.

Note

This is a penalty-based (soft) formulation. Infeasible solutions have higher energy, but the exact cut-off depends on penalty strength. For hard enforcement use slack binary variables (adds ceil(log2(W+1)) ancilla qubits).

References

Lucas, A. (2014). Ising formulations of many NP problems. Frontiers in Physics, 2, 5. doi:10.3389/fphy.2014.00005

Classes

Functions

knapsack(weights, values, capacity, penalty=None)

Build the QUBO for the 0/1 Knapsack problem.

Parameters:

Name	Type	Description	Default
weights	(ndarray, shape(n))	Item weights.	required
values	(ndarray, shape(n))	Item values (non-negative).	required
capacity	float	Maximum total weight W.	required
penalty	float	Lagrange multiplier for the capacity constraint. Defaults to max(values) + 1, which is tight enough to enforce feasibility for integer-valued instances.	None

Returns:

Type	Description
QUBOResult	Variable x_i = 1 means item i is selected.

Raises:

Type	Description
ValueError	If weights and values have different lengths.

Source code in quprep/qubo/problems/knapsack.py

def knapsack(
    weights: np.ndarray,
    values: np.ndarray,
    capacity: float,
    penalty: float | None = None,
) -> QUBOResult:
    """
    Build the QUBO for the 0/1 Knapsack problem.

    Parameters
    ----------
    weights : np.ndarray, shape (n,)
        Item weights.
    values : np.ndarray, shape (n,)
        Item values (non-negative).
    capacity : float
        Maximum total weight W.
    penalty : float, optional
        Lagrange multiplier for the capacity constraint.
        Defaults to max(values) + 1, which is tight enough to enforce
        feasibility for integer-valued instances.

    Returns
    -------
    QUBOResult
        Variable x_i = 1 means item i is selected.

    Raises
    ------
    ValueError
        If ``weights`` and ``values`` have different lengths.
    """
    w = np.asarray(weights, dtype=float)
    v = np.asarray(values, dtype=float)
    n = len(w)
    if len(v) != n:
        raise ValueError("weights and values must have the same length.")

    if penalty is None:
        penalty = float(np.max(v)) + 1.0

    Q = np.zeros((n, n))
    offset = 0.0

    # Objective: minimize -v_i * x_i
    np.fill_diagonal(Q, Q.diagonal() - v)

    # Capacity penalty: penalty * (sum w_i x_i - W)^2
    # Diagonal: penalty * (w_i^2 - 2*W*w_i)
    np.fill_diagonal(Q, Q.diagonal() + penalty * (w ** 2 - 2.0 * capacity * w))
    # Off-diagonal: penalty * 2 * w_i * w_j
    outer = np.outer(w, w)
    Q += penalty * 2.0 * np.triu(outer, k=1)
    offset += penalty * capacity ** 2

    var_map = {f"x{i}": i for i in range(n)}
    return QUBOResult(Q=Q, offset=offset, variable_map=var_map)

quprep.qubo.problems.tsp

Travelling Salesman Problem (TSP) QUBO formulation.

TSP: given n cities with distance matrix D, find the shortest Hamiltonian cycle visiting every city exactly once.

QUBO formulation uses \(n^2\) binary variables \(x_{i,t}\) where \(x_{i,t}=1\) means city \(i\) is visited at time step \(t\).

Objective (minimize total distance):

\[\min \sum_{i,j,t} D_{ij}\, x_{i,t}\, x_{j,(t+1)\bmod n}\]

Constraints (enforced via quadratic penalties):

• \(C_1\): each city visited once: \(\sum_t x_{i,t} = 1\) for all \(i\)
• \(C_2\): one city per time slot: \(\sum_i x_{i,t} = 1\) for all \(t\)

Variable index: \(v(i, t) = i \cdot n + t\).

References

Lucas, A. (2014). Ising formulations of many NP problems. Frontiers in Physics, 2, 5. doi:10.3389/fphy.2014.00005

Classes

Functions

tsp(distance_matrix, penalty=None)

Build the QUBO for the Travelling Salesman Problem.

Parameters:

Name	Type	Description	Default
distance_matrix	(ndarray, shape(n, n))	Pairwise distances between cities. Asymmetric matrices are supported. Self-distances (diagonal) are ignored.	required
penalty	float	Lagrange multiplier for both city and time-slot constraints. Defaults to max(D) * n, which is typically large enough to enforce feasibility while keeping constraint violations expensive.	None

Returns:

Type	Description
QUBOResult	\(n^2\) binary variables. variable_map["x_i_t"] gives the index of \(x_{i,t}\). A feasible solution has exactly n variables equal to 1, one per city and one per time step.

Raises:

Type	Description
ValueError	If distance_matrix is not a square 2-D array.

Source code in quprep/qubo/problems/tsp.py

def tsp(distance_matrix: np.ndarray, penalty: float | None = None) -> QUBOResult:
    """
    Build the QUBO for the Travelling Salesman Problem.

    Parameters
    ----------
    distance_matrix : np.ndarray, shape (n, n)
        Pairwise distances between cities. Asymmetric matrices are supported.
        Self-distances (diagonal) are ignored.
    penalty : float, optional
        Lagrange multiplier for both city and time-slot constraints.
        Defaults to max(D) * n, which is typically large enough to enforce
        feasibility while keeping constraint violations expensive.

    Returns
    -------
    QUBOResult
        $n^2$ binary variables. ``variable_map["x_i_t"]`` gives the index of $x_{i,t}$.
        A feasible solution has exactly n variables equal to 1, one per city
        and one per time step.

    Raises
    ------
    ValueError
        If ``distance_matrix`` is not a square 2-D array.
    """
    D = np.asarray(distance_matrix, dtype=float)
    n = D.shape[0]
    if D.shape != (n, n):
        raise ValueError(f"distance_matrix must be square, got {D.shape}.")

    if penalty is None:
        penalty = float(np.max(D)) * n + 1.0

    N = n * n  # total variables
    Q = np.zeros((N, N))
    offset = 0.0

    def idx(city: int, step: int) -> int:
        return city * n + step

    # ---- Objective: minimize total tour distance ----
    for i in range(n):
        for j in range(n):
            if i == j:
                continue
            d = D[i, j]
            for t in range(n):
                t_next = (t + 1) % n
                a = idx(i, t)
                b = idx(j, t_next)
                if a < b:
                    Q[a, b] += d
                elif a > b:
                    Q[b, a] += d
                else:
                    Q[a, a] += d  # shouldn't happen when i != j

    # ---- Constraint 1: each city visited exactly once ----
    # penalty * (sum_t x[i,t] - 1)^2  for each city i
    for i in range(n):
        for t1 in range(n):
            a = idx(i, t1)
            Q[a, a] += penalty * (1.0 - 2.0)  # diagonal: penalty*(1 - 2*1)
            for t2 in range(t1 + 1, n):
                b = idx(i, t2)
                Q[a, b] += penalty * 2.0
    offset += penalty * n  # n cities x penalty * 1^2

    # ---- Constraint 2: each time slot has exactly one city ----
    # penalty * (sum_i x[i,t] - 1)^2  for each time step t
    for t in range(n):
        for i1 in range(n):
            a = idx(i1, t)
            Q[a, a] += penalty * (1.0 - 2.0)
            for i2 in range(i1 + 1, n):
                b = idx(i2, t)
                Q[a, b] += penalty * 2.0
    offset += penalty * n  # n time slots x penalty * 1^2

    var_map = {f"x_{i}_{t}": idx(i, t) for i in range(n) for t in range(n)}
    return QUBOResult(Q=Q, offset=offset, variable_map=var_map)

quprep.qubo.problems.portfolio

Markowitz Portfolio Optimization QUBO formulation.

Portfolio optimization: given n assets with expected returns mu_i and covariance matrix Sigma, select exactly K assets to maximize risk-adjusted return.

QUBO formulation (minimization):

\[\min -\sum_i \mu_i x_i + \lambda_r\, x^T \Sigma x + \lambda_b\!\left(\sum_i x_i - K\right)^2\]

References

Mugel et al. (2022). Dynamic portfolio optimization with real datasets using quantum processors and quantum-inspired tensor networks. Physical Review Research, 4(1), 013006. doi:10.1103/PhysRevResearch.4.013006

Classes

Functions

portfolio(returns, covariance, budget, risk_penalty=1.0, budget_penalty=None)

Build the QUBO for Markowitz portfolio optimization.

Parameters:

Name	Type	Description	Default
returns	(ndarray, shape(n))	Expected return for each asset.	required
covariance	(ndarray, shape(n, n))	Return covariance matrix (positive semi-definite).	required
budget	int	Number of assets to select (K).	required
risk_penalty	float	Lagrange multiplier for the risk (variance) term. Higher values favour lower-risk portfolios. Default is 1.0.	1.0
budget_penalty	float	Lagrange multiplier enforcing the budget constraint sum(x) = K. Defaults to max(|returns|) * n, which is generally strong enough.	None

Returns:

Type	Description
QUBOResult	Variable x_i = 1 means asset i is selected.

Raises:

Type	Description
ValueError	If covariance shape does not match the number of assets in returns.

Source code in quprep/qubo/problems/portfolio.py

def portfolio(
    returns: np.ndarray,
    covariance: np.ndarray,
    budget: int,
    risk_penalty: float = 1.0,
    budget_penalty: float | None = None,
) -> QUBOResult:
    """
    Build the QUBO for Markowitz portfolio optimization.

    Parameters
    ----------
    returns : np.ndarray, shape (n,)
        Expected return for each asset.
    covariance : np.ndarray, shape (n, n)
        Return covariance matrix (positive semi-definite).
    budget : int
        Number of assets to select (K).
    risk_penalty : float
        Lagrange multiplier for the risk (variance) term. Higher values
        favour lower-risk portfolios. Default is 1.0.
    budget_penalty : float, optional
        Lagrange multiplier enforcing the budget constraint sum(x) = K.
        Defaults to max(|returns|) * n, which is generally strong enough.

    Returns
    -------
    QUBOResult
        Variable x_i = 1 means asset i is selected.

    Raises
    ------
    ValueError
        If ``covariance`` shape does not match the number of assets in ``returns``.
    """
    mu = np.asarray(returns, dtype=float)
    Sigma = np.asarray(covariance, dtype=float)
    n = len(mu)
    if Sigma.shape != (n, n):
        raise ValueError(f"covariance must be ({n}, {n}), got {Sigma.shape}.")

    if budget_penalty is None:
        budget_penalty = float(np.max(np.abs(mu)) * n) + 1.0

    Q = np.zeros((n, n))
    offset = 0.0

    # Objective term: -mu_i * x_i
    np.fill_diagonal(Q, Q.diagonal() - mu)

    # Risk term: risk_penalty * x^T Sigma x
    # Diagonal: risk_penalty * Sigma[i,i]
    np.fill_diagonal(Q, Q.diagonal() + risk_penalty * np.diag(Sigma))
    # Off-diagonal: risk_penalty * 2 * Sigma[i,j] (symmetrize)
    Sigma_sym = (Sigma + Sigma.T) / 2.0
    Q += risk_penalty * 2.0 * np.triu(Sigma_sym, k=1)

    # Budget constraint: budget_penalty * (sum x_i - K)^2
    # Diagonal: budget_penalty * (1 - 2K)
    np.fill_diagonal(Q, Q.diagonal() + budget_penalty * (1.0 - 2.0 * budget))
    # Off-diagonal: budget_penalty * 2
    ones = np.ones((n, n))
    Q += budget_penalty * 2.0 * np.triu(ones, k=1)
    offset += budget_penalty * float(budget) ** 2

    var_map = {f"x{i}": i for i in range(n)}
    return QUBOResult(Q=Q, offset=offset, variable_map=var_map)

quprep.qubo.problems.graph_color

Graph Coloring QUBO formulation.

Graph coloring: assign one of \(K\) colors to each node in graph \(G=(V, E)\) such that no two adjacent nodes share the same color.

QUBO uses \(n \cdot K\) binary variables \(x_{i,c}\) where \(x_{i,c}=1\) means node \(i\) is assigned color \(c\).

Constraints (both enforced as penalty terms):

• \(C_1\) (one color per node): \(\sum_c x_{i,c} = 1\) for all \(i\)
• \(C_2\) (valid coloring): \(x_{i,c} \cdot x_{j,c} = 0\) for all edges \((i,j)\), all \(c\)

Variable index: \(v(i, c) = i \cdot K + c\)

References

Lucas, A. (2014). Ising formulations of many NP problems. Frontiers in Physics, 2, 5. doi:10.3389/fphy.2014.00005

Classes

Functions

graph_color(adjacency, n_colors, penalty=10.0)

Build the QUBO for the Graph Coloring problem.

Parameters:

Name	Type	Description	Default
adjacency	(ndarray, shape(n, n))	Adjacency matrix of the graph (weighted or 0/1). Symmetric; diagonal is ignored.	required
n_colors	int	Number of colors K available for assignment.	required
penalty	float	Lagrange multiplier for both constraints. Should be large enough that any constraint violation costs more than the best feasible objective. Typically >= n * max_edge_weight.	10.0

Returns:

Type	Description
QUBOResult	n * n_colors binary variables. variable_map["x_{i}{c}"] gives the index of variable x. A feasible solution has exactly n variables set to 1, one per node.

Notes

A valid K-coloring may not exist for all graphs and K values. If the minimum energy solution has non-zero cost, the graph is not K-colorable.

Source code in quprep/qubo/problems/graph_color.py

def graph_color(
    adjacency: np.ndarray,
    n_colors: int,
    penalty: float = 10.0,
) -> QUBOResult:
    """
    Build the QUBO for the Graph Coloring problem.

    Parameters
    ----------
    adjacency : np.ndarray, shape (n, n)
        Adjacency matrix of the graph (weighted or 0/1).
        Symmetric; diagonal is ignored.
    n_colors : int
        Number of colors K available for assignment.
    penalty : float
        Lagrange multiplier for both constraints.
        Should be large enough that any constraint violation costs more than
        the best feasible objective. Typically >= n * max_edge_weight.

    Returns
    -------
    QUBOResult
        n * n_colors binary variables. variable_map["x_{i}_{c}"] gives the
        index of variable x_{i,c}. A feasible solution has exactly n variables
        set to 1, one per node.

    Notes
    -----
    A valid K-coloring may not exist for all graphs and K values. If the
    minimum energy solution has non-zero cost, the graph is not K-colorable.
    """
    adj = np.asarray(adjacency, dtype=float)
    n = adj.shape[0]
    W = (adj + adj.T) / 2.0
    np.fill_diagonal(W, 0.0)

    N = n * n_colors
    Q = np.zeros((N, N))
    offset = 0.0

    def idx(node: int, color: int) -> int:
        return node * n_colors + color

    # C1: each node has exactly one color
    # penalty * (sum_c x[i,c] - 1)^2 for each node i
    for i in range(n):
        for c1 in range(n_colors):
            a = idx(i, c1)
            Q[a, a] += penalty * (1.0 - 2.0)   # -penalty per variable
            for c2 in range(c1 + 1, n_colors):
                b = idx(i, c2)
                Q[a, b] += penalty * 2.0
    offset += penalty * n  # n constraints, each contributes penalty * 1^2

    # C2: no two adjacent nodes share a color
    # penalty * x[i,c] * x[j,c] for each edge (i,j) and each color c
    for i in range(n):
        for j in range(i + 1, n):
            w = W[i, j]
            if w == 0.0:
                continue
            for c in range(n_colors):
                a = idx(i, c)
                b = idx(j, c)
                if a > b:
                    a, b = b, a
                Q[a, b] += penalty * w

    var_map = {f"x_{i}_{c}": idx(i, c) for i in range(n) for c in range(n_colors)}
    return QUBOResult(Q=Q, offset=offset, variable_map=var_map)

quprep.qubo.problems.scheduling

Job Scheduling QUBO formulation.

Scheduling: assign n jobs to m machines to minimize total squared load (a proxy for balanced makespan minimization).

QUBO uses \(n \cdot m\) binary variables \(x_{i,k}\) where \(x_{i,k}=1\) means job \(i\) is assigned to machine \(k\).

Objective (minimize load imbalance):

\[\min \sum_k \left(\sum_i p_i\, x_{i,k}\right)^2\]

where \(p_i\) is the processing time of job \(i\).

Constraint: each job assigned to exactly one machine: \(\sum_k x_{i,k} = 1\) for all \(i\).

Variable index: \(v(i, k) = i \cdot m + k\).

Notes

This is a load-balancing formulation. It does not model job dependencies, deadlines, or preemption. For those, see the Travelling Salesman or QUBO references for richer scheduling models.

References

Lucas, A. (2014). Ising formulations of many NP problems. Frontiers in Physics, 2, 5. doi:10.3389/fphy.2014.00005

Classes

Functions

scheduling(processing_times, n_machines, penalty=None)

Build the QUBO for the job scheduling (load balancing) problem.

Parameters:

Name	Type	Description	Default
processing_times	(ndarray, shape(n))	Processing time of each job (non-negative).	required
n_machines	int	Number of machines m available.	required
penalty	float	Lagrange multiplier for the assignment constraint (each job assigned to exactly one machine). Defaults to sum(processing_times)^2 + 1, which is always large enough to enforce feasibility.	None

Returns:

Type	Description
QUBOResult	n * n_machines binary variables. variable_map["x_{i}{k}"] gives the index of x. A feasible solution has exactly n variables set to 1, one per job.

Source code in quprep/qubo/problems/scheduling.py

def scheduling(
    processing_times: np.ndarray,
    n_machines: int,
    penalty: float | None = None,
) -> QUBOResult:
    """
    Build the QUBO for the job scheduling (load balancing) problem.

    Parameters
    ----------
    processing_times : np.ndarray, shape (n,)
        Processing time of each job (non-negative).
    n_machines : int
        Number of machines m available.
    penalty : float, optional
        Lagrange multiplier for the assignment constraint (each job assigned
        to exactly one machine). Defaults to sum(processing_times)^2 + 1,
        which is always large enough to enforce feasibility.

    Returns
    -------
    QUBOResult
        n * n_machines binary variables. variable_map["x_{i}_{k}"] gives the
        index of x_{i,k}. A feasible solution has exactly n variables set to 1,
        one per job.
    """
    p = np.asarray(processing_times, dtype=float)
    n = len(p)
    m = n_machines

    if penalty is None:
        penalty = float(np.sum(p) ** 2) + 1.0

    N = n * m
    Q = np.zeros((N, N))
    offset = 0.0

    def idx(job: int, machine: int) -> int:
        return job * m + machine

    # Objective: minimize sum_k (sum_i p_i x[i,k])^2
    # Expanding: sum_k sum_{i,j} p_i p_j x[i,k] x[j,k]
    # For i == j (same job, same machine): p_i^2 x[i,k] (diagonal)
    # For i != j (diff jobs, same machine): 2*p_i*p_j x[i,k]*x[j,k] (off-diagonal)
    for k in range(m):
        for i in range(n):
            Q[idx(i, k), idx(i, k)] += p[i] ** 2
        for i in range(n):
            for j in range(i + 1, n):
                a = idx(i, k)
                b = idx(j, k)
                Q[a, b] += 2.0 * p[i] * p[j]

    # Constraint: each job assigned to exactly one machine
    # penalty * (sum_k x[i,k] - 1)^2 for each job i
    for i in range(n):
        for k1 in range(m):
            a = idx(i, k1)
            Q[a, a] += penalty * (1.0 - 2.0)   # -penalty per variable
            for k2 in range(k1 + 1, m):
                b = idx(i, k2)
                Q[a, b] += penalty * 2.0
    offset += penalty * n   # n constraints, each contributes penalty * 1^2

    var_map = {f"x_{i}_{k}": idx(i, k) for i in range(n) for k in range(m)}
    return QUBOResult(Q=Q, offset=offset, variable_map=var_map)

quprep.qubo.problems.number_partition

Number Partitioning QUBO formulation.

Number partitioning: given n positive numbers v_i, split them into two subsets A and B such that the difference |sum(A) - sum(B)| is minimised (ideally zero for a perfect partition).

Let \(x_i = 1\) if value \(i\) is in subset \(A\), \(0\) if in subset \(B\). Define \(s_i = 2x_i - 1 \in \{-1, +1\}\). A perfect partition requires:

\[\sum_i v_i s_i = 0\]

The objective to minimize is the squared difference:

\[\left(\sum_i v_i s_i\right)^2 = \left(\sum_i v_i (2x_i - 1)\right)^2\]

Expanding (using \(x_i^2 = x_i\) for binary \(x\)), with \(S = \sum_i v_i\):

\[Q_{ii} = v_i(v_i - S), \quad Q_{ij} = 2v_iv_j \; (i<j), \quad \text{offset} = S^2\]

Note: the penalty argument scales the entire objective — useful when combining with other QUBO terms via add_qubo().

References

Lucas, A. (2014). Ising formulations of many NP problems. Frontiers in Physics, 2, 5. doi:10.3389/fphy.2014.00005

Classes

Functions

number_partition(values, penalty=1.0)

Build the QUBO for the Number Partitioning problem.

Parameters:

Name	Type	Description	Default
values	(ndarray, shape(n))	Positive numbers to partition.	required
penalty	float	Global scale factor. Set > 1 when combining with other QUBO objectives to ensure the partition constraint dominates. Default is 1.0.	1.0

Returns:

Type	Description
QUBOResult	Variable x_i = 1 means value i goes into subset A. Minimum energy = 0 indicates a perfect partition exists. Minimum energy > 0 means no perfect partition; the solution with lowest energy is the most balanced achievable split.

Examples:

>>> import numpy as np
>>> from quprep.qubo.problems.number_partition import number_partition
>>> from quprep.qubo.solver import solve_brute
>>> v = np.array([3.0, 1.0, 1.0, 2.0, 2.0, 1.0])  # sum=10, perfect split at 5
>>> sol = solve_brute(number_partition(v))
>>> sol.energy   # should be 0.0 for a perfect partition
0.0

Source code in quprep/qubo/problems/number_partition.py

def number_partition(values: np.ndarray, penalty: float = 1.0) -> QUBOResult:
    """
    Build the QUBO for the Number Partitioning problem.

    Parameters
    ----------
    values : np.ndarray, shape (n,)
        Positive numbers to partition.
    penalty : float
        Global scale factor. Set > 1 when combining with other QUBO objectives
        to ensure the partition constraint dominates. Default is 1.0.

    Returns
    -------
    QUBOResult
        Variable x_i = 1 means value i goes into subset A.
        Minimum energy = 0 indicates a perfect partition exists.
        Minimum energy > 0 means no perfect partition; the solution with
        lowest energy is the most balanced achievable split.

    Examples
    --------
    >>> import numpy as np
    >>> from quprep.qubo.problems.number_partition import number_partition
    >>> from quprep.qubo.solver import solve_brute
    >>> v = np.array([3.0, 1.0, 1.0, 2.0, 2.0, 1.0])  # sum=10, perfect split at 5
    >>> sol = solve_brute(number_partition(v))
    >>> sol.energy   # should be 0.0 for a perfect partition
    0.0
    """
    v = np.asarray(values, dtype=float)
    n = len(v)
    S = float(np.sum(v))

    Q = np.zeros((n, n))
    offset = 0.0

    # (sum_i v_i s_i)^2 where s_i = 2x_i - 1
    # = (2 sum_i v_i x_i - S)^2
    # = 4 (sum_i v_i x_i)^2 - 4S sum_i v_i x_i + S^2
    #
    # Expanding (sum_i v_i x_i)^2:
    #   = sum_i v_i^2 x_i^2 + 2 sum_{i<j} v_i v_j x_i x_j
    #   = sum_i v_i^2 x_i   + 2 sum_{i<j} v_i v_j x_i x_j   (binary)
    #
    # Diagonal: 4 v_i^2 - 4S v_i = 4 v_i (v_i - S)
    np.fill_diagonal(Q, Q.diagonal() + penalty * 4.0 * v * (v - S))

    # Off-diagonal: 4 * 2 * v_i v_j / 2 ... wait let me redo:
    # 4*(sum v_i x_i)^2 = 4*sum_i v_i^2 x_i + 4*2*sum_{i<j} v_i v_j x_i x_j
    #                   = 4 v_i^2 x_i diagonal + 8 v_i v_j x_i x_j off-diag
    # Combined with -4S sum v_i x_i (diagonal only):
    # Q[i,i] = 4v_i^2 - 4S v_i = 4v_i(v_i - S)
    # Q[i,j] = 8 v_i v_j   for i < j
    # offset  = S^2
    #
    # (The fill_diagonal above already set Q[i,i] = penalty * 4 v_i(v_i-S))
    outer = np.outer(v, v)
    Q += penalty * 8.0 * np.triu(outer, k=1)
    offset += penalty * S ** 2

    var_map = {f"x{i}": i for i in range(n)}
    return QUBOResult(Q=Q, offset=offset, variable_map=var_map)

---

Classical reference solvers

Import from quprep.qubo.solver — not part of the quprep.qubo public namespace.

from quprep.qubo.solver import solve_brute, solve_sa, SolveResult

quprep.qubo.solver.SolveResult(x, energy, n_evaluated) dataclass

Result of a QUBO solve.

Attributes:

Name	Type	Description
x	(ndarray, shape(n))	Optimal binary assignment vector (values in {0, 1}).
energy	float	Optimal objective value including the constant offset.
n_evaluated	int	Number of binary strings evaluated.

quprep.qubo.solver.solve_brute(qubo, max_n=20)

Find the exact minimum of a QUBO by exhaustive enumeration.

Evaluates all \(2^n\) binary strings and returns the one with the lowest objective value \(x^T Q x + ext{offset}\).

Parameters:

Name	Type	Description	Default
qubo	QUBOResult	The QUBO problem to solve.	required
max_n	int	Safety limit on problem size. Raises ValueError if n > max_n. Default is 20 (2^20 = ~1M evaluations, runs in < 1s).	20

Returns:

Type	Description
SolveResult	Best solution found, its energy, and the number of states evaluated. SolveResult.x uses LSB-first ordering: x[0] corresponds to bit 0 (the least-significant bit of each enumerated integer). Ensure your QUBO matrix Q is built with the same variable ordering.

Raises:

Type	Description
ValueError	If n > max_n.

Examples:

>>> from quprep.qubo.problems.maxcut import max_cut
>>> from quprep.qubo.solver import solve_brute
>>> import numpy as np
>>> adj = np.array([[0,1,1],[1,0,1],[1,1,0]], dtype=float)
>>> result = solve_brute(max_cut(adj))
>>> result.energy   # max cut of triangle = -1 (one node vs two)
-1.0

Source code in quprep/qubo/solver.py

def solve_brute(qubo, max_n: int = 20) -> SolveResult:
    """
    Find the exact minimum of a QUBO by exhaustive enumeration.

    Evaluates all $2^n$ binary strings and returns the one with the lowest
    objective value $x^T Q x + \text{offset}$.

    Parameters
    ----------
    qubo : QUBOResult
        The QUBO problem to solve.
    max_n : int
        Safety limit on problem size. Raises ValueError if n > max_n.
        Default is 20 (2^20 = ~1M evaluations, runs in < 1s).

    Returns
    -------
    SolveResult
        Best solution found, its energy, and the number of states evaluated.
        ``SolveResult.x`` uses LSB-first ordering: ``x[0]`` corresponds to
        bit 0 (the least-significant bit of each enumerated integer). Ensure
        your QUBO matrix Q is built with the same variable ordering.

    Raises
    ------
    ValueError
        If n > max_n.

    Examples
    --------
    >>> from quprep.qubo.problems.maxcut import max_cut
    >>> from quprep.qubo.solver import solve_brute
    >>> import numpy as np
    >>> adj = np.array([[0,1,1],[1,0,1],[1,1,0]], dtype=float)
    >>> result = solve_brute(max_cut(adj))
    >>> result.energy   # max cut of triangle = -1 (one node vs two)
    -1.0
    """
    Q = qubo.Q
    n = Q.shape[0]
    if n > max_n:
        raise ValueError(
            f"n={n} exceeds max_n={max_n}. "
            "Pass max_n=<larger value> or use a heuristic solver."
        )

    best_x = np.zeros(n)
    best_energy = float("inf")

    for bits in range(1 << n):
        x = np.array([(bits >> i) & 1 for i in range(n)], dtype=float)
        energy = float(x @ Q @ x) + qubo.offset
        if energy < best_energy:
            best_energy = energy
            best_x = x.copy()

    return SolveResult(x=best_x, energy=best_energy, n_evaluated=1 << n)

quprep.qubo.solver.solve_sa(qubo, n_steps=10000, T_start=None, T_end=0.01, seed=None, restarts=1)

Find a near-optimal QUBO solution via simulated annealing.

Uses an incremental O(n)-per-step energy update and a geometric cooling schedule. Suitable for problems where solve_brute is impractical (n > 20).

Parameters:

Name	Type	Description	Default
qubo	QUBOResult	The QUBO problem to solve.	required
n_steps	int	Number of single-bit-flip proposals per restart. Default 10 000.	10000
T_start	float or None	Initial temperature. If None (default), auto-set to max(|Q|) * n, which is a reasonable scale for most problems.	None
T_end	float	Final temperature. Default 0.01.	0.01
seed	int or None	Random seed for reproducibility.	None
restarts	int	Number of independent restarts. The best result is returned. Default 1.	1

Returns:

Type	Description
SolveResult	n_evaluated is restarts * n_steps (proposals, not full evaluations).

Examples:

>>> from quprep.qubo.problems.maxcut import max_cut
>>> from quprep.qubo.solver import solve_sa
>>> import numpy as np
>>> adj = np.array([[0,1,1],[1,0,1],[1,1,0]], dtype=float)
>>> sol = solve_sa(max_cut(adj), seed=0)
>>> sol.energy
-1.0

Source code in quprep/qubo/solver.py

def solve_sa(
    qubo,
    n_steps: int = 10_000,
    T_start: float | None = None,
    T_end: float = 0.01,
    seed: int | None = None,
    restarts: int = 1,
) -> SolveResult:
    """
    Find a near-optimal QUBO solution via simulated annealing.

    Uses an incremental O(n)-per-step energy update and a geometric
    cooling schedule. Suitable for problems where ``solve_brute`` is
    impractical (n > 20).

    Parameters
    ----------
    qubo : QUBOResult
        The QUBO problem to solve.
    n_steps : int
        Number of single-bit-flip proposals per restart. Default 10 000.
    T_start : float or None
        Initial temperature. If None (default), auto-set to
        ``max(|Q|) * n``, which is a reasonable scale for most problems.
    T_end : float
        Final temperature. Default 0.01.
    seed : int or None
        Random seed for reproducibility.
    restarts : int
        Number of independent restarts. The best result is returned.
        Default 1.

    Returns
    -------
    SolveResult
        ``n_evaluated`` is ``restarts * n_steps`` (proposals, not full
        evaluations).

    Examples
    --------
    >>> from quprep.qubo.problems.maxcut import max_cut
    >>> from quprep.qubo.solver import solve_sa
    >>> import numpy as np
    >>> adj = np.array([[0,1,1],[1,0,1],[1,1,0]], dtype=float)
    >>> sol = solve_sa(max_cut(adj), seed=0)
    >>> sol.energy
    -1.0
    """
    Q = qubo.Q
    n = Q.shape[0]

    if T_start is None:
        max_q = np.max(np.abs(Q))
        T_start = float(max_q * n) if max_q > 0 else 1.0

    # Symmetric Q for O(n) delta computation: delta = flip * (Q_sym[k,:] @ x) + Q[k,k]
    Q_sym = Q + Q.T  # Q_sym[i,j] = Q[i,j] + Q[j,i]; diagonal = 2*Q[i,i]

    rng = np.random.default_rng(seed)
    alpha = (T_end / T_start) ** (1.0 / max(n_steps - 1, 1))

    best_x = np.zeros(n)
    best_energy = float("inf")

    for _ in range(restarts):
        x = rng.integers(0, 2, size=n).astype(float)
        fields = Q_sym @ x  # fields[k] = Q_sym[k,:] @ x
        energy = float(x @ Q @ x) + qubo.offset

        T = T_start
        for _ in range(n_steps):
            k = int(rng.integers(0, n))
            flip = 1.0 - 2.0 * x[k]
            delta = flip * fields[k] + Q[k, k]

            if delta < 0.0 or rng.random() < np.exp(-delta / T):
                fields += flip * Q_sym[:, k]
                x[k] = 1.0 - x[k]
                energy += delta

            T *= alpha

        if energy < best_energy:
            best_energy = energy
            best_x = x.copy()

    return SolveResult(x=best_x, energy=best_energy, n_evaluated=restarts * n_steps)

---

QAOA

quprep.qubo.qaoa.qaoa_circuit(qubo, p=1, gamma=None, beta=None)

Generate a QAOA circuit for a QUBO problem as OpenQASM 3.0 string.

The QUBO is first converted to Ising form (h, J) via qubo.to_ising(), then a p-layer QAOA ansatz is constructed.

Parameters:

Name	Type	Description	Default
qubo	QUBOResult	The QUBO problem. Use any of the problem library functions (max_cut, knapsack, tsp, portfolio) or to_qubo().	required
p	int	Number of QAOA layers (circuit depth). Default is 1.	1
gamma	list of float	Cost unitary angles, length p. Defaults to [pi/4] * p. Optimal values depend on the problem; use a classical optimizer (e.g. scipy.optimize.minimize) to tune them.	None
beta	list of float	Mixer unitary angles, length p. Defaults to [pi/8] * p.	None

Returns:

Type	Description
str	OpenQASM 3.0 circuit string. Can be run directly on Qiskit, Cirq, or any QASM-compatible backend.

Examples:

>>> from quprep.qubo.problems.maxcut import max_cut
>>> from quprep.qubo.qaoa import qaoa_circuit
>>> import numpy as np
>>> adj = np.array([[0,1,1],[1,0,1],[1,1,0]], dtype=float)
>>> qasm = qaoa_circuit(max_cut(adj), p=2)
>>> print(qasm[:60])
OPENQASM 3.0;
include "stdgates.inc";

Source code in quprep/qubo/qaoa.py

def qaoa_circuit(
    qubo,
    p: int = 1,
    gamma: list[float] | np.ndarray | None = None,
    beta: list[float] | np.ndarray | None = None,
) -> str:
    """
    Generate a QAOA circuit for a QUBO problem as OpenQASM 3.0 string.

    The QUBO is first converted to Ising form (h, J) via ``qubo.to_ising()``,
    then a p-layer QAOA ansatz is constructed.

    Parameters
    ----------
    qubo : QUBOResult
        The QUBO problem. Use any of the problem library functions
        (max_cut, knapsack, tsp, portfolio) or to_qubo().
    p : int
        Number of QAOA layers (circuit depth). Default is 1.
    gamma : list of float, optional
        Cost unitary angles, length p. Defaults to [pi/4] * p.
        Optimal values depend on the problem; use a classical optimizer
        (e.g. scipy.optimize.minimize) to tune them.
    beta : list of float, optional
        Mixer unitary angles, length p. Defaults to [pi/8] * p.

    Returns
    -------
    str
        OpenQASM 3.0 circuit string. Can be run directly on Qiskit, Cirq,
        or any QASM-compatible backend.

    Examples
    --------
    >>> from quprep.qubo.problems.maxcut import max_cut
    >>> from quprep.qubo.qaoa import qaoa_circuit
    >>> import numpy as np
    >>> adj = np.array([[0,1,1],[1,0,1],[1,1,0]], dtype=float)
    >>> qasm = qaoa_circuit(max_cut(adj), p=2)
    >>> print(qasm[:60])
    OPENQASM 3.0;
    include "stdgates.inc";
    """
    if gamma is None:
        gamma = [np.pi / 4] * p
    if beta is None:
        beta = [np.pi / 8] * p

    gamma = list(gamma)
    beta = list(beta)
    if len(gamma) != p or len(beta) != p:
        raise ValueError(f"gamma and beta must each have length p={p}.")

    # Convert QUBO -> Ising
    ising = qubo.to_ising()
    h = ising.h
    J = ising.J
    n = len(h)

    lines: list[str] = []
    lines.append("OPENQASM 3.0;")
    lines.append('include "stdgates.inc";')
    lines.append(f"qubit[{n}] q;")
    lines.append(f"bit[{n}] c;")
    lines.append("")
    lines.append("// Initialize |+>^n")
    for i in range(n):
        lines.append(f"h q[{i}];")

    for layer in range(p):
        g = gamma[layer]
        b = beta[layer]
        lines.append("")
        lines.append(f"// Layer {layer}: cost unitary U_C(gamma={g:.6f})")

        # Single-qubit Z terms: exp(-i*gamma*h_i*Z_i) = RZ(2*gamma*h_i)
        for i in range(n):
            if abs(h[i]) > 1e-12:
                angle = 2.0 * g * h[i]
                lines.append(f"rz({angle:.6f}) q[{i}];")

        # Two-qubit ZZ terms: exp(-i*gamma*J_ij*Z_i*Z_j)
        # = CNOT_{i,j} RZ(2*gamma*J_ij)_j CNOT_{i,j}
        for i in range(n):
            for j in range(i + 1, n):
                if abs(J[i, j]) > 1e-12:
                    angle = 2.0 * g * J[i, j]
                    lines.append(f"cx q[{i}], q[{j}];")
                    lines.append(f"rz({angle:.6f}) q[{j}];")
                    lines.append(f"cx q[{i}], q[{j}];")

        lines.append("")
        lines.append(f"// Layer {layer}: mixer U_B(beta={b:.6f})")
        for i in range(n):
            lines.append(f"rx({2.0 * b:.6f}) q[{i}];")

    lines.append("")
    lines.append("// Measurement")
    lines.append("c = measure q;")

    return "\n".join(lines)

---

Constraints

quprep.qubo.constraints.equality_penalty(A, b, penalty)

Encode linear equality constraints \(Ax = b\) as a QUBO penalty matrix.

Each row of \(A\) defines one constraint \(a^T x = b_i\), penalised as:

\[\lambda (a^T x - b_i)^2\]

Expanding and collecting QUBO terms:

\[ Q_{ii} \mathrel{+}= \lambda(a_i^2 - 2b_i a_i), \quad Q_{ij} \mathrel{+}= 2\lambda a_i a_j \; (i<j), \quad \text{offset} \mathrel{+}= \lambda b_i^2 \]

Parameters:

Name	Type	Description	Default
A	(ndarray, shape(m, n) or (n,))	Constraint matrix. Each row is one constraint. A 1-D array is treated as a single constraint.	required
b	(ndarray, shape(m) or scalar)	RHS vector. A scalar is broadcast across all constraints.	required
penalty	float	Lagrange multiplier \(\lambda\). Must be large enough to enforce constraints; a common heuristic is 10x the largest cost coefficient.	required

Returns:

Name	Type	Description
Q_penalty	(ndarray, shape(n, n))	Upper-triangular QUBO penalty matrix to add to the cost matrix.
offset	float	Constant offset contributed by the penalty terms.

Raises:

Type	Description
ValueError	If the number of rows in A does not match the length of b.

Source code in quprep/qubo/constraints.py

def equality_penalty(A: np.ndarray, b: np.ndarray, penalty: float) -> tuple[np.ndarray, float]:
    r"""
    Encode linear equality constraints $Ax = b$ as a QUBO penalty matrix.

    Each row of $A$ defines one constraint $a^T x = b_i$, penalised as:

    $$\lambda (a^T x - b_i)^2$$

    Expanding and collecting QUBO terms:

    $$
    Q_{ii} \mathrel{+}= \lambda(a_i^2 - 2b_i a_i), \quad
    Q_{ij} \mathrel{+}= 2\lambda a_i a_j \; (i<j), \quad
    \text{offset} \mathrel{+}= \lambda b_i^2
    $$

    Parameters
    ----------
    A : np.ndarray, shape (m, n) or (n,)
        Constraint matrix. Each row is one constraint.
        A 1-D array is treated as a single constraint.
    b : np.ndarray, shape (m,) or scalar
        RHS vector. A scalar is broadcast across all constraints.
    penalty : float
        Lagrange multiplier $\lambda$. Must be large enough to enforce
        constraints; a common heuristic is 10x the largest cost coefficient.

    Returns
    -------
    Q_penalty : np.ndarray, shape (n, n)
        Upper-triangular QUBO penalty matrix to add to the cost matrix.
    offset : float
        Constant offset contributed by the penalty terms.

    Raises
    ------
    ValueError
        If the number of rows in ``A`` does not match the length of ``b``.
    """
    A = np.atleast_2d(np.asarray(A, dtype=float))
    b = np.atleast_1d(np.asarray(b, dtype=float))
    if b.shape[0] == 1:
        b = np.broadcast_to(b, (A.shape[0],))

    m, n = A.shape
    if b.shape[0] != m:
        raise ValueError(f"A has {m} rows but b has {b.shape[0]} elements.")

    Q = np.zeros((n, n))
    offset = 0.0

    for k in range(m):
        a = A[k]
        bk = b[k]
        # Diagonal: penalty * (a_i^2 - 2*bk*a_i)
        np.fill_diagonal(Q, Q.diagonal() + penalty * (a ** 2 - 2.0 * bk * a))
        # Off-diagonal (upper-triangular): penalty * 2 * a_i * a_j
        outer = np.outer(a, a)
        Q += penalty * 2.0 * np.triu(outer, k=1)
        # Constant
        offset += penalty * bk ** 2

    return Q, offset

quprep.qubo.constraints.inequality_penalty(A, b, penalty)

Encode linear inequality constraints \(Ax \leq b\) as a QUBO penalty matrix.

Each constraint \(a^T x \leq b_i\) is converted to an equality by introducing \(K = \lceil \log_2(\text{max\_slack}+1) \rceil\) binary slack variables \(z_0, \ldots, z_{K-1}\):

\[a^T x + \sum_{k=0}^{K-1} 2^k z_k = b_i\]

This expands the variable space from \(n\) to \(n + n_{\text{slack}}\).

Parameters:

Name	Type	Description	Default
A	(ndarray, shape(m, n) or (n,))	Constraint matrix. Each row is one constraint.	required
b	(ndarray, shape(m) or scalar)	RHS values. Must satisfy \(b_i \geq \min(a^T x)\) for feasibility.	required
penalty	float	Lagrange multiplier \(\lambda\).	required

Returns:

Name	Type	Description
Q_penalty	(ndarray, shape(n + n_slack, n + n_slack))	Upper-triangular QUBO penalty matrix. Rows/columns \(0 \ldots n-1\) correspond to original variables; \(n \ldots n+n_{\text{slack}}-1\) are slack variables.
offset	float	Constant offset.
n_slack	int	Number of slack binary variables added.

Raises:

Type	Description
ValueError	If the number of rows in A does not match the length of b, or if any constraint is detected as infeasible (b_i < min(a^T x)).

Notes

The slack encoding covers slack values \(0 \ldots 2^K - 1\). If \(\text{max\_slack}\) is not a power-of-two minus one, some slack assignments are infeasible but are penalised naturally by the equality term.

Source code in quprep/qubo/constraints.py

def inequality_penalty(
    A: np.ndarray,
    b: np.ndarray,
    penalty: float,
) -> tuple[np.ndarray, float, int]:
    r"""
    Encode linear inequality constraints $Ax \leq b$ as a QUBO penalty matrix.

    Each constraint $a^T x \leq b_i$ is converted to an equality by introducing
    $K = \lceil \log_2(\text{max\_slack}+1) \rceil$ binary slack variables
    $z_0, \ldots, z_{K-1}$:

    $$a^T x + \sum_{k=0}^{K-1} 2^k z_k = b_i$$

    This expands the variable space from $n$ to $n + n_{\text{slack}}$.

    Parameters
    ----------
    A : np.ndarray, shape (m, n) or (n,)
        Constraint matrix. Each row is one constraint.
    b : np.ndarray, shape (m,) or scalar
        RHS values. Must satisfy $b_i \geq \min(a^T x)$ for feasibility.
    penalty : float
        Lagrange multiplier $\lambda$.

    Returns
    -------
    Q_penalty : np.ndarray, shape (n + n_slack, n + n_slack)
        Upper-triangular QUBO penalty matrix. Rows/columns $0 \ldots n-1$
        correspond to original variables; $n \ldots n+n_{\text{slack}}-1$
        are slack variables.
    offset : float
        Constant offset.
    n_slack : int
        Number of slack binary variables added.

    Raises
    ------
    ValueError
        If the number of rows in ``A`` does not match the length of ``b``,
        or if any constraint is detected as infeasible
        (``b_i < min(a^T x)``).

    Notes
    -----
    The slack encoding covers slack values $0 \ldots 2^K - 1$. If
    $\text{max\_slack}$ is not a power-of-two minus one, some slack assignments
    are infeasible but are penalised naturally by the equality term.
    """
    A = np.atleast_2d(np.asarray(A, dtype=float))
    b = np.atleast_1d(np.asarray(b, dtype=float))
    if b.shape[0] == 1:
        b = np.broadcast_to(b, (A.shape[0],))

    m, n = A.shape
    if b.shape[0] != m:
        raise ValueError(f"A has {m} rows but b has {b.shape[0]} elements.")

    # Determine slack bits needed per constraint
    slack_specs: list[tuple[int, int, list[int]]] = []  # (start, K, powers)
    total_slack = 0
    for k in range(m):
        a = A[k]
        bk = b[k]
        # Minimum possible value of a^T x (binary x)
        min_val = float(np.sum(a[a < 0]))
        max_slack = int(np.ceil(bk - min_val))
        if max_slack < 0:
            raise ValueError(
                f"Constraint {k} appears infeasible: "
                f"b={bk:.4f} < min(a^T x)={min_val:.4f}."
            )
        K = max(1, int(np.ceil(np.log2(max_slack + 1 + 1e-10))))
        powers = [2 ** j for j in range(K)]
        slack_specs.append((total_slack, K, powers))
        total_slack += K

    n_slack = total_slack
    N = n + n_slack
    Q = np.zeros((N, N))
    offset = 0.0

    for k in range(m):
        a = A[k]
        bk = b[k]
        start, K, powers = slack_specs[k]

        # Augmented constraint: [a | 2^0 2^1 ... 2^{K-1}] * [x | z] = b
        a_aug = np.zeros(N)
        a_aug[:n] = a
        for j, p in enumerate(powers):
            a_aug[n + start + j] = float(p)

        Q_pen, off = equality_penalty(a_aug.reshape(1, -1), np.array([bk]), penalty)
        Q += Q_pen
        offset += off

    return Q, offset, n_slack

---

Utilities

quprep.qubo.utils.add_qubo(q1, q2, weight=1.0)

Add two QUBO problems of the same size.

Useful for multi-objective problems where you want to combine an objective term (e.g. max_cut) with a constraint term (e.g. equality penalty) that was built separately.

Parameters:

Name	Type	Description	Default
q1	QUBOResult	First QUBO.	required
q2	QUBOResult	Second QUBO. Must have the same number of variables as q1.	required
weight	float	Scalar multiplier applied to q2 before addition. Default is 1.0.	1.0

Returns:

Type	Description
QUBOResult	Combined QUBO with Q = q1.Q + weight * q2.Q and merged variable maps.

Raises:

Type	Description
ValueError	If q1 and q2 have different Q matrix shapes.

Examples:

>>> from quprep.qubo.problems.maxcut import max_cut
>>> from quprep.qubo.constraints import equality_penalty
>>> from quprep.qubo.converter import QUBOResult
>>> from quprep.qubo.utils import add_qubo
>>> import numpy as np
>>> q_cut = max_cut(np.array([[0,1],[1,0]], dtype=float))
>>> Q_pen, off = equality_penalty(np.array([[1.0, 1.0]]), np.array([1.0]), 5.0)
>>> q_pen = QUBOResult(Q=Q_pen, offset=off)
>>> combined = add_qubo(q_cut, q_pen)

Source code in quprep/qubo/utils.py

def add_qubo(q1, q2, weight: float = 1.0):
    """
    Add two QUBO problems of the same size.

    Useful for multi-objective problems where you want to combine an
    objective term (e.g. max_cut) with a constraint term (e.g. equality
    penalty) that was built separately.

    Parameters
    ----------
    q1 : QUBOResult
        First QUBO.
    q2 : QUBOResult
        Second QUBO. Must have the same number of variables as q1.
    weight : float
        Scalar multiplier applied to q2 before addition. Default is 1.0.

    Returns
    -------
    QUBOResult
        Combined QUBO with ``Q = q1.Q + weight * q2.Q`` and merged variable maps.

    Raises
    ------
    ValueError
        If q1 and q2 have different Q matrix shapes.

    Examples
    --------
    >>> from quprep.qubo.problems.maxcut import max_cut
    >>> from quprep.qubo.constraints import equality_penalty
    >>> from quprep.qubo.converter import QUBOResult
    >>> from quprep.qubo.utils import add_qubo
    >>> import numpy as np
    >>> q_cut = max_cut(np.array([[0,1],[1,0]], dtype=float))
    >>> Q_pen, off = equality_penalty(np.array([[1.0, 1.0]]), np.array([1.0]), 5.0)
    >>> q_pen = QUBOResult(Q=Q_pen, offset=off)
    >>> combined = add_qubo(q_cut, q_pen)
    """
    from quprep.qubo.converter import QUBOResult

    if q1.Q.shape != q2.Q.shape:
        raise ValueError(
            f"Q shapes must match: {q1.Q.shape} vs {q2.Q.shape}. "
            "Use to_qubo(cost, constraints=[...]) to build combined QUBOs directly."
        )

    Q = q1.Q + weight * q2.Q
    offset = q1.offset + weight * q2.offset
    n_original = q1.n_original

    # Merge variable maps — q1 takes precedence on conflict
    var_map = {**q2.variable_map, **q1.variable_map}

    return QUBOResult(Q=Q, offset=offset, variable_map=var_map, n_original=n_original)

quprep.qubo.visualize.draw_qubo(qubo, title='QUBO Matrix', cmap='RdBu_r', ax=None)

Draw a heatmap of the QUBO Q matrix.

Positive entries (coupling penalties) are shown in one colour, negative entries (incentives) in another. The diagonal encodes linear terms; off-diagonal entries encode quadratic interactions.

Parameters:

Name	Type	Description	Default
qubo	QUBOResult	The QUBO problem to visualize.	required
title	str	Plot title. Default is "QUBO Matrix".	'QUBO Matrix'
cmap	str	Matplotlib colormap. Default is "RdBu_r" (blue=negative, red=positive).	'RdBu_r'
ax	Axes	Existing axes to draw on. Creates a new figure if None.	None

Returns:

Type	Description
Axes	Axes object with the Q matrix heatmap drawn on it.

Raises:

Type	Description
ImportError	If matplotlib is not installed. Install with: pip install quprep[viz]

quprep.qubo.visualize.draw_ising(ising, title='Ising Model', ax=None)

Draw the Ising model as a weighted graph.

Nodes are arranged in a circle. Node colour represents the bias (h_i): blue = negative bias (prefers s=-1), red = positive (prefers s=+1). Edge colour and thickness represent coupling strength (J_ij).

Parameters:

Name	Type	Description	Default
ising	IsingResult	The Ising model to visualize.	required
title	str	Plot title.	'Ising Model'
ax	Axes	Existing axes. Creates a new figure if None.	None

Returns:

Type	Description
Axes	Axes object with the Ising graph drawn on it.

Raises:

Type	Description
ImportError	If matplotlib is not installed.