Order Partitioning for A/B Testing - QUBONotebooks: Quantum Integer Programming

Maintained by the JuliaQUBO organization
SECQUOIA · PSR Energy

An original Julia case study with a corrected grouped-risk objective and exact checks.

Setup¶

Local installation¶

From the repository root, instantiate the shared Julia environment before opening this notebook:

julia --project=notebooks_jl -e 'using Pkg; Pkg.instantiate()'

Run this notebook from a clean kernel with:

make verify-order-partitioning-julia

The default workflow is credential-free, uses no external data, and solves the six-order instance locally.

Google Colab¶

Open the badge above, select a Julia runtime, and run the setup cells. The bootstrap clones this repository only when Colab does not already have it, then activates the same checked-in notebook project used locally.

function load_qubonotebooks_bootstrap()
    candidates = (
        joinpath(pwd(), "scripts", "notebook_bootstrap.jl"),
        joinpath(pwd(), "..", "scripts", "notebook_bootstrap.jl"),
        joinpath(pwd(), "QUBONotebooks", "scripts", "notebook_bootstrap.jl"),
        joinpath("/content", "QUBONotebooks", "scripts", "notebook_bootstrap.jl"),
    )

    for candidate in candidates
        if isfile(candidate)
            include(candidate)
            return nothing
        end
    end

    in_colab = haskey(ENV, "COLAB_RELEASE_TAG") || haskey(ENV, "COLAB_JUPYTER_IP") || isdir(joinpath("/content", "sample_data"))
    if in_colab
        repo_dir = get(ENV, "QUBONOTEBOOKS_REPO_DIR", joinpath(pwd(), "QUBONotebooks"))
        if !isdir(repo_dir)
            println("[bootstrap] Cloning JuliaQUBO/QUBONotebooks into $repo_dir")
            run(Cmd(["git", "clone", "--quiet", "--depth", "1", "https://github.com/JuliaQUBO/QUBONotebooks.git", repo_dir]))
        end
        include(joinpath(repo_dir, "scripts", "notebook_bootstrap.jl"))
        return nothing
    end

    error("Could not locate scripts/notebook_bootstrap.jl from $(pwd()).")
end

load_qubonotebooks_bootstrap()

BOOTSTRAP = Base.invokelatest(QUBONotebooksBootstrap.bootstrap_notebook, "8-OrderPartitioning")
QUBONOTEBOOKS_REPO_DIR = BOOTSTRAP.repo_dir
JULIA_NOTEBOOKS_DIR = BOOTSTRAP.notebooks_dir
JULIA_PROJECT_DIR = BOOTSTRAP.project_dir
IN_COLAB = BOOTSTRAP.in_colab;

import Pkg

python_warning_filter = "ignore:invalid escape sequence:SyntaxWarning"
python_warning_filters = String.(filter(!isempty, split(get(ENV, "PYTHONWARNINGS", ""), ",")))
if python_warning_filter ∉ python_warning_filters
    push!(python_warning_filters, python_warning_filter)
    ENV["PYTHONWARNINGS"] = join(python_warning_filters, ",")
end

if @isdefined(JULIA_PROJECT_DIR)
    Pkg.activate(JULIA_PROJECT_DIR; io = devnull)
else
    Pkg.activate(@__DIR__; io = devnull)
end
Pkg.instantiate(; io = devnull, allow_autoprecomp = false)

Learning objectives¶

By the end of this notebook you will be able to:

formulate value and factor-risk balance as a QUBO with explicit weights;
explain why each risk-factor square must enclose the sum over orders;
verify all 64 assignments independently of the optimizer;
decode group values and risk exposures from an exact optimum; and
explain complement symmetry and the trade-off induced by the two weights.

Prerequisites¶

Prior notebooks: Notebook 2 introduces JuMP and QUBO models; Notebook 7 develops exhaustive checking helpers. This notebook restates the helpers it needs so it remains useful on its own.

Mathematical background: Binary variables, finite sums, squared deviations, and basic A/B testing terminology.

Software: Julia 1.10+ with the shared notebook project instantiated.

Accounts required: None.

