Optimization

eigenkit.optimization solves combinatorial problems. In these problems, the answer is an assignment of bits, and the cost is a quadratic function of the bits. You write the problem as a QUBO (quadratic unconstrained binary optimisation) or an Ising model, and pass it to QAOA (the Quantum Approximate Optimization Algorithm). You get back the best bitstring that it found, together with the distribution that it sampled from.

Like chemistry and reservoir, this is a domain package of the core. It uses only NumPy and SciPy, without PyTorch or the package of any backend. QAOA also uses provider="numpy" by default, so nothing here needs an optional dependency.

The theory is on the QAOA page of Algorithms Explained.

MaxCut in five lines

MaxCut looks for the split of the nodes of a graph that cuts the most edges. max_cut builds the Ising model for you, written as an energy to minimise:

from eigenkit.optimization import QAOA, max_cut

problem = max_cut([(0, 1), (1, 2), (2, 0), (2, 3)])
result = QAOA(layers=2, random_state=0).solve(problem)

print(result.best_bitstring)  # (1, 1, 0, 1)
print(result.best_energy)  # -3.0
print(result.success)  # True

The graph has four nodes and four edges, and the best cut has three edges. On a problem this small, you can check the answer by trying every bitstring. It is worth doing once, to see what "the answer" means here:

import itertools

best = min(itertools.product([0, 1], repeat=4), key=problem.energy)

print(best)  # (0, 0, 1, 0)
print(problem.energy(best))  # -3.0

The bitstring is different, but the energy is the same, and each bit of one bitstring is the opposite of the same bit in the other. Moving every node to the other side of the cut describes the same partition, so the solutions of MaxCut always come in pairs. Compare energies, not bitstrings.

Stating the problem

QUBO minimises xᵀQx + constant over x ∈ {0,1}ⁿ. energy gives the cost of any candidate. This is the quickest way to check a formulation before you spend time on a solve:

import numpy as np
from eigenkit.optimization import QUBO

Q = np.array([[-2.0, 1.0], [0.0, -2.0]])
qubo = QUBO(Q, constant=0.5)

print(qubo.energy([0, 0]))  # 0.5
print(qubo.energy([1, 0]))  # -1.5
print(qubo.energy([1, 1]))  # -2.5   the minimum

print(QAOA(layers=2, random_state=0).solve(qubo).best_bitstring)  # (1, 1)

Ising is the same problem in spin variables. It minimises Σᵢhᵢsᵢ + Σᵢ<ⱼJᵢⱼsᵢsⱼ + constant over s ∈ {−1,+1}ⁿ. The two forms are equivalent, and QUBO.to_ising() converts a QUBO to an Ising model:

from eigenkit.optimization import Ising

ising = qubo.to_ising()
print(ising.h)  # [0.75 0.75]
print(ising.constant)  # -1.25

direct = Ising(np.array([0.5, -0.5]), np.array([[0.0, 1.0], [0.0, 0.0]]))
print(direct.energy([0, 0]))  # 1.0
print(direct.energy([1, 0]))  # -2.0

Note that Ising.energy takes bits, not spins: a 0 means the spin +1, and a 1 means the spin −1. This keeps energy the same for both types of problem, but it will surprise you if you pass -1 directly.

max_cut(edges, weights=None, num_nodes=None) is a shortcut that builds the Ising model of a weighted graph. Pass weights for weighted edges, and num_nodes when the graph has isolated nodes that the list of edges never mentions.

Reading the result

solve returns an OptimizationResult, which holds more than the best bitstring:

Field	What it holds
`best_bitstring`	The best assignment found, as a tuple of `0` and `1` values.
`best_energy`	Its exact classical cost. This is the answer you would act on.
`expectation`	The mean cost over the whole distribution. This is the quantity that QAOA minimises.
`probabilities`	The whole distribution over all `2**n` bitstrings.
`optimal_parameters`	The optimised angles. There are `2 * layers` of them.
`n_iterations`	The number of steps of the optimiser.
`success`	Whether the optimiser reported convergence.

result = QAOA(layers=2, random_state=0).solve(problem)

