Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Environment Setup

Graver Augmentation Multiseed Algorithm

David E. Bernal Neira
Davidson School of Chemical Engineering, Purdue University

Pedro Maciel Xavier
Davidson School of Chemical Engineering, Purdue University

Azain Khalid
Department of Computer Science, Purdue University
Undergraduate Researcher

Open In Colab SECQUOIA

Environment Setup

If you open this notebook in Google Colab, make sure the runtime is set to Julia before running the next cell. The setup cell will clone SECQUOIA/QuIP into the Colab workspace when needed, activate notebooks_jl/Project.toml, and install the Julia packages used in this notebook.

import Pkg

IN_COLAB = haskey(ENV, "COLAB_RELEASE_TAG") || haskey(ENV, "COLAB_JUPYTER_IP")

function detect_quip_repo_dir()
    candidates = (pwd(), normpath(pwd(), ".."))

    for candidate in candidates
        if isfile(joinpath(candidate, "notebooks_jl", "Project.toml"))
            return candidate
        end
    end

    error("Could not locate the QuIP repository root from $(pwd()).")
end

QUIP_REPO_DIR = if IN_COLAB
    repo_dir = get(ENV, "QUIP_REPO_DIR", joinpath(pwd(), "QuIP"))

    if !isdir(repo_dir)
        run(`git clone --depth 1 https://github.com/SECQUOIA/QuIP.git $repo_dir`)
    end

    repo_dir
else
    detect_quip_repo_dir()
end

JULIA_NOTEBOOKS_DIR = joinpath(QUIP_REPO_DIR, "notebooks_jl")

Pkg.activate(JULIA_NOTEBOOKS_DIR)
Pkg.instantiate(; io = devnull)

cd(JULIA_NOTEBOOKS_DIR)

About this notebook

This notebook makes simple computations of Graver basis. Because of the complexity of these computation, we suggest that for more complicated problems you install the excellent 4ti2 software, an open-source implementation of several routines useful for the study of integer programming through algebraic geometry. It can be used as a stand-alone library and call it from C++ or directly from Julia. In Julia, a binding is provided by lib4ti2_jll.

Introduction to GAMA

The Graver Augmentation Multiseed Algorithm (GAMA) was proposed by two papers by Alghassi, Dridi, and Tayur from the CMU Quantum Computing group. The three main ingredients of this algorithm, designed to solve integer programs with linear constraints and nonlinear objective, are:

  • Computing the Graver basis (or a subset of it) of an integer program.

  • Performing an augmentation.

  • Given that only for certain objective functions, the Graver augmentation is guaranteed to find a globally optimal solution, the algorithm is initialized in several points.

This algorithm can be adapted to take advantage of Quantum Computers by leveraging them as black-box Ising/QUBO problem solvers. In particular, obtaining several feasible solution points for the augmentation and computing the Kernel of the constraint matrix can be posed as QUBO problems. After obtaining these solutions, other routines implemented in classical computers are used to solve the optimization problems, making this a hybrid quantum-classical algorithm.

Introduction to Graver basis computation

A Graver basis is defined as

G(A)=⋃jHj(A)\mathcal{G}(\mathbf{A}) = \bigcup_{j} \mathcal{H}_{j}(\mathbf{A})
(1)

where Hj(A)\mathcal{H}_{j}(\mathbf{A}) are the minimal Hilbert basis of A\mathbf{A} in each orthant.

Equivalently we can define the Graver basis as the ⊑\sqsubseteq-minimal set of a lattice

L(A)={x:Ax=0,x∈Zn}∖{0}=ker⁡A∩Zn\mathcal{L}(\mathbf{A}) = \left\lbrace{}{\mathbf{x} : \mathbf{A} \mathbf{x} = 0, \mathbf{x} \in \mathbb{Z}^{n}}\right\rbrace{} \setminus \left\lbrace{}{0}\right\rbrace{} = \ker A \cap \mathbb{Z}^{n}
(2)

where the partial ordering x⊑y\mathbf{x} \sqsubseteq \mathbf{y} holds whenever xiyi≥0x_i y_i \geq 0 and ∣xi∣≤∣yi∣\left\vert x_i \right\vert \leq \left\vert y_i \right\vert for all ii.