# QUBONOTEBOOKS_COLAB_IMPORT_CELL
Base.invokelatest(
    QUBONotebooksBootstrap.warm_notebook_packages!,
    "8-OrderPartitioning",
);

Order-partitioning model¶

Suppose a desk wants to assign indivisible orders to two labeled groups for an educational A/B comparison. Bit $x_j=0$ places order $j$ in group A and $x_j=1$ places it in group B. The value term minimizes the signed group-value difference, while each row of $p$ represents one factor exposure to balance.

The fixture below is independently constructed for this notebook. Order values are in units of USD 100,000; the two risk rows are dimensionless synthetic factor-exposure scores. They are not observed trades, an empirical trading result, or investment advice.

order_names = ["ORD-1", "ORD-2", "ORD-3", "ORD-4", "ORD-5", "ORD-6"]
order_values = [2, 3, 4, 5, 6, 8]  # units of USD 100,000
risk_factor_names = ["market beta proxy", "liquidity proxy"]
risk_exposures = [
    1   1  2  5   4  5
    4  -2  3  6  -3  0
]
value_weight = 2
risk_weight = 1

n_orders = length(order_values)
n_risk_factors = size(risk_exposures, 1)
total_value = sum(order_values)

@assert size(risk_exposures, 2) == n_orders
@assert length(risk_factor_names) == n_risk_factors
println("Synthetic fixture: $n_orders orders, total value = \$$(total_value / 10)M, risk factors = $risk_factor_names")

Synthetic fixture: 6 orders, total value = $2.8M, risk factors = ["market beta proxy", "liquidity proxy"]

For explicit positive weights $a$ and $b$ , we minimize

Q(x)=a\left(T-2\sum_j q_jx_j\right)^2 +b\sum_i\left(\sum_j p_{ij}(2x_j-1)\right)^2.

(1)

The parentheses are essential. For factor $i$ , the inner sum is group B exposure minus group A exposure, so its square penalizes imbalance. Moving the square onto each $(2x_j-1)$ makes that term constant because $(2x_j-1)^2=1$ for binary $x_j$ . The direct functions below are separate from the JuMP expression and serve as the independent reference calculation.

function all_binary_states(n::Integer)
    states = Vector{Vector{Int}}()
    for state in Iterators.product(ntuple(_ -> (0, 1), n)...)
        push!(states, collect(state))
    end
    return states
end

complement_bits(bits) = 1 .- bits

signed_value_imbalance(values, bits) =
    sum(values[j] * (1 - 2 * bits[j]) for j in eachindex(values))

function risk_imbalances(exposures, bits)
    return [
        sum(exposures[i, j] * (2 * bits[j] - 1) for j in axes(exposures, 2))
        for i in axes(exposures, 1)
    ]
end

value_term(values, bits) = signed_value_imbalance(values, bits)^2
risk_term(exposures, bits) = sum(abs2, risk_imbalances(exposures, bits))

function order_partitioning_energy(values, exposures, bits; a, b)
    return a * value_term(values, bits) + b * risk_term(exposures, bits)
end

function incorrect_risk_term_square_inside(exposures, bits)
    return sum(
        exposures[i, j] * (2 * bits[j] - 1)^2
        for i in axes(exposures, 1), j in axes(exposures, 2)
    )
end

function decoded_metrics(values, exposures, bits)
    group_a = findall(==(0), bits)
    group_b = findall(==(1), bits)
    group_a_value = sum(values[group_a])
    group_b_value = sum(values[group_b])
    group_a_risk = [sum(exposures[i, group_a]) for i in axes(exposures, 1)]
    group_b_risk = [sum(exposures[i, group_b]) for i in axes(exposures, 1)]
    signed_risk_imbalance = group_b_risk .- group_a_risk

    return (
        group_a = group_a,
        group_b = group_b,
        group_a_value = group_a_value,
        group_b_value = group_b_value,
        signed_value_imbalance = group_a_value - group_b_value,
        absolute_value_imbalance = abs(group_a_value - group_b_value),
        group_a_risk = group_a_risk,
        group_b_risk = group_b_risk,
        signed_risk_imbalance = signed_risk_imbalance,
        absolute_risk_imbalance = abs.(signed_risk_imbalance),
        aggregate_risk_imbalance = sqrt(sum(abs2, signed_risk_imbalance)),
    )
