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.

Graver Augmentation Multiseed Algorithm (Julia)

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

For local execution from the repository root, run uv sync --group qubo and make setup-julia NOTEBOOK=notebooks_jl/3-GAMA.ipynb before launching Jupyter. This notebook reuses the repo-local Python environment for the D-Wave Ocean stack instead of relying on Julia’s CondaPkg resolver.

This notebook should open directly with the Julia runtime in Google Colab. If it does not, switch the runtime before running the next cell. The setup cell will clone SECQUOIA/QuIP into the Colab workspace when needed, activate notebooks_jl/envs/3-GAMA/Project.toml, install the Python D-Wave Ocean packages into the Colab Python runtime, and print progress for clone, activate, instantiate, and package-load steps. Package precompilation still happens as needed during setup. The first Colab run can still take several minutes.

function load_notebook_bootstrap()
    candidates = (
        joinpath(pwd(), "scripts", "notebook_bootstrap.jl"),
        joinpath(pwd(), "..", "scripts", "notebook_bootstrap.jl"),
        joinpath(pwd(), "QuIP", "scripts", "notebook_bootstrap.jl"),
        joinpath("/content", "QuIP", "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, "QUIP_REPO_DIR", joinpath(pwd(), "QuIP"))
        if !isdir(repo_dir)
            println("[bootstrap] Cloning SECQUOIA/QuIP into $repo_dir")
            run(`git clone --depth 1 https://github.com/SECQUOIA/QuIP.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_notebook_bootstrap()

BOOTSTRAP = QuIPNotebookBootstrap.bootstrap_notebook("3-GAMA")

QUIP_REPO_DIR = BOOTSTRAP.repo_dir
JULIA_NOTEBOOKS_DIR = BOOTSTRAP.notebooks_dir
JULIA_PROJECT_DIR = BOOTSTRAP.project_dir
IN_COLAB = BOOTSTRAP.in_colab

Graver Augmentation Multiseed Algorithm (Julia)

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

if !@isdefined(JuMP)
    using JuMP
end

if !@isdefined(DWave)
    using DWave
end

if !@isdefined(LinearAlgebra)
    using LinearAlgebra
end

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.

if !@isdefined(Measures)
    using Measures
end

if !@isdefined(Random)
    using Random
end

if !@isdefined(Plots)
    using Plots
end

if !@isdefined(StatsBase)
    using StatsBase
end

if !@isdefined(StatsPlots)
    using StatsPlots
end

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
Graver Augmentation Multiseed Algorithm (Python)
D-Wave
Quantum Annealing via D-Wave (Python)