print(result.best_energy)  # -3.0        the best sample
print(result.expectation)  # -2.903547   the mean over the distribution
print(result.probabilities.shape)  # (16,)       one entry per bitstring
print(result.n_iterations)  # 18
print(result.optimal_parameters)  # 4 angles: gamma and beta per layer

The two energies answer different questions. best_energy is the cost of the best single bitstring among those that the optimised state makes most likely. This is the answer you would act on. expectation is the average cost over the whole distribution, and it is the quantity that QAOA actually minimises. The difference between them comes from the circuit not putting all its amplitude on one solution. More layers make it smaller.

probabilities is the whole distribution over all 2ⁿ bitstrings, so you can also look at the solutions that come after the best one. To choose best_bitstring, solve scores the most probable bitstrings with the energy of the problem. By default, it scores the 64 most probable ones, or all of them when there are fewer. Pass top_k= to solve to change how many bitstrings are scored, which matters once n is large.

The circuit behind it

QAOA builds the cost Hamiltonian and the ansatz itself, but you can also get both separately, to inspect or reuse them. to_hamiltonian() turns a problem into the pair (Observable, offset), with the same convention as the chemistry mappings. qaoa_ansatz builds the standard circuit: |+⟩ⁿ, followed by layers repetitions of the cost unitary and the mixer.

from eigenkit import Circuit
from eigenkit.optimization import qaoa_ansatz

hamiltonian, offset = problem.to_hamiltonian()
print(len(hamiltonian.terms))  # 4
print(offset)  # -2.0
print(hamiltonian.terms[0])  # PauliTerm(coeff=0.5, paulis=(('Z', 0), ('Z', 1)))

spec = qaoa_ansatz(hamiltonian, 4, layers=2)
print(spec.n_weights)  # 4

circuit = Circuit(spec, observables=hamiltonian, provider="numpy")
print(float(circuit(result.optimal_parameters, [])) + offset)  # -2.903547

The last number is result.expectation, computed again from the parts. It is a useful check when you start to use your own mixer or a different schedule of parameters.

One limitation: qaoa_ansatz expects an Observable, which is the type that to_hamiltonian returns. It does not accept a PauliExpr built with the X, Y and Z helpers, although Circuit(observables=...) does. Passing one raises an AttributeError about terms.

API reference

`eigenkit.optimization`

Combinatorial optimisation with the core of the library.

Write a problem as an Ising or QUBO model, then solve it with QAOA, the Quantum Approximate Optimization Algorithm. The ansatz uses the gates of the core, the angles are optimised with minimize_expectation, and the solution is read from the distribution of bitstrings of the optimised state. The package uses only NumPy and SciPy, needs no PyTorch, and works with any backend.

Example

from eigenkit.optimization import max_cut, QAOA

problem = max_cut([(0, 1), (1, 2), (2, 0), (2, 3)])
result = QAOA(layers=2).solve(problem)
print(result.best_bitstring, result.best_energy)

`Ising` `dataclass`

An Ising problem, which minimises Σ h_i s_i + Σ_{i<j} J_ij s_i s_j + constant.

The minimum is taken over s ∈ {+1, -1}^n. Only the strict upper triangle of J (i < j) is used.

Attributes:

Name	Type	Description
`h`	`ndarray`	The linear field, of shape `(n,)`, with one bias per spin.
`J`	`ndarray`	The coupling matrix, of shape `(n, n)`. Only the strict upper triangle is read, so a symmetric matrix and its upper half give the same problem, and the entries on or below the diagonal are ignored.
`constant`	`float`	An energy offset added to every configuration. It shifts the reported energy without changing which configuration wins. Default: `0.0`.

`num_vars` `property`

The number of spins, read from the length of h.

`brute_force()`

Find the exact minimum by trying every assignment.

This is practical only for a small n.

Returns:

Name	Type	Description
`bitstring`	`tuple[int, ...]`	The assignment with the lowest energy.
`energy`	`float`	Its energy.

`energy(bits)`

Return the energy of one assignment.

The bits are converted to spins with s = 1 - 2b, so a 0 is the spin +1 and a 1 is the spin -1.

Parameters:

