From 7db5188afd6cd5f1b8738121ebc248b4c6f03ea4 Mon Sep 17 00:00:00 2001 From: Luis-Varona Date: Mon, 27 Oct 2025 17:42:45 -0300 Subject: [PATCH] Finish epsilon-optimization solvers; add skeleton structure This commit does two things: it finishes the epsilon-optimization solvers in the EpsilonOptimization helper submodule, and it adds some skeleton files to the root src folder for inclusion in the QuantumStateTransfer module definition. Next steps are to actually apply the optimization methods to various state transfer objective functions. (Of course, several documentation steps remain 'TODO', but we shall focus on that in later stages of the project.) --- Project.toml | 4 +- .../EpsilonOptimization.jl | 39 +++++ src/EpsilonOptimization/core.jl | 46 ++++++ .../solvers/alpha_branch_and_bound.jl | 147 ++++++++++++++++++ .../solvers/lipschitz_branch_and_bound.jl | 144 +++++++++++++++++ src/EpsilonOptimization/types.jl | 80 ++++++++++ src/QuantumStateTransfer.jl | 22 ++- src/core.jl | 7 + src/startup.jl | 7 + src/types.jl | 86 ++++++++++ src/utils.jl | 7 + 11 files changed, 587 insertions(+), 2 deletions(-) create mode 100644 src/EpsilonOptimization/EpsilonOptimization.jl create mode 100644 src/EpsilonOptimization/core.jl create mode 100644 src/EpsilonOptimization/solvers/alpha_branch_and_bound.jl create mode 100644 src/EpsilonOptimization/solvers/lipschitz_branch_and_bound.jl create mode 100644 src/EpsilonOptimization/types.jl create mode 100644 src/core.jl create mode 100644 src/startup.jl create mode 100644 src/types.jl create mode 100644 src/utils.jl diff --git a/Project.toml b/Project.toml index 58be85a..d9cbf8d 100644 --- a/Project.toml +++ b/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 ", "Nathaniel Johnston "] description = "A Julia toolbox for modeling state transfer on quantum networks." homepage = "https://graphquantum.github.io/QuantumStateTransfer.jl/" maintainers = ["Luis M. B. Varona "] 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" diff --git a/src/EpsilonOptimization/EpsilonOptimization.jl b/src/EpsilonOptimization/EpsilonOptimization.jl new file mode 100644 index 0000000..a880a00 --- /dev/null +++ b/src/EpsilonOptimization/EpsilonOptimization.jl @@ -0,0 +1,39 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . 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 diff --git a/src/EpsilonOptimization/core.jl b/src/EpsilonOptimization/core.jl new file mode 100644 index 0000000..e2285e6 --- /dev/null +++ b/src/EpsilonOptimization/core.jl @@ -0,0 +1,46 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . 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 diff --git a/src/EpsilonOptimization/solvers/alpha_branch_and_bound.jl b/src/EpsilonOptimization/solvers/alpha_branch_and_bound.jl new file mode 100644 index 0000000..c003c32 --- /dev/null +++ b/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 . 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 diff --git a/src/EpsilonOptimization/solvers/lipschitz_branch_and_bound.jl b/src/EpsilonOptimization/solvers/lipschitz_branch_and_bound.jl new file mode 100644 index 0000000..df52346 --- /dev/null +++ b/src/EpsilonOptimization/solvers/lipschitz_branch_and_bound.jl @@ -0,0 +1,144 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . This file may not be copied, modified, or +# distributed except according to those terms. + +""" + LipschitzBranchAndBound <: AbstractEpsilonSolver + +[TODO: Write here] +""" +struct LipschitzBranchAndBound <: AbstractEpsilonSolver + epsilon::Real + threshold::Union{<:Real,Nothing} + max_iterations::Union{<:Integer,Nothing} + lipschitz_constant::Real +end + +""" + LBBHyperrectangle{Tx,Tf} + +[TODO: Write here] +""" +struct LBBHyperrectangle{Tx<:AbstractVector{<:Real},Tf<:Real} + lower::Tx + upper::Tx + center::Tx + f_center::Tf + lower_bound::Tf + + function LBBHyperrectangle( + lower::Tx, upper::Tx, f::Function, lipschitz_constant::Real + ) where {Tx<:AbstractVector{<:Real}} + center = lower .+ (upper .- lower) ./ 2 + f_center = f(center) + diameter = norm(upper .- lower) + lower_bound = f_center - (lipschitz_constant / 2) * diameter + return new{Tx,typeof(f_center)}(lower, upper, center, f_center, lower_bound) + end +end + +function Base.isless(rect1::LBBHyperrectangle, rect2::LBBHyperrectangle) + return rect1.lower_bound < rect2.lower_bound +end + +function _epsilon_minimize_impl( + f::Function, lower::Tx, upper::Tx, solver::LipschitzBranchAndBound +) where {Tx<:AbstractVector{<:Real}} + epsilon = solver.epsilon + lipschitz_constant = solver.lipschitz_constant + + 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 = LBBHyperrectangle(lower, upper, f, solver.lipschitz_constant) + rects_cand = BinaryMinHeap{LBBHyperrectangle{Tx,typeof(rect_init.f_center)}}() + push!(rects_cand, rect_init) + + minimizer = rect_init.center + minimum = rect_init.f_center + 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) + + if rect.f_center < minimum + minimizer = rect.center + minimum = rect.f_center + end + + children = _lbb_split_hyperrectangle(rect, f, lipschitz_constant) + + for child in children + if child.f_center < minimum + minimizer = child.center + minimum = child.f_center + 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 _lbb_split_hyperrectangle( + rect::LBBHyperrectangle, f::Function, lipschitz_constant::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 = LBBHyperrectangle(lower1, upper1, f, lipschitz_constant) + + lower2 = copy(lower) + upper2 = copy(upper) + lower2[dim_split] = mid + child2 = LBBHyperrectangle(lower2, upper2, f, lipschitz_constant) + + return child1, child2 +end diff --git a/src/EpsilonOptimization/types.jl b/src/EpsilonOptimization/types.jl new file mode 100644 index 0000000..5de5ad7 --- /dev/null +++ b/src/EpsilonOptimization/types.jl @@ -0,0 +1,80 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . This file may not be copied, modified, or +# distributed except according to those terms. + +""" + AbstractEpsilonSolver + +Abstract base type for all ``ϵ``-optimization solvers. + +# Interface +Concrete subtypes of `AbstractEpsilonSolver` must implement the following methods: +- `Base.summary(::T) where {T<:AbstractAlgorithm}`: returns a `String` indicating the name + of the algorithm (e.g., `"Lipschitz branch-and-bound"`). + +Concrete subtypes of `AbstractEpsilonSolver` must also have the following fields: +- `epsilon::Real`: The tolerance for ``ϵ``-convergence. +- `threshold::Union{<:Real,Nothing}`: An optional threshold value resulting in termination + once either (1) a function value less than `threshold + epsilon` is found or (2) it is + determined that no function value less than or equal to `threshold` exists. +- `max_iterations::Union{<:Integer,Nothing}`: An optional maximum number of iterations after + which the algorithm will terminate. +""" +abstract type AbstractEpsilonSolver end + +""" + EpsilonMinimizationResult{Tx,Tf} + +Output struct for ``ϵ``-optimization results. + +# Fields +- `algorithm::AbstractEpsilonSolver`: The algorithm used for the optimization. +- `lower::Tx<:AbstractVector{<:Real}}`: The lower bounds of the hyperrectangular domain. +- `upper::Tx<:AbstractVector{<:Real}}`: The upper bounds of the hyperrectangular domain. +- `epsilon::Real`: The tolerance provided for ``ϵ``-convergence. +- `minimizer::Tx<:AbstractVector{<:Real}}`: The estimated location of the global minimizer. +- `minimum::Tf<:Real`: The estimated global minimum function value. +- `stopped_by::NamedTuple`: The reason(s) for termination, namely: + - `:epsilon_optimal::Bool`: Whether ``ϵ``-optimality was achieved. + - `:threshold_reached::Bool`: Whether ta function value less than `threshold + epsilon` + was found for some optionally provided `threshold`. + - `:threshold_unreachable::Bool`: Whether it was determined that no function value less + than or equal to `threshold` existed for some optionally provided `threshold`. + - `:iterations::Bool`: Whether termination occurred due to reaching the maximum number + of iterations for some optionally provided `max_iterations`. +""" +struct EpsilonMinimizationResult{Tx<:AbstractVector{<:Real},Tf<:Real} + algorithm::AbstractEpsilonSolver + lower::Tx + upper::Tx + epsilon::Real + minimizer::Tx + minimum::Tf + stopped_by::NamedTuple{ + (:epsilon_optimal, :threshold_reached, :threshold_unreachable, :iterations), + Tuple{Bool,Bool,Bool,Bool}, + } +end + +# The `Base.show` override here takes heavy inspiration from the Optim.jl package +function Base.show(io::IO, res::EpsilonMinimizationResult) + println(io, "Results of Epsilon Minimization Algorithm") + print(io, " * Status: ") + + if res.stopped_by.epsilon_optimal || res.stopped_by.threshold_reached + println(io, "success") + else + println(io, "failure") + end + + println(io, " * Algorithm: $(summary(res.algorithm))") + println(io, " * Lower bounds: $(res.lower)") + println(io, " * Upper bounds: $(res.upper)") + println(io, " * Epsilon tolerance: $(res.epsilon)") + println(io, " * Estimated minimizer: $(res.minimizer)") + println(io, " * Estimated minimum: $(res.minimum)") + + return nothing +end diff --git a/src/QuantumStateTransfer.jl b/src/QuantumStateTransfer.jl index fb37123..8b36be3 100644 --- a/src/QuantumStateTransfer.jl +++ b/src/QuantumStateTransfer.jl @@ -4,8 +4,28 @@ # http://opensource.org/licenses/MIT>. This file may not be copied, modified, or # distributed except according to those terms. +""" + QuantumStateTransfer + +A Julia toolbox for modeling state transfer on quantum networks. +""" module QuantumStateTransfer -# Write your package code here. +using DataStructures +using Graphs +using LinearAlgebra +using PrecompileTools: @setup_workload, @compile_workload + +include("utils.jl") +include("types.jl") + +include("EpsilonOptimization/EpsilonOptimization.jl") +using .EpsilonOptimization + +include("core.jl") + +# TODO: Exports + +include("startup.jl") end diff --git a/src/core.jl b/src/core.jl new file mode 100644 index 0000000..fa351df --- /dev/null +++ b/src/core.jl @@ -0,0 +1,7 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . This file may not be copied, modified, or +# distributed except according to those terms. + +# TODO: Write here diff --git a/src/startup.jl b/src/startup.jl new file mode 100644 index 0000000..fa351df --- /dev/null +++ b/src/startup.jl @@ -0,0 +1,7 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . This file may not be copied, modified, or +# distributed except according to those terms. + +# TODO: Write here diff --git a/src/types.jl b/src/types.jl new file mode 100644 index 0000000..018e7fc --- /dev/null +++ b/src/types.jl @@ -0,0 +1,86 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . This file may not be copied, modified, or +# distributed except according to those terms. + +""" + NotImplementedError{Nothing}(f, subtype, abstracttype) + NotImplementedError{Symbol}(f, arg, subtype, abstracttype) + +An exception indicating that a function lacks dispatch to handle a specific argument type. + +Semantically, this differs from `MethodError` in that it connotes a developer-side failure +to implement a method rather than erroneous user input. Throughout this package, it is often +used to warn when an existing function with multiple dispatch on some abstract type is +called on a newly created subtype for which no method has been defined. + +# Fields +- `f::Function`: the function called. +- `arg::Symbol`: the name of the argument with the unsupported type, if the function has + multiple arguments. If the function has only one argument, this field should be set to + `nothing`. +- `subtype::Type`: the type of the argument. May be the actual concrete type or some + intermediate supertype. (For instance, if the relevant input has concrete type `A` with + hierarchy `A <: B <: C` and the `abstracttype` field is `C`, then both `A` and `B` are + perfectly valid choices for `subtype`.) +- `abstracttype::Type`: the abstract type under which the argument is meant to fall. + +# Constructors +- `NotImplementedError(::Function, ::Type, ::Type)`: constructs a new `NotImplementedError` + instance for a single-argument function. Throws an error if the second type is not + abstract or the first type is not a subtype of the second. +- `NotImplementedError(::Function, ::Symbol, ::Type, ::Type)`: constructs a new + `NotImplementedError` instance for a multi-argument function. Throws an error if the + second type is not abstract or the first type is not a subtype of the second. + +# Supertype Hierarchy +`NotImplementedError` <: `Exception` + +# Notes +This struct was taken from one of the authors' other packages, MatrixBandwidth.jl. +""" +struct NotImplementedError{T<:Union{Nothing,Symbol}} <: Exception + f::Function + arg::T + subtype::Type + abstracttype::Type + + function NotImplementedError(f::Function, subtype::Type, abstracttype::Type) + return NotImplementedError(f, nothing, subtype, abstracttype) + end + + function NotImplementedError( + f::Function, arg::T, subtype::Type, abstracttype::Type + ) where {T<:Union{Nothing,Symbol}} + if !isabstracttype(abstracttype) + throw(ArgumentError("Expected an abstract type, got $abstracttype")) + end + + if !(subtype <: abstracttype) + throw(ArgumentError("Expected a subtype of $abstracttype, got $subtype")) + end + + return new{T}(f, arg, subtype, abstracttype) + end +end + +function Base.showerror(io::IO, e::NotImplementedError{Nothing}) + print( + io, + """NotImplementedError with $(e.subtype): + $(e.f) is not yet implemented for this subtype of $(e.abstracttype). + Try defining method dispatch manually if this is a newly created subtype.""", + ) + return nothing +end + +function Base.showerror(io::IO, e::NotImplementedError{Symbol}) + print( + io, + """NotImplementedError with argument $(e.arg)::$(e.subtype): + $(e.f) is not yet implemented for this subtype of $(e.abstracttype). + Try defining method dispatch manually if this is a newly created subtype.""", + ) + return nothing +end diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..fa351df --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,7 @@ +# Copyright 2025 Luis M. B. Varona and Nathaniel Johnston +# +# Licensed under the MIT license . This file may not be copied, modified, or +# distributed except according to those terms. + +# TODO: Write here