Finish epsilon-optimization solvers; add skeleton structure

Project.toml

@@ -2,23 +2,25 @@ name = "QuantumStateTransfer"
 uuid = "5afda8b7-781f-4fb2-840f-bdb502253c46"
 keywords = ["quantum computing", "quantum information", "linear algebra", "graph theory", "optimization"]
 license = "MIT"
+version = "0.1.0-dev"
 authors = ["Luis M. B. Varona <lm.varona@outlook.com>", "Nathaniel Johnston <nathaniel.johnston@gmail.com>"]
 description = "A Julia toolbox for modeling state transfer on quantum networks."
 homepage = "https://graphquantum.github.io/QuantumStateTransfer.jl/"
 maintainers = ["Luis M. B. Varona <lm.varona@outlook.com>"]
 readme = "README.md"
 repository = "https://github.com/GraphQuantum/QuantumStateTransfer.jl"
-version = "0.1.0-dev"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 
 [compat]
 DataStructures = "0.18.15 - 0.19"
 Graphs = "1.10"
 LinearAlgebra = "1.10"
+Optim = "1.13.2"
 PrecompileTools = "1.2"
 julia = "1.10"

src/EpsilonOptimization/EpsilonOptimization.jl

@@ -0,0 +1,39 @@
+# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston
+#
+# Licensed under the MIT license <LICENSE or
+# http://opensource.org/licenses/MIT>. This file may not be copied, modified, or
+# distributed except according to those terms.
+
+"""
+    EpsilonOptimization
+
+Optimization algorithms with guaranteed finite ``ϵ``-convergence to the true global minimum
+of ``ℝⁿ → ℝ`` functions over hyperrectangular domains.
+"""
+module EpsilonOptimization
+
+using QuantumStateTransfer: NotImplementedError
+
+using DataStructures: BinaryMinHeap
+using LinearAlgebra: norm
+using Optim
+
+export
+    # Types
+    AbstractEpsilonSolver,
+    EpsilonMinimizationResult,
+
+    # Core functions
+    epsilon_minimize,
+
+    # Solvers
+    LipschitzBranchAndBound,
+    AlphaBranchAndBound
+
+include("types.jl")
+include("core.jl")
+
+include("solvers/lipschitz_branch_and_bound.jl")
+include("solvers/alpha_branch_and_bound.jl")
+
+end

src/EpsilonOptimization/core.jl

@@ -0,0 +1,46 @@
+# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston
+#
+# Licensed under the MIT license <LICENSE or
+# http://opensource.org/licenses/MIT>. This file may not be copied, modified, or
+# distributed except according to those terms.
+
+"""
+    epsilon_minimize(f, lower, upper, solver)
+
+Converge to within an arbitrary ``ϵ`` of the true global minimum of an ``ℝⁿ → ℝ`` function
+`f` over the hyperrectangular domain defined by the bounds `lower` and `upper`, using the
+specified ``ϵ``-optimization `solver`.
+
+# Arguments
+- `f::Function`: The objective function to be minimized. Must map from `ℝⁿ → ℝ`.
+- `lower::Tx<:AbstractVector{<:Real}`: The lower bounds of the hyperrectangular domain.
+- `upper::Tx<:AbstractVector{<:Real}`: The upper bounds of the hyperrectangular domain.
+- `solver::AbstractEpsilonSolver`: The ``ϵ``-optimization solver to use.
+
+# Returns
+- `::EpsilonMinimizationResult{Tx,<:Real}}`: The result of the ``ϵ``-optimization.
+
+# Examples
+[TODO: Write here]
+"""
+function epsilon_minimize(
+    f::Function, lower::Tx, upper::Tx, solver::S
+) where {Tx<:AbstractVector{<:Real},S<:AbstractEpsilonSolver}
+    if length(lower) != length(upper)
+        throw(ArgumentError("Lower and upper bounds must have the same dimension"))
+    end
+
+    if any(lower .> upper)
+        throw(
+            ArgumentError("Lower bound must be entrywise less than or equal to upper bound")
+        )
+    end
+
+    return _epsilon_minimize_impl(f, lower, upper, solver)
+end
+
+function _epsilon_minimize_impl(
+    ::Function, ::Tx, ::Tx, ::S
+) where {Tx<:AbstractVector{<:Real},S<:AbstractEpsilonSolver}
+    throw(NotImplementedError(_epsilon_minimize_impl, :solver, S, AbstractEpsilonSolver))
+end

src/EpsilonOptimization/solvers/alpha_branch_and_bound.jl