end

decoded_metrics (generic function with 1 method)

function evaluate_jump_objective(model::Model, variables, bits)
    assignment = Dict(variables[i] => Float64(bits[i]) for i in eachindex(variables))
    return JuMP.value(variable -> assignment[variable], JuMP.objective_function(model))
end

function exact_sampler_optima!(model::Model, variables; atol = 1e-9)
    set_optimizer(model, QUBO.ExactSampler.Optimizer)
    optimize!(model)
    @assert termination_status(model) in (JuMP.MOI.OPTIMAL, JuMP.MOI.LOCALLY_SOLVED)

    energies = [objective_value(model; result = result) for result in 1:result_count(model)]
    best_energy = minimum(energies)
    optima = [
        round.(Int, value.(variables; result = result))
        for result in eachindex(energies)
        if isapprox(energies[result], best_energy; atol = atol, rtol = 0)
    ]
    return best_energy, optima
end

same_states(left, right) = Set(Tuple.(left)) == Set(Tuple.(right))

same_states (generic function with 1 method)

order_model = Model()
@variable(order_model, order_x[1:n_orders], Bin)
@objective(
    order_model,
    Min,
    value_weight * (total_value - 2 * sum(order_values[j] * order_x[j] for j in 1:n_orders))^2 +
    risk_weight * sum(
        (sum(risk_exposures[i, j] * (2 * order_x[j] - 1) for j in 1:n_orders))^2
        for i in 1:n_risk_factors
    ),
)

Exact verification¶

Six binary orders give only $2^6=64$ assignments. We compare the expanded JuMP objective with the direct displayed formula for every one, verify that the correct risk term is nonconstant, and check complement symmetry. The deliberately named incorrect helper exists only as a regression witness for the square-inside-the-sum defect.

states = all_binary_states(n_orders)
@assert length(states) == 64

jump_energies = [evaluate_jump_objective(order_model, order_x, bits) for bits in states]
direct_energies = [
    order_partitioning_energy(
        order_values,
        risk_exposures,
        bits;
        a = value_weight,
        b = risk_weight,
    )
    for bits in states
]
@assert all(isapprox.(jump_energies, direct_energies; atol = 1e-9, rtol = 0))

correct_risk_terms = [risk_term(risk_exposures, bits) for bits in states]
incorrect_risk_terms = [incorrect_risk_term_square_inside(risk_exposures, bits) for bits in states]
@assert length(unique(correct_risk_terms)) > 1
@assert length(unique(incorrect_risk_terms)) == 1
@assert all(
    direct_energies[index] == order_partitioning_energy(
        order_values,
        risk_exposures,
        complement_bits(bits);
        a = value_weight,
        b = risk_weight,
    )
    for (index, bits) in enumerate(states)
)

println("All 64 assignments match the direct grouped-square formula.")
println("Correct risk term values: $(sort(unique(correct_risk_terms))); square-inside regression witness: $(only(unique(incorrect_risk_terms))) for every assignment.")

All 64 assignments match the direct grouped-square formula.
Correct risk term values: [4, 16, 20, 32, 36, 40, 80, 100, 104, 116, 128, 136, 160, 180, 196, 200, 212, 256, 296, 328, 388, 400]; square-inside regression witness: 26 for every assignment.

enumerated_best_energy = minimum(direct_energies)
enumerated_optima = states[findall(==(enumerated_best_energy), direct_energies)]

exact_best_energy, exact_optima = exact_sampler_optima!(order_model, order_x)
@assert enumerated_best_energy == 28
@assert isapprox(exact_best_energy, enumerated_best_energy; atol = 1e-9, rtol = 0)
@assert same_states(exact_optima, enumerated_optima)
@assert length(enumerated_optima) == 2
@assert any(Tuple(bits) in Set(Tuple.(enumerated_optima)) for bits in exact_optima)