Name	Type	Description	Default
`bits`	`Iterable[int]`	One bit per variable, in the order of the variables.	required

Returns:

Type	Description
`float`	The energy of the assignment for this problem.

`to_hamiltonian()`

Return the qubit Hamiltonian of the problem, with the constant kept apart.

The Hamiltonian is Σ h_i Z_i + Σ_{i<j} J_ij Z_i Z_j.

Returns:

Name	Type	Description
`observable`	`Observable`	The Hamiltonian, as an `Observable`.
`offset`	`float`	The constant of the problem.

Raises:

Type	Description
`ValueError`	If every coefficient of the model is `0`.

`OptimizationResult` `dataclass`

Bases: ValueEquality

The result of a run of the Quantum Approximate Optimization Algorithm (QAOA).

The result is read-only. Assigning to a field raises FrozenInstanceError, and writing into an array field raises ValueError. This prevents a mistyped = (instead of ==) or an in-place edit from leaving a result whose scalar fields describe a different run from the arrays it holds. To get a writable copy of an array, use np.array(result.field).

Results compare by value, including their arrays. They are not hashable, because a result is a collection of numbers and not a key, and a correct hash of a float array is expensive to compute.

Attributes:

Name	Type	Description
`best_bitstring`	`tuple[int, ...]`	The configuration with the lowest energy that was found, as a tuple of bits.
`best_energy`	`float`	Its energy for the problem being solved. This is the answer to report: QAOA samples, so the best sample is better than the average.
`expectation`	`float`	The mean energy of the final state. Compare it with `best_energy` to see how concentrated the distribution is.
`optimal_parameters`	`ndarray`	The angles `γ` and `β` that the optimiser settled on.
`probabilities`	`ndarray`	The probability of every basis state, of length `2**n`.
`n_iterations`	`int`	The number of iterations of the optimiser.
`success`	`bool`	Whether the optimiser reported convergence.

`__post_init__()`

Make the array fields read-only.

Freezing the dataclass prevents a field from being reassigned, but it does not prevent writing into the array that the field holds, and a result that has been written into still reports the scalar values of the original run. The result therefore keeps its own copies, so the arrays it was built from can still be changed by the caller. To get a writable copy, use np.array(result.optimal_parameters).

`QAOA`

Bases: BaseQuantumSolver

The Quantum Approximate Optimization Algorithm (QAOA) solver.

It solves an Ising or QUBO problem.

Parameters:

Name	Type	Description	Default
`layers`	`int`	The depth `p` of QAOA. Default: `1`.	`1`
`provider`	`ProviderSpec`	The backend. Default: `"numpy"`, the reference NumPy simulator.	`'numpy'`
`optimizer`	`Optimizer`	A SciPy method name, `"spsa"`, or a custom callable (see `minimize_expectation`). Default: `"L-BFGS-B"`.	`'L-BFGS-B'`
`random_state`	`int \| None`	The seed of the random initial angles. With `None`, they are not seeded. Default: `42`.	`42`

`solve(problem, *, initial_parameters=None, top_k=None, **minimize_kwargs)`

Optimise the QAOA angles and read the answer from the final state.

Parameters:

Name	Type	Description	Default
`problem`	`QUBO \| Ising`	The problem to minimise.	required
`initial_parameters`	`ndarray \| None`	The initial angles `γ` and `β`, of shape `(2 * layers,)`. They are random when they are not given. Default: `None`.	`None`
`top_k`	`int \| None`	How many of the most probable bitstrings to score with the energy of the problem. The best of them is returned as `best_bitstring`, and the others are not returned. Default: `None`, which scores the `64` most probable bitstrings, or all of them when there are fewer.	`None`
`**minimize_kwargs`	`Any`	Other keyword arguments, passed to `minimize_expectation`, for example `max_iterations`.	`{}`

Returns:

Type	Description
`OptimizationResult`	An `OptimizationResult` with the best bitstring, its energy and the final distribution.

Raises:

Type	Description
`ValueError`	If `top_k` is smaller than `1`.

`QUBO` `dataclass`

A quadratic unconstrained binary optimisation (QUBO) problem.

It minimises xᵀQx + constant over x ∈ {0, 1}^n.