@@ -0,0 +1,147 @@
+# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston
+#
+# Licensed under the MIT license <LICENSE or
+# http://opensource.org/licenses/MIT>. This file may not be copied, modified, or
+# distributed except according to those terms.
+
+"""
+    AlphaBranchAndBound <: AbstractEpsilonSolver
+
+[TODO: Write here]
+"""
+struct AlphaBranchAndBound <: AbstractEpsilonSolver
+    epsilon::Real
+    threshold::Union{<:Real,Nothing}
+    max_iterations::Union{<:Integer,Nothing}
+    alpha::Real
+end
+
+"""
+    ABBHyperrectangle{Tx,Tf}
+
+[TODO: Write here]
+"""
+struct ABBHyperrectangle{Tx<:AbstractVector{<:Real},Tf<:Real}
+    lower::Tx
+    upper::Tx
+    x_min::Tx
+    lower_bound::Tf
+
+    function ABBHyperrectangle(
+        lower::Tx, upper::Tx, f::Function, alpha::Real
+    ) where {Tx<:AbstractVector{<:Real}}
+        f_convex(x::Tx) = f(x) + alpha * sum((x .- lower) .* (upper .- x))
+        x0 = lower .+ (upper .- lower) ./ 2
+
+        res = optimize(f_convex, lower, upper, x0, Fminbox(LBFGS()))
+        x_min = res.minimizer
+        lower_bound = res.minimum
+
+        return new{Tx,typeof(lower_bound)}(lower, upper, x_min, lower_bound)
+    end
+end
+
+function Base.isless(rect1::ABBHyperrectangle, rect2::ABBHyperrectangle)
+    return rect1.lower_bound < rect2.lower_bound
+end
+
+function _epsilon_minimize_impl(
+    f::Function, lower::Tx, upper::Tx, solver::AlphaBranchAndBound
+) where {Tx<:AbstractVector{<:Real}}
+    epsilon = solver.epsilon
+    alpha = solver.alpha
+
+    if isnothing(solver.threshold)
+        threshold = -Inf
+    else
+        threshold = solver.threshold
+    end
+
+    if isnothing(solver.max_iterations)
+        max_iterations = Inf
+    else
+        max_iterations = solver.max_iterations
+    end
+
+    rect_init = ABBHyperrectangle(lower, upper, f, alpha)
+    rects_cand = BinaryMinHeap{ABBHyperrectangle{Tx,typeof(rect_init.lower_bound)}}()
+    push!(rects_cand, rect_init)
+
+    minimizer = rect_init.x_min
+    minimum = f(rect_init.x_min)
+    lower_bound = Inf
+    iterations = 0
+
+    while (
+        iterations < max_iterations &&
+        !isempty(rects_cand) &&
+        min(lower_bound, first(rects_cand).lower_bound) <= threshold &&
+        minimum - threshold >= epsilon
+    )
+        rect = pop!(rects_cand)
+        f_val = f(rect.x_min)
+
+        if f_val < minimum
+            minimizer = rect.x_min
+            minimum = f_val
+        end
+
+        children = _abb_split_hyperrectangle(rect, f, alpha)
+
+        for child in children
+            f_val_child = f(child.x_min)
+
+            if f_val_child < minimum
+                minimizer = child.x_min
+                minimum = f_val_child
+            end
+
+            if minimum - child.lower_bound >= epsilon
+                push!(rects_cand, child)
+            else
+                lower_bound = min(lower_bound, child.lower_bound)
+            end
+        end
+
+        iterations += 1
+    end
+
+    if !isempty(rects_cand)
+        lower_bound = min(lower_bound, first(rects_cand).lower_bound)
+    end
+
+    return EpsilonMinimizationResult(
+        solver,
+        epsilon,
+        lower,
+        upper,
+        minimizer,
+        minimum,
+        (
+            epsilon_optimal=minimum - lower_bound < epsilon,
+            threshold_reached=minimum - threshold < epsilon,
+            threshold_unreachable=lower_bound > threshold,
+            iterations=iterations >= max_iterations,
+        ),
+    )
+end
+
+function _abb_split_hyperrectangle(rect::ABBHyperrectangle, f::Function, alpha::Real)
+    lower = rect.lower
+    upper = rect.upper
+
+    dim_split = argmax(upper .- lower)
+    mid = lower[dim_split] + (upper[dim_split] - lower[dim_split]) / 2
+
+    lower1 = copy(lower)
+    upper1 = copy(upper)
+    upper1[dim_split] = mid
+    child1 = ABBHyperrectangle(lower1, upper1, f, alpha)
+
+    lower2 = copy(lower)
+    upper2 = copy(upper)
+    lower2[dim_split] = mid
+    child2 = ABBHyperrectangle(lower2, upper2, f, alpha)
+
+    return child1, child2
+end