println("Enumerated optimum energy: $(round(Int, enumerated_best_energy))")
println("ExactSampler optimum energy: $(round(Int, exact_best_energy))")
println("Global minima: $enumerated_optima")

Enumerated optimum energy: 28
ExactSampler optimum energy: 28

Global minima: [[1, 1, 1, 0, 1, 0], [0, 0, 0, 1, 0, 1]]

Decoded balances¶

We select the optimum whose first order is in group A. Values are converted from the fixture’s USD 100,000 units to millions of dollars. Risk exposure is reported factor by factor; the aggregate risk imbalance is the Euclidean norm of the signed factor differences. The raw QUBO energy remains the weighted sum of squared imbalances.

representative = first(filter(bits -> bits[1] == 0, enumerated_optima))
metrics = decoded_metrics(order_values, risk_exposures, representative)

@assert metrics.group_a == [1, 2, 3, 5]
@assert metrics.group_b == [4, 6]
@assert metrics.group_a_value == sum(order_values[metrics.group_a]) == 15
@assert metrics.group_b_value == sum(order_values[metrics.group_b]) == 13
@assert metrics.signed_value_imbalance == 2
@assert metrics.absolute_value_imbalance == 2
@assert metrics.group_a_risk == [8, 2]
@assert metrics.group_b_risk == [10, 6]
@assert metrics.signed_risk_imbalance == [2, 4]
@assert metrics.absolute_risk_imbalance == [2, 4]
@assert isapprox(metrics.aggregate_risk_imbalance, sqrt(20); atol = 1e-12)
@assert risk_term(risk_exposures, representative) == sum(abs2, metrics.signed_risk_imbalance)

println("x = $representative")
println("Group A orders: $(join(order_names[metrics.group_a], ", "))")
println("Group B orders: $(join(order_names[metrics.group_b], ", "))")
@printf("Group A total: \$%.1fM\n", metrics.group_a_value / 10)
@printf("Group B total: \$%.1fM\n", metrics.group_b_value / 10)
@printf("Signed value imbalance (A-B): \$%+.1fM; absolute: \$%.1fM\n", metrics.signed_value_imbalance / 10, metrics.absolute_value_imbalance / 10)
for i in eachindex(risk_factor_names)
    @printf(
        "%s: A=%d, B=%d, signed imbalance (B-A)=%+d, absolute=%d\n",
        risk_factor_names[i],
        metrics.group_a_risk[i],
        metrics.group_b_risk[i],
        metrics.signed_risk_imbalance[i],
        metrics.absolute_risk_imbalance[i],
    )
end
@printf("Aggregate risk imbalance (L2): %.6f\n", metrics.aggregate_risk_imbalance)
println("Raw QUBO energy: $(order_partitioning_energy(order_values, risk_exposures, representative; a = value_weight, b = risk_weight))")

x = [0, 0, 0, 1, 0, 1]
Group A orders: ORD-1, ORD-2, ORD-3, ORD-5
Group B orders: ORD-4, ORD-6
Group A total: $1.5M
Group B total: $1.3M

Signed value imbalance (A-B): $+0.2M; absolute: $0.2M
market beta proxy: A=8, B=10, signed imbalance (B-A)=+2, absolute=2

liquidity proxy: A=2, B=6, signed imbalance (B-A)=+4, absolute=4
Aggregate risk imbalance (L2): 4.472136
Raw QUBO energy: 28

Complement symmetry¶

Replacing $x$ by $1-x$ swaps the two group labels. Every signed imbalance changes sign, while every squared term and the raw energy remain unchanged.

complement = complement_bits(representative)
complement_metrics = decoded_metrics(order_values, risk_exposures, complement)
complement_energy = order_partitioning_energy(
    order_values,
    risk_exposures,
    complement;
    a = value_weight,
    b = risk_weight,
)