Attributes:

Name	Type	Description
`Q`	`ndarray`	The cost matrix, of shape `(n, n)`. The diagonal entries are the linear terms, because `x_i² = x_i` for binary variables.
`constant`	`float`	An energy offset added to every configuration. Default: `0.0`.

`num_vars` `property`

The number of binary variables, read from the shape of Q.

`brute_force()`

Find the exact minimum by trying every assignment.

This is practical only for a small n.

Returns:

Name	Type	Description
`bitstring`	`tuple[int, ...]`	The assignment with the lowest energy.
`energy`	`float`	Its energy.

`energy(bits)`

Return the energy of one assignment.

Parameters:

Name	Type	Description	Default
`bits`	`Iterable[int]`	One bit per variable, in the order of the variables.	required

Returns:

Type	Description
`float`	The energy of the assignment for this problem.

`to_hamiltonian()`

Return the qubit Hamiltonian of the problem, through to_ising().

Returns:

Name	Type	Description
`observable`	`Observable`	The Hamiltonian, as an `Observable`.
`offset`	`float`	The constant of the problem.

`to_ising()`

Convert the problem exactly to an Ising model, with x_i = (1 - s_i) / 2.

Returns:

Type	Description
`Ising`	The equivalent `Ising` model.

`max_cut(edges, weights=None, num_nodes=None)`

Build the Ising model of the MaxCut problem for a weighted graph.

The problem is written as an energy to minimise. The spin assignment with the lowest energy is a maximum cut, and the cut value of a bitstring is -ising.energy(bits).

Parameters:

Name	Type	Description	Default
`edges`	`Sequence[tuple[int, int]]`	The edges of the graph, as pairs `(i, j)`.	required
`weights`	`Sequence[float] \| None`	One weight per edge. Default: `None`, which gives every edge the weight `1`.	`None`
`num_nodes`	`int \| None`	The number of nodes. Default: `None`, which takes it from the edges.	`None`

Returns:

Type	Description
`Ising`	The Ising problem whose minimum is a maximum cut.

`qaoa_ansatz(hamiltonian, num_qubits, layers)`

Build the standard QAOA circuit.

QAOA is the Quantum Approximate Optimization Algorithm. The circuit starts in |+⟩^n, then applies layers repetitions of the cost unitary exp(-i γ_l H_C) and the mixer exp(-i β_l Σ_i X_i). It has 2 * layers weights, in the order [γ_0, β_0, γ_1, β_1, ...].

Parameters:

Name	Type	Description	Default
`hamiltonian`	`Observable`	The cost Hamiltonian.	required
`num_qubits`	`int`	The number of qubits, one per variable.	required
`layers`	`int`	The depth of QAOA, which is the number of pairs of a cost layer and a mixer layer.	required

Returns:

Type	Description
`CircuitSpec`	The ansatz, as a `CircuitSpec` with `2 * layers` trainable angles.

Optimization

MaxCut in five lines

Stating the problem

Reading the result

The circuit behind it

API reference

eigenkit.optimization

Ising dataclass

num_vars property

brute_force()

energy(bits)

to_hamiltonian()

OptimizationResult dataclass

__post_init__()

QAOA

solve(problem, *, initial_parameters=None, top_k=None, **minimize_kwargs)

QUBO dataclass

num_vars property

brute_force()

energy(bits)

to_hamiltonian()

to_ising()

max_cut(edges, weights=None, num_nodes=None)

qaoa_ansatz(hamiltonian, num_qubits, layers)

`eigenkit.optimization`

`Ising` `dataclass`

`num_vars` `property`

`brute_force()`

`energy(bits)`

`to_hamiltonian()`

`OptimizationResult` `dataclass`

`__post_init__()`

`QAOA`

`solve(problem, *, initial_parameters=None, top_k=None, **minimize_kwargs)`

`QUBO` `dataclass`

`num_vars` `property`

`brute_force()`

`energy(bits)`

`to_hamiltonian()`

`to_ising()`

`max_cut(edges, weights=None, num_nodes=None)`

`qaoa_ansatz(hamiltonian, num_qubits, layers)`