API reference

MFAOA

A solver for binary quadratic programs without constraints (UBQP in the operations-research literature, QUBO elsewhere), built on the mean-field approximate optimization algorithm of Misra-Spieldenner, Bode, Schuhmacher, Stollenwerk, Bagrets and Wilhelm, PRX Quantum 4, 030335 (2023).

The solver minimizes (or maximizes) a quadratic function of binary variables. It accepts the problem either as a QUBO, x' * Q * x + c' * x + constant over x ∈ {0,1}ⁿ (QUBOProblem), or in the Ising form the algorithm works with natively, h' * σ + Σ_{i<j} J[i,j] σ[i] σ[j] over σ ∈ {-1,+1}ⁿ (IsingProblem). The matrices can be dense or sparse. Call solve with a Schedule, or with the number of steps p and the step size τ of the default linear schedule. From JuMP, use Model(MFAOA.Optimizer) with binary variables and a quadratic objective.

QUBOProblem(Q; c = zeros(n), constant = 0, sense = :min)

A binary quadratic program in QUBO form: minimize (or, with sense = :max, maximize)

x' * Q * x + c' * x + constant      over x[i] ∈ {0, 1}.

Q is used exactly as given: a symmetric matrix counts each off-diagonal pair twice, an upper-triangular matrix once. Diagonal entries are allowed; since x[i]^2 = x[i] they act as linear terms. Q can be dense or sparse, and is kept without copying when its element type is already the problem's. solve converts the problem to Ising form with x = (1 + σ) / 2 (see IsingProblem(::QUBOProblem)) and reports the solution both as x and as σ.

IsingProblem(h, J; Δ = nothing, offset = 0, sense = :min, scale = nothing)

A binary quadratic program in Ising form, in variables σ[i] that are -1 or +1: minimize (or, with sense = :max, maximize)

h' * σ + Σ_{i<j} J[i,j] σ[i] σ[j] + offset.

J must be symmetric with a zero diagonal. Since σ[i]^2 = 1, a diagonal entry would only add a constant; put it into offset instead. J can be a dense Matrix, a SparseMatrixCSC, which is much faster when most entries are zero, or any other matrix that LinearAlgebra.mul! can multiply with a vector, such as a Symmetric view. h and J are used as they are, without a copy, when their element type is already the problem's. Changing them afterwards changes the problem, except for its fieldscale, which is computed when the problem is built.

Two keywords tune the algorithm and rarely need to be set:

• Δ: one positive factor per variable on the angles β of the Schedule. The default, nothing, means 1 for every variable and takes a faster code path.
• scale: the unit in which the schedule measures the slope of the objective; see fieldscale. The default makes the solution independent of the units of h and J; scale = 1 uses those units themselves.

IsingProblem(prob::QUBOProblem; Δ = nothing, scale = nothing)

The same problem in Ising form, with x = (1 + σ) / 2. A pair term Q[i,j] x[i] x[j] + Q[j,i] x[j] x[i] becomes the coupling J[i,j] = (Q[i,j] + Q[j,i]) / 4 and contributes the same amount to both fields and to the offset. The weight of x[i] alone, c[i] + Q[i,i], goes half into h[i] and half into the offset. J is sparse if Q is.

objective(prob, σ)

The objective value of prob at σ, a vector with entries -1 and +1, including the offset.

objective(prob::QUBOProblem, x)

The objective value x' * Q * x + c' * x + constant at x (entries 0 or 1). Maximization problems report the value as stated, without sign change.

isflipsymmetric(prob) -> Bool

Whether prob has a group of variables that MFAOA cannot solve for: flipping all variables of the group at once (σ → -σ, or x → 1 - x in a QUBOProblem) leaves the objective unchanged. In such a group, no term of the objective links a variable to one outside the group, and h is zero on every variable of the group, up to rounding errors relative to the largest coupling. When h is zero everywhere, as in Max-Cut, the whole problem is such a group. A single variable that appears in no term of the objective does not count: any value is optimal for it, and the solution sets it to σ = +1.

The algorithm starts every variable at σ = 0, halfway between -1 and +1. There the objective has zero slope in the variables of such a group, so they would never move. solve throws an ArgumentError instead, which names the groups. Fix one variable of each group to either value and solve for the others. This loses nothing: flipping the group turns every solution into one with the same objective value, so of each such pair, one solution remains.

To fix σ[i] = s, remove variable i from the problem, add s * J[j,i] to h[j] for every other variable j, and add s * h[i] to the offset. In a QUBOProblem, fixing x[i] = v adds v * (Q[i,j] + Q[j,i]) to c[j] and v * (Q[i,i] + c[i]) to the constant.

fieldscale(prob::IsingProblem) -> s

The unit in which the schedule measures the slope of the objective of prob.

The algorithm moves each variable according to the slope of the objective in that variable, m[i] = h[i] + Σⱼ J[i,j] z[j], where z[j], between -1 and +1, is the current value of variable j. In step k, the angle that this slope sets is 2γ[k] m[i] / s (see Schedule). Unless scale was given when the problem was built, s is the root mean square of m[i] over all variables i and over random values ±1 of z:

s = sqrt((Σᵢ hᵢ² + Σᵢⱼ Jᵢⱼ²) / n).

Multiplying h and J by a constant multiplies s by the same constant, so the algorithm, and the solution, do not depend on the units of the objective.

Schedule(β, γ)