@assert complement in enumerated_optima
@assert complement_metrics.group_a == metrics.group_b
@assert complement_metrics.group_b == metrics.group_a
@assert complement_metrics.group_a_value == metrics.group_b_value
@assert complement_metrics.group_b_value == metrics.group_a_value
@assert complement_metrics.group_a_risk == metrics.group_b_risk
@assert complement_metrics.group_b_risk == metrics.group_a_risk
@assert complement_metrics.signed_value_imbalance == -metrics.signed_value_imbalance
@assert complement_metrics.signed_risk_imbalance == -metrics.signed_risk_imbalance
@assert complement_energy == enumerated_best_energy

println("Complement x = $complement swaps A and B.")
println("Complement energy: $(round(Int, complement_energy))")

Complement x = [1, 1, 1, 0, 1, 0] swaps A and B.
Complement energy: 28

Weight sensitivity¶

Weights $a$ and $b$ express a modeling priority; they do not reveal one universally best partition. Increasing $a$ favors value balance, while increasing $b$ favors factor-risk balance. We enumerate the same fixture under three choices so that the trade-off is transparent rather than inferred from one solver run.

function best_assignment_for_weights(states, values, exposures; a, b)
    energies = [
        order_partitioning_energy(values, exposures, bits; a = a, b = b)
        for bits in states
    ]
    best_energy = minimum(energies)
    optima = states[findall(==(best_energy), energies)]
    representative = first(filter(bits -> bits[1] == 0, optima))
    return best_energy, representative, optima
end

weight_cases = [
    (label = "balanced", a = 2, b = 1),
    (label = "value priority", a = 8, b = 1),
    (label = "risk priority", a = 1, b = 8),
]
weight_results = [
    begin
        energy, bits, optima = best_assignment_for_weights(
            states,
            order_values,
            risk_exposures;
            a = case.a,
            b = case.b,
        )
        (
            label = case.label,
            a = case.a,
            b = case.b,
            energy = energy,
            bits = bits,
            value_imbalance = abs(signed_value_imbalance(order_values, bits)),
            risk_imbalances = risk_imbalances(risk_exposures, bits),
            risk_score = risk_term(risk_exposures, bits),
            degeneracy = length(optima),
        )
    end
    for case in weight_cases
]

@assert weight_results[1].energy == 28
@assert weight_results[2].energy == 40
@assert weight_results[2].value_imbalance == 0
@assert weight_results[2].risk_score == 40
@assert weight_results[3].energy == 68
@assert weight_results[3].value_imbalance == 6
@assert weight_results[3].risk_score == 4
@assert all(result.degeneracy == 2 for result in weight_results)

println("case            a  b  |value A-B|  risk B-A  risk score  energy  representative x")
for result in weight_results
    @printf(
        "%-15s %2d %2d %12d  %-10s %10d %7d  %s\n",
        result.label,
        result.a,
        result.b,
        result.value_imbalance,
        string(result.risk_imbalances),
        result.risk_score,
        result.energy,
        string(result.bits),
    )
end

case            a  b  |value A-B|  risk B-A  risk score  energy  representative x

balanced         2  1            2  [2, 4]             20      28  [0, 0, 0, 1, 0, 1]
value priority   8  1            0  [2, -6]            40      40  [0, 1, 0, 1, 1, 0]
risk priority    1  8            6  [0, -2]             4      68  [0, 0, 0, 1, 1, 0]

The value-priority case reaches equal group value but accepts a larger factor-risk score. The risk-priority case reduces the factor-risk score from 20 to 4 but accepts a value difference of six fixture units (USD 600,000). The balanced case lies between those outcomes. Appropriate weights depend on units, risk governance, and the purpose of the experiment; this tutorial does not prescribe them.

Practice checkpoints¶

Set $a=b=1$ . Predict the optimal value and risk imbalances, then confirm them by enumerating all states.
Change the liquidity exposure of ORD-6 from 0 to 2. Before solving, identify which assertions must be recomputed and why hard-coded solver output is not a correctness proof.
Replace the grouped risk expression in a scratch model with the square-inside expression. Confirm that its risk contribution is identical for all 64 assignments, then restore the corrected model.

For every change, rerun the all-state comparison before trusting the decoded result.

# EXERCISE: start with a = b = 1 and reuse best_assignment_for_weights.
# Keep the exhaustive JuMP-versus-direct comparison as the acceptance check.
nothing