Here we won’t interact with the Quantum Computer. However, we will obtain the Graver basis of a problem using package 4ti2. This notebook studies the behavior of the search algorithm in the case that we only have available a subset of the Graver basis.

Problem statement

We will be solving EXAMPLE 4 in the code, which corresponds to Case 2 in the original GAMA paper. The problem is derived from finance and deals with the maximization of expected returns on investments and the minimization of the variance.

min⁡−∑i=1nμixi+1−εε∑i=1nσi2xi2s.t.Ax=b x∈{−2,−1,0,1,2}n\begin{array}{rll} \displaystyle \min & \displaystyle -\sum_{i = 1}^{n} \mu_{i} x_{i} + \sqrt{\frac{1 - \varepsilon}{\varepsilon} \sum_{i = 1}^{n} \sigma_{i}^2 x_{i}^2 } \\ \textrm{s.t.} & A \mathbf{x} = \mathbf{b} \\ ~ & \mathbf{x} \in \left\lbrace{}{-2, -1, 0, 1, 2}\right\rbrace{}^{n} \end{array}
(3)
f(x) = μ'x + sqrt(((1-ε)/ε) * (σ .^ 2)' * (x .^ 2))
f (generic function with 1 method)

Example

Let

A=[11111111110101010111010101111010100100010010111111010001010110110110010011100000001011101111001000000111110001000100000101010]A = \begin{bmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 1 & 0 & 1 & 0 & 1 & 0 & 1 & 1 & 1 & 0 & 1 & 0 & 1 & 0 \\ 1 & 1 & 1 & 1 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 1 & 1 & 1 & 1 & 1 & 1 \\ 0 & 1 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 1 & 1 & 0 & 1 & 1 & 0 & 1 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 1 & 1 & 0 & 1 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 1 & 0 \\ \end{bmatrix}
(4)

This particular instance of convex INLP has m=5m = 5, n=25n = 25, ε=0.01\varepsilon = 0.01, μi=rand[0,1]\mu_{i} = \textrm{rand}[0, 1], σi=rand[0,μi]\sigma_{i} = \textrm{rand}[0, \mu_{i}]. A∈Bm×nA \in \mathbb{B}^{m \times n} and each bjb_{j} is half the sum of the jj-th row of AA. In this example, b=(9,8,7,5,5)′\mathbf{b} = \left({9, 8, 7, 5, 5}\right)'.

First we would write this problem as a an unconstrained one by penalizing the linear constraints as quadratics in the objective. Let’s first define the problem parameters.

A = [
    1 1 1 1 1 1 1 1 1 1 0 1 0 1 0 1 0 1 1 1 0 1 0 1 0
    1 1 1 1 0 1 0 1 0 0 1 0 0 0 1 0 0 1 0 1 1 1 1 1 1
    0 1 0 0 0 1 0 1 0 1 1 0 1 1 0 1 1 0 0 1 0 0 1 1 1
    0 0 0 0 0 0 0 1 0 1 1 1 0 1 1 1 1 0 0 1 0 0 0 0 0
    0 1 1 1 1 1 0 0 0 1 0 0 0 1 0 0 0 0 0 1 0 1 0 1 0
]

m, n = size(A)

b = vec(ceil.(Int, sum(A; dims = 2) / 2));
x0 = [1, 1, 1, 1, -1, 1, 1, 1, 1, 1, 1, 1, -2, 1, 0, -1, 0, 1, -1, 1, -2, -2, 1, 1, 1]

xl = fill(-2, n)
xu = fill( 2, n);
ε = 0.01
μ = rand(n)       # ~ [0, 1]
σ = rand(n) .* μ; # ~ [0, μ]
using DelimitedFiles
import lib4ti2_jll
import BinaryWrappers
using NPZ

const lib4ti2_bin = BinaryWrappers.@generate_wrappers(lib4ti2_jll)

function has_graver()::Bool
    try
        return success(`$(lib4ti2_bin)/graver --help`)
    catch
        return false
    end
end

function graver()
    run(`$(lib4ti2_bin)/graver --help`)

    return nothing
end

function graver(proj_path::AbstractString; silent::Bool = true)
    if silent
        run(`$(lib4ti2_bin)/graver -q $(proj_path)`)
    else
        run(`$(lib4ti2_bin)/graver $(proj_path)`)
    end

    return nothing
end

function write_mat(path::AbstractString, A)
    m, n = size(A)

    open(path, "w") do io
        println(io, "$m $n")

        join(io, (join(@view(A[i, :]), " ") for i = 1:m), "\n")
    end

    return nothing
end

function read_mat(path::AbstractString, type = Int)
    m, n = parse.(Int, split(readline(path)))

    A = Matrix{type}(undef, m, n)

    open(path, "r") do io
        readline(io)

        for (i, line) in enumerate(eachline(io))
            A[i, :] .= parse.(Int, split(line))
        end
    end

    return A
end

# Compute graver basis locally
function compute_graver_basis_local(A, xl, xu)
    G = nothing

    mktempdir() do path
        proj_path = joinpath(path, "proj")

        write_mat("$(proj_path).mat", A)
        write_mat("$(proj_path).lb", xl)
        write_mat("$(proj_path).ub", xu)

        graver(proj_path)

        G = read_mat("$(proj_path).gra")
    end

    return G    
end

# Fallback: load/download the pre-computed graver basis and read with NPZ.jl  
function download_graver_basis()
    npy_path = joinpath(JULIA_NOTEBOOKS_DIR, "graver.npy")
    G = NPZ.npzread(npy_path) 
    return Array{Int}(G)
end


function graver_basis(A, xl, xu)
    if has_graver()
        return compute_graver_basis_local(A, xl, xu)
    else
        return download_graver_basis()
    end
end
graver_basis (generic function with 1 method)
G = graver_basis(A, xl', xu')
29789×25 Matrix{Int64}: 0 0 0 0 0 0 1 0 0 0 … 0 0 0 0 -1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 -1 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 -1 0 0 0 0 0 0 … 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 1 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 -1 0 0 0 0 -1 0 -2 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 … 0 0 0 0 0 0 -1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 -1 0 -2 0 0 0 0 0 0 0 0 0 0 0 1 0 0 -1 0 0 0 -1 0 0 0 0 0 0 ⋮ ⋮ ⋱ ⋮ ⋮ 2 0 0 0 0 -3 0 0 0 3 -1 -2 0 0 0 0 0 0 1 0 1 2 0 0 0 0 -3 0 0 0 3 -1 -2 0 0 0 0 0 0 2 0 0 2 0 0 0 0 -3 0 0 0 2 -1 -2 0 0 0 0 0 0 2 0 0 2 0 0 0 0 -3 0 0 0 2 … -1 -2 0 0 0 0 0 0 1 0 1 2 0 0 0 0 -3 0 0 0 1 -1 -2 0 0 0 0 0 0 1 0 1 2 0 0 0 0 -3 0 0 0 1 -1 -2 0 0 0 0 0 0 2 0 0 2 0 0 0 0 -3 0 0 0 0 -1 -2 0 0 0 0 0 0 1 0 1 2 0 0 0 0 -3 0 0 0 0 -1 -2 0 0 0 0 0 0 2 0 0 2 0 0 0 0 -3 0 0 0 3 … -1 -2 0 0 0 0 0 0 0 0 2 2 0 0 0 0 -3 0 0 0 2 -1 -2 0 0 0 0 0 0 0 0 2 2 0 0 0 0 -3 0 0 0 1 -1 -2 0 0 0 0 0 0 0 0 2 2 0 0 0 0 -3 0 0 0 0 -1 -2 0 0 0 0 0 0 0 0 2
# Define rules to choose augmentation element, either the best one (argmin) or the first one that is found
function argmin_rule(a)
    local i, y, α = (nothing, NaN, Inf)

    for (_i, (_y, _α)) in enumerate(a)
        if _α < α
            i, y, α = (_i, _y, _α)
        end
    end

    return (i, (y, α))
end

function greedy_rule(a)
    local i, y, α = (nothing, NaN, Inf)

    for (_i, (_y, _α)) in enumerate(a)
        if _α ≉ 0.0
            return (_i, (_y, _α))
        else
            i, y, α = (_i, _y, _α)
        end
    end

    return (i, (y, α))
end
greedy_rule (generic function with 1 method)
# Bisection rules for finding best step size
function bisection_rule(f::Function, g, x, xl = nothing, xu = nothing, laststep = nothing)
    if g == laststep
        return (f(x), 0)
    end

    if xl === nothing
        xl = zeros(length(x))
    end

    if xu === nothing
        xu = fill(2 * maximum(x),length(x))
    end

    u = maximum(xu) - minimum(xl)
    l = -u

    for (i, gi) in enumerate(g)
        if gi >= 1
            u = min(u, floor(Int, (xu[i] - x[i]) / gi))
            l = max(l, ceil(Int, (xl[i] - x[i]) / gi))
        elseif gi <= -1
            u = min(u, floor(Int, (xl[i] - x[i]) / gi))
            l = max(l, ceil(Int, (xu[i] - x[i]) / gi))
        end
    end

    α = u

    while u - l > 1
        α = ifelse(f(x + l*g) < f(x + u*g), l, u)

        p1 = floor(Int, (l+u)/2) - 1
        p2 = floor(Int, (l+u)/2)
        p3 = floor(Int, (l+u)/2) + 1

        if f(x + p1*g) < f(x + p2*g)
            u = floor(Int, (l+u)/2)
        elseif f(x + p3*g) < f(x + p2*g)
            l = floor(Int, (l+u)/2) + 1
        else
            α = p2
            break
        end
    end

    if f(x + l*g) < f(x + u*g) && f(x + l*g) < f(x + α*g)
        α = l
    elseif f(x + u*g) < f(x + α*g)
        α = u
    end

    return (f(x + α*g), α)
end

# We can just have a single step move (works well with greedy approach)
function single_move_rule(f::Function, g, x, xl = nothing, xu = nothing, laststep = nothing)
    if xl === nothing
        xl = zeros(length(x))
    end

    if xu === nothing
        xu = fill(2 * maximum(x),length(x))
    end

    α = 0

    if all(x + g .<= xu) && all(x + g .>= xl)
        if f(x + g) < f(x)
            α = 1
        end
    elseif all(x - g .<= xu) && all(x - g .>= xl)
        if f(x - g) < f(x) && f(x - g) < f(x + g)
            α = -1
        end
    end

    return (f(x + α*g), α)
end
single_move_rule (generic function with 4 methods)
function augmentation(
    f::Function,
    A,
    b,
    x0,
    xl,
    xu;
    silent::Bool = true,
    max_iter::Integer = 1_000,
    step_rule::Function = bisection_rule,
    choice_rule::Function = argmin_rule,
    )
    G = graver_basis(A, xl', xu')

    return augmentation(
        f,
        G,
        A,
        b,
        x0,
        xl,
        xu;
        silent,
        max_iter,
        step_rule,
        choice_rule,
    )
end

function augmentation(
    f::Function,
    G,
    A,
    b,
    x0,
    xl,
    xu;
    silent::Bool = true,
    max_iter::Integer = 1_000,
    step_rule::Function = bisection_rule,
    choice_rule::Function = argmin_rule,
    )
    # Let's perform the augmentation and return the number of steps and the best solution

    # Constraints definition
    isfeasible = (x) -> (A * x == b)

    Δ = 1
    n = length(x0)
    k = 1

    if !silent
        println("Initial point: $(x0)")
        println("Objective function: $(f(x0))")
    end

    
    y = nothing
    x = copy(x0)
    δ = Vector{Int}(undef, n)
    gprev = Vector{Int}(undef, n)

    while Δ != 0 && k < max_iter
        i, (y, Δ) = choice_rule((step_rule(f, g, x, xl, xu, gprev) for g in eachrow(G)))

        δ .= @view(G[i,:]) .* Δ 
        x .+= δ
        gprev .= @view(G[i,:])

        if !silent
            println("Iteration ", k)
            println(i, (y, Δ))
            println("Augmentation direction:", gprev)
            println("Distanced moved:", Δ)
            println("Step taken:", δ)
            println("yective function:", y)
            println(f(x))
            println("Current point:", x)
            println("Are constraints satisfied?", isfeasible(x))
        end

        k += 1
    end

    return (k, y, x)
end
augmentation (generic function with 2 methods)

First, we will prove our augmentation strategies, either best or greedy, and for that last case, either computing the best step or a single move. In the order that was mentioned, the augmentation will take more iterations, but each one of the augmentation steps or iterations is going to be cheaper.

println("Best-augmentation: Choosing among the best step that each element of G can do (via bisection), the one that reduces the most the objective")

iter, f_obj, xf = @time augmentation(
    f, G, A, b, x0, xl, xu;
    step_rule = bisection_rule, choice_rule = argmin_rule
)

println("$(iter), iterations")
println("solution: $(xf)")
println("objective: $(f_obj)")

println("Greedy-best-augmentation: Choosing among the best step that each element of G can do (via bisection), the first one encountered that reduces the objective")

iter,f_obj,xf = @time augmentation(
    f, G, A, b, x0, xl, xu;
    step_rule = bisection_rule, choice_rule = greedy_rule,
)

println("$(iter), iterations")
println("solution: $(xf)")
println("objective: $(f_obj)")

println("Greedy-augmentation: Choosing among the first element of G that with a single step reduces the objective")

iter,f_obj,xf = @time augmentation(
    f, G, A, b, x0, xl, xu;
    step_rule = single_move_rule, choice_rule = greedy_rule,
)

println("$(iter), iterations")
println("solution: $(xf)")
println("objective: $(f_obj)")
Best-augmentation: Choosing among the best step that each element of G can do (via bisection), the one that reduces the most the objective
  6.514303 seconds (74.90 M allocations: 4.122 GiB, 6.52% gc time, 13.41% compilation time)
19, iterations
solution: [-2, 0, 2, 0, 0, 1, 0, 1, 0, 0, -1, 0, -2, 1, 2, 2, 0, 2, 1, 0, -2, 0, 2, 1, 2]
objective: 12.85110679499682
Greedy-best-augmentation: Choosing among the best step that each element of G can do (via bisection), the first one encountered that reduces the objective
  0.600559 seconds (5.38 M allocations: 302.309 MiB, 4.85% gc time, 26.37% compilation time)
36, iterations
solution: [2, 0, 2, 1, 0, 1, 0, 1, 0, 0, 1, 0, -2, 0, 0, 2, 1, -1, 0, 0, -2, 0, 0, 1, 2]
objective: 10.066006763945154
Greedy-augmentation: Choosing among the first element of G that with a single step reduces the objective
  0.870849 seconds (3.07 M allocations: 162.861 MiB, 2.64% gc time, 80.84% compilation time)
33, iterations
solution: [2, 0, 2, 1, 0, 1, 0, 1, 0, 0, 1, 0, -2, 0, 0, 2, 1, -1, 0, 0, -2, 0, 1, 1, 1]
objective: 10.310919445961854

using JuMP, DWave, LinearAlgebra

function get_feasible(A, b; num_reads=1_000)
    m, n = size(A)

    Q = A'A - 2 * diagm(A'b)
    β = b'b

    # Define Binary Quadratic Model using new DWave.jl
    model = Model(DWave.Neal.Optimizer)

    @variable(model, x[1:n], Bin)

    @objective(model, Min, x' * Q * x + β)

    set_optimizer_attribute(model, "num_reads", num_reads)

    optimize!(model)

    X = [
        round.(Int, value.(x; result=i))
        for i = 1:result_count(model)
        if objective_value(model; result=i) ≈ 0.0
    ]

    println("$(length(X)) feasible solutions found.")

    return X
end
get_feasible (generic function with 1 method)
X_feas = get_feasible(A,b; num_reads = 20);
20 feasible solutions found.

We take 20 samples using simulated annealing and notice that most (if not all of them) are feasible and different. Let’s now apply the augmentation procedure to each one of them and record the final objective and the number of iterations it takes. Here we will use the 3rd augmentation strategy (Greedy) because of runtime.

function augmentation_experiment(f, G, A, b, X, xl, xu)
    K = length(X)

    Y = Vector{Float64}(undef, K)
    I = Vector{Int}(undef, K)
    T = Vector{Float64}(undef, K)

    for (i, x) in enumerate(X)
        result = @timed augmentation(
            f, G, A, b, x, xl, xu;
            step_rule   = single_move_rule,
            choice_rule = greedy_rule,
        )

        num_iter, obj_val, _ = result.value

        Y[i] = obj_val
        I[i] = num_iter
        T[i] = result.time
    end

    return (Y, I, T)
end

Y_feas = f.(X_feas)

Y_aug, I_aug, T_aug = augmentation_experiment(f, G, A, b, X_feas, xl, xu);

We record the initial objective function, the one after doing the augmentation, and the number of augmentation steps.

using Measures, Random, Plots, StatsBase, StatsPlots

function plot_augmentation(Y_feas, Y_aug, I_aug)
    plt = plot(;
        plot_title = "Augmentation",
        bottom_margin = 15mm,
    )

    plot!(
        plt, Y_feas;
        label       = "Initial",
        color       = :blue,
        markershape = :square,
    )
    plot!(
        plt, Y_aug;
        label       = "Augmented",
        color       = :orange,
        markershape = :square,
        ylabel      = "Objective Value",
        xlabel      = "Solutions",
        legend      = (0.1, -0.25),
    )

    scatter!(
        twinx(plt), I_aug;
        label       = "Iterations",
        color       = :green,
        legend      = (0.85, -0.25),
        markershape = :diamond,
        ylabel      = "Iterations",
        ylims       = (0, maximum(I_aug) + 1),
    )

    return plt
end

plot_augmentation(Y_feas, Y_aug, I_aug)
Loading...

Notice that we reach the globally optimal solution regardless of the initial point and that even if the initial objective was closer to the optimal objective function, it might take more iterations to reach the optimum.

Now let’s try an extreme case, where we only have 10 of the elements of the Graver basis.

function partial_augmentation_experiment(f, G, A, b, X, xl, xu; seed = nothing, num_samples = 10)
    K = length(X)

    Y = Vector{Float64}(undef, K)
    I = Vector{Int}(undef, K)
    T = Vector{Float64}(undef, K)

    index = sample(1:size(G,1), num_samples, replace = false)
    G_idx = G[index,:]

    for (i, x) in enumerate(X)
        result = @timed augmentation(
            f, G_idx, A, b, x, xl, xu;
            step_rule   = single_move_rule,
            choice_rule = greedy_rule,
        )

        num_iter, obj_val, _ = result.value

        Y[i] = obj_val
        I[i] = num_iter
        T[i] = result.time
    end

    return (Y, I, T)
end

Y_feas = f.(X_feas)

Y_paug, I_paug, T_paug = partial_augmentation_experiment(f, G, A, b, X_feas, xl, xu);
plot_augmentation(Y_feas, Y_paug, I_paug)
Loading...

Here we can barely improve the objective and can only perform a few iterations before we cannot improve the solution. But if we compare the runtimes in both cases, we find that...

function plot_augmentation_runtime(T_aug, T_paug)
    plt = plot(;
        plot_title = "Augmentation Runtime",
        xlabel     = "Solutions",
        ylabel     = "Runtime (s)",
        legend     = (0.8, 0.25),
    )

    plot!(T_aug; label="Complete Basis", markershape=:diamond)
    plot!(T_paug; label="Partial Basis", markershape=:diamond)

    return plt
end

plot_augmentation_runtime(T_aug, T_paug)
Loading...

...the time to do augmentation only having 10 choices is minimal. We can search for a sweet spot in between, with good solutions and little time.

function multiple_partial_augmentation_experiment(f, G, A, b, X, xl, xu; seed = nothing)
    N = 10 # Discretization of the fractions of Graver considered
    M = length(X)
    K = size(G, 1)

    Y = Matrix{Float64}(undef, N, M)
    I = Matrix{Int}(undef, N, M)
    T = Matrix{Float64}(undef, N, M)
    
    for j in 1:N
        Y_p, I_p, T_p = partial_augmentation_experiment(
            f, G, A, b, X, xl, xu; num_samples = trunc(Int, K / N * j)
        )

        Y[j, :] .= Y_p
        I[j, :] .= I_p
        T[j, :] .= T_p
    end

    return (Y, I, T)
end

Y_mpaug, I_mpaug, T_mpaug = multiple_partial_augmentation_experiment(f, G, A, b, X_feas, xl, xu);
function plot_multiple_partial_augmentation(Y_feas, Y_mpaug)
    Y = hcat(Y_feas, Y_mpaug')
    X = permutedims(["Initial";["\$ $(10i) \\%|G| \$" for i = 1:10]])

    plt = boxplot(X, Y;
        plot_title = "Multiple Partial Agumentation Samples",
        size       = (800, 400),
        margin     = 5mm,
        ylabel     = "Objective Value",
        xlabel     = "Sample Sizes",
        legend     = false
    )

    return plt
end

plot_multiple_partial_augmentation(Y_feas, Y_mpaug)
Loading...

References

GAMA
Introduction to GAMA
D-Wave
Quantum Annealing via D-Wave