The angles that drive the algorithm, one pair β[k], γ[k] per step k.

The algorithm represents each variable σ[i] by a unit vector in three dimensions. The vector's z component, between -1 and +1, is the variable's current value. Every vector starts at (1, 0, 0), so every variable starts at 0, halfway between -1 and +1. Each step k turns every vector i twice:

1. about the z axis by the angle 2γ[k] m[i] / s, where m[i] is the slope of the objective in variable i at the current values, and s is the problem's fieldscale (the opposite angle for a maximization);
2. about the x axis by the angle 2β[k] Δ[i], where Δ[i] is 1 unless the problem sets it.

After the last step, σ[i] is the sign of the z component. The two vectors β and γ must have the same length, the number of steps p. linear_schedule builds the usual choice; any other pair of vectors can be passed directly.

Schedule(β::AbstractVector, γ::AbstractVector)

Build a schedule from any two real vectors, promoted to a common floating-point type.

linear_schedule(p, τ = 0.5; T = Float64)

The usual schedule: over p steps, the angles γ grow and the angles β shrink linearly,

γ[k] = τ k / p,    β[k] = τ (1 - (k - 1) / p),    k = 1, …, p.

The schedule is Eq. 13 of Misra-Spieldenner et al., PRX Quantum 4, 030335 (2023). τ sets the size of the angles. Because γ is measured in units of the problem's fieldscale, the same τ suits problems of any size and units; values between 0.25 and 0.5 work well on typical instances. With more steps, the angles change more slowly from one step to the next. That usually gives a better solution, and the time of a solve grows in proportion to p.

solve(prob, sched::Schedule) -> Result
solve(prob; p = 1000, τ = 0.5) -> Result

Solve the IsingProblem or QUBOProblem prob with the schedule sched, or with the linear_schedule of p steps and step size τ. The run is deterministic. See solve! for the allocation-free form and Result for what comes back.

Result

The result of solve, with the fields

• σ::Vector{Int8}: the solution in Ising form, entries -1 and +1.
• x::Vector{Int8}: the same solution in QUBO form, x = (1 + σ) / 2, entries 0 and 1.
• objective: the objective value of the solution, in the problem's sense and including its constant.
• objective_mf: the objective at the values between -1 and +1 that the algorithm ends with, before they are rounded to σ. It equals objective when every value has reached -1 or +1. A large gap means that variables were still undecided at the end; more steps usually close it.
• state::SpinState: the final vectors that represent the variables (see Schedule).

solve!(ws::Workspace, prob::IsingProblem, sched::Schedule) -> ws

Run the algorithm on prob with the schedule sched, using only the memory in ws. Afterwards ws.state holds the final vectors that represent the variables, ws.m the slopes of the objective there, and ws.σ the solution. Allocates nothing.

Throws an ArgumentError when isflipsymmetric(prob), naming the variables that the algorithm cannot move, or when ws was made for a different number of variables.

Workspace(prob)

All memory that solve! needs for prob: the vectors state that represent the variables, the slopes m of the objective in the variables, the solution σ, and group, which the check of the problem before a solve uses (see isflipsymmetric). Allocate it once and reuse it across solves of problems with the same number of variables.

evolve!(ws::Workspace, prob::IsingProblem, sched::Schedule) -> ws

Run every step of sched with step!, starting from the vectors in ws.state as they are. Allocates nothing.

step!(ws::Workspace, prob::IsingProblem, β_k, γ_k) -> ws

Apply one step of the algorithm (see Schedule) to the vectors in ws.state. Compute the slope m[i] of the objective in every variable into ws.m. Then turn every vector i about the z axis by the angle 2γ_k * m[i] / s, where s is the problem's fieldscale (the opposite angle for a maximization), and about the x axis by 2β_k * Δ[i]. Allocates nothing.

SpinState{T}(n)

The unit vectors in three dimensions by which the algorithm represents n variables (see Schedule), stored as three vectors x, y, z of their components. The z component of a vector is the current value of its variable, between -1 and +1. A fresh state has every vector at (1, 0, 0), where the algorithm starts: every variable at 0.

MFAOA.Optimizer(; p = 1000, tau = 0.5)

The MFAOA solver for JuMP (as a MathOptInterface optimizer):

using JuMP, MFAOA
model = Model(MFAOA.Optimizer)
@variable(model, x[1:3], Bin)
@objective(model, Min, 2x[1]*x[2] - 3x[2]*x[3] - x[1] - x[2] + x[3])
optimize!(model)
value.(x), objective_value(model)

The model may contain binary variables, a linear or quadratic objective in either sense, and fixed variables (fix(x[i], 1)). Bounds on a variable are accepted when they allow 0 and 1, as in 0 <= x <= 1, and integer variables when their bounds allow only 0 and 1. A bound that allows only one of the two values holds the variable there, as fix would. Nothing else is accepted. Options, set with set_attribute(model, name, value):

• "p": the number of steps of the linear schedule (default 1000);
• "tau": its step size (default 0.5);
• "beta", "gamma": an explicit Schedule as two vectors of equal length, used instead of the linear schedule when both are set.

A solve ends with LOCALLY_SOLVED: the solver returns a solution, without a proof that no better one exists. INVALID_MODEL means that MFAOA cannot solve the model, and raw_status(model) says why and what to do. Most often, the objective does not change when every variable x is replaced by 1 - x, as in Max-Cut; fixing one variable, say with fix(x[1], 0), solves that.