# SOLUTION (hidden in workshop version):
equal_weight_energy, equal_weight_bits, equal_weight_optima = best_assignment_for_weights(
    states,
    order_values,
    risk_exposures;
    a = 1,
    b = 1,
)
@assert equal_weight_energy == 24
@assert abs(signed_value_imbalance(order_values, equal_weight_bits)) == 2
@assert risk_term(risk_exposures, equal_weight_bits) == 20
@assert length(equal_weight_optima) == 2
println("Equal weights: energy=$equal_weight_energy, x=$equal_weight_bits")

Equal weights: energy=24, x=[0, 0, 0, 1, 0, 1]

# EXERCISE: change the liquidity exposure of ORD-6 from 0 to 2.
# Re-enumerate before making any claim about the new optimum.
nothing

# SOLUTION (hidden in workshop version):
changed_exposures = copy(risk_exposures)
changed_exposures[2, 6] = 2
changed_energy, changed_bits, changed_optima = best_assignment_for_weights(
    states,
    order_values,
    changed_exposures;
    a = value_weight,
    b = risk_weight,
)
@assert changed_energy == 28
@assert length(changed_optima) == 4
@assert all(
    order_partitioning_energy(
        order_values,
        changed_exposures,
        bits;
        a = value_weight,
        b = risk_weight,
    ) == changed_energy
    for bits in changed_optima
)
println("Changed exposure: energy=$changed_energy, degeneracy=$(length(changed_optima)), representative x=$changed_bits")

Changed exposure: energy=28, degeneracy=4, representative x=[0, 0, 1, 1, 1, 0]

# EXERCISE: compare the grouped-square risk term with the square-inside term.
# Count the distinct values each expression takes across all 64 states.
nothing

# SOLUTION (hidden in workshop version):
grouped_values = unique(risk_term(risk_exposures, bits) for bits in states)
square_inside_values = unique(
    incorrect_risk_term_square_inside(risk_exposures, bits) for bits in states
)
@assert length(grouped_values) > 1
@assert length(square_inside_values) == 1
println("Grouped square has $(length(grouped_values)) distinct values; square-inside has $(length(square_inside_values)).")

Grouped square has 22 distinct values; square-inside has 1.

Summary¶

Learning objectives met:

Value balance and each factor-risk balance are squared only after summing signed order contributions.
The independently constructed six-order fixture has a nonconstant risk term, unlike the square-inside regression witness.
The direct formula and expanded JuMP objective agree for every assignment.
ExactSampler and independent enumeration return the same two complementary global minima.
Decoded group totals, factor exposures, signed and absolute imbalances, aggregate risk imbalance, and raw energy agree with direct calculations.
Weight changes expose a real trade-off; no one pair of weights is universally best.

Next steps: Later notebooks may compare local or hardware-oriented solvers, but those comparisons are intentionally outside this modeling case study.

Further reading:

The Five Starter Problems paper motivates order partitioning and derives the grouped-square objective.
The QUBO.jl documentation explains the JuMP modeling and exact-sampling interfaces used here.

References¶

A. R. Mazumder and S. Tayur, Five Starter Problems: Solving Quadratic Unconstrained Binary Optimization Models on Quantum Computers, TutORials in Operations Research (2025), pp. 145–183, Mazumder & Tayur (2025).
Companion materials: https://github.com/arulrhikm/Solving-QUBOs-on-Quantum-Computers.
QUBO.jl public API and ExactSampler: https://github.com/JuliaQUBO/QUBO.jl.

This notebook contains original Julia code, prose, and an independently constructed synthetic fixture. The companion repository is cited as context; no source cells, prose, saved output, figures, or other implementation material were copied from it. The example is educational and is not investment advice or an empirical trading result.

References¶

Mazumder, A. R., & Tayur, S. (2025). Five Starter Problems: Solving Quadratic Unconstrained Binary Optimization Models on Quantum Computers. In Tutorials in Operations Research: Advances in Analytics and Operations Research: Improving Decisions to Secure the Future (pp. 145–183). INFORMS. 10.1287/educ.2025.0288