Fix/add objfun deriv bounds; rename vars (#8)

We do two things: first, we fix an erroneous bound on the n'th order derivatives of the transfer fidelity objective function (which previously looser than intended), concurrently adding a bound on the n'th order derivatives of the mixing uniformity objective function. Second, we rename a few variables across the existing codebase for clarity/brevity (e.g., lipschitz_constant -> lipschitz_const).

Author: Luis-Varona

src/EpsilonOptimization/solvers/lipschitz_branch_and_bound.jl

@@ -13,22 +13,20 @@ struct LipschitzBranchAndBound <: AbstractEpsilonSolver
     epsilon::Real
     target::Union{<:Real,Nothing}
     max_iterations::Union{<:Integer,Nothing}
-    lipschitz_constant::Real
+    lipschitz_const::Real
 
     function LipschitzBranchAndBound(
         epsilon::Real,
-        lipschitz_constant::Real;
+        lipschitz_const::Real;
         target::Union{<:Real,Nothing}=nothing,
         max_iterations::Union{<:Integer,Nothing}=nothing,
     )
-        solver = new(epsilon, target, max_iterations, lipschitz_constant)
+        solver = new(epsilon, target, max_iterations, lipschitz_const)
         validate_solver_params(solver)
 
-        if lipschitz_constant <= 0
+        if lipschitz_const <= 0
             throw(
-                ArgumentError(
-                    "Lipschitz constant must be positive, got $lipschitz_constant"
-                ),
+                ArgumentError("Lipschitz constant must be positive, got $lipschitz_const")
             )
         end
 
@@ -49,12 +47,12 @@ struct LBBHyperrectangle{Tx<:AbstractVector{<:Real},Tf<:Real}
     lower_bound::Tf
 
     function LBBHyperrectangle(
-        lower::Tx, upper::Tx, f::Function, lipschitz_constant::Real
+        lower::Tx, upper::Tx, f::Function, lipschitz_const::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
+        lower_bound = f_center - (lipschitz_const / 2) * diameter
         return new{Tx,typeof(f_center)}(lower, upper, center, f_center, lower_bound)
     end
 end
@@ -67,7 +65,7 @@ function _epsilon_minimize_impl(
     f::Function, lower::Tx, upper::Tx, solver::LipschitzBranchAndBound
 ) where {Tx<:AbstractVector{<:Real}}
     epsilon = solver.epsilon
-    lipschitz_constant = solver.lipschitz_constant
+    lipschitz_const = solver.lipschitz_const
 
     if isnothing(solver.target)
         target = -Inf
@@ -81,7 +79,7 @@ function _epsilon_minimize_impl(
         max_iterations = solver.max_iterations
     end
 
-    rect_init = LBBHyperrectangle(lower, upper, f, solver.lipschitz_constant)
+    rect_init = LBBHyperrectangle(lower, upper, f, solver.lipschitz_const)
     rects_cand = BinaryMinHeap{LBBHyperrectangle{Tx,typeof(rect_init.f_center)}}()
     push!(rects_cand, rect_init)
 
@@ -107,7 +105,7 @@ function _epsilon_minimize_impl(
             minimum = rect.f_center
         end
 
-        children = _lbb_split_hyperrectangle(rect, f, lipschitz_constant)
+        children = _lbb_split_hyperrectangle(rect, f, lipschitz_const)
 
         for child in children
             if child.f_center < minimum
@@ -146,7 +144,7 @@ function _epsilon_minimize_impl(
 end
 
 function _lbb_split_hyperrectangle(
-    rect::LBBHyperrectangle, f::Function, lipschitz_constant::Real
+    rect::LBBHyperrectangle, f::Function, lipschitz_const::Real
 )
     lower = rect.lower
     upper = rect.upper
@@ -157,12 +155,12 @@ function _lbb_split_hyperrectangle(
     lower1 = copy(lower)
     upper1 = copy(upper)
     upper1[dim_split] = mid
-    child1 = LBBHyperrectangle(lower1, upper1, f, lipschitz_constant)
+    child1 = LBBHyperrectangle(lower1, upper1, f, lipschitz_const)
 
     lower2 = copy(lower)
     upper2 = copy(upper)
     lower2[dim_split] = mid
-    child2 = LBBHyperrectangle(lower2, upper2, f, lipschitz_constant)
+    child2 = LBBHyperrectangle(lower2, upper2, f, lipschitz_const)
 
     return child1, child2
 end

src/QuantumStateTransfer.jl

@@ -27,7 +27,7 @@ include("core/uniform_mixing.jl")
 include("core/fractional_revival.jl")
 
 # TODO: Exports (add more later)
-export max_state_transfer, check_state_transfer
+export maximize_state_transfer, check_state_transfer
 
 include("startup.jl")
 

src/core/state_transfer.jl

@@ -24,7 +24,7 @@ struct StateTransferMaximizationResult{
 end
 
 function Base.show(io::IO, res::StateTransferMaximizationResult{Tn,Int}) where {Tn}
-    println(io, "Vertex State Transfer Maximization:")
+    println(io, "Results of Vertex State Transfer Maximization")
     println(io, " * Network: $(summary(res.network))")
     println(io, " * Source vertex: $(res.src)")
     println(io, " * Destination vertex: $(res.dst)")
@@ -39,7 +39,7 @@ end
 function Base.show(
     io::IO, res::StateTransferMaximizationResult{Tn,Tuple{Int,Int}}
 ) where {Tn}
-    println(io, "Pair State Transfer Maximization:")
+    println(io, "Results of Pair State Transfer Maximization")
     println(io, " * Network: $(summary(res.network))")
     println(io, " * Source vertices: $(res.src)")
     println(io, " * Destination vertices: $(res.dst)")
@@ -72,7 +72,7 @@ struct StateTransferRecognitionResult{
 end
 
 function Base.show(io::IO, res::StateTransferRecognitionResult{Tn,Int}) where {Tn}
-    println(io, "Vertex State Transfer Recognition:")
+    println(io, "Results of Vertex State Transfer Recognition")
     println(io, " * Network: $(summary(res.network))")
     println(io, " * Source vertex: $(res.src)")
     println(io, " * Destination vertex: $(res.dst)")
@@ -92,7 +92,7 @@ end
 function Base.show(
     io::IO, res::StateTransferRecognitionResult{Tn,Tuple{Int,Int}}
 ) where {Tn}
-    println(io, "Pair State Transfer Recognition:")
+    println(io, "Results of Pair State Transfer Recognition")
     println(io, " * Network: $(summary(res.network))")
     println(io, " * Source vertices: $(res.src)")
     println(io, " * Destination vertices: $(res.dst)")
@@ -110,8 +110,8 @@ function Base.show(
 end
 
 """
-    max_state_transfer(g::AbstractGraph, args...) -> StateTransferMaximizationResult
-    max_state_transfer(A::AbstractMatrix{<:Real}, args...) -> StateTransferMaximizationResult
+    maximize_state_transfer(g::AbstractGraph, args...) -> StateTransferMaximizationResult
+    maximize_state_transfer(A::AbstractMatrix{<:Real}, args...) -> StateTransferMaximizationResult
 
 [TODO: Write here]
 
@@ -133,15 +133,15 @@ end
 # Notes
 [TODO: Refer to [`transfer_fidelity_deriv_bound`](@ref) for proof sketch of bounds]
 """
-function max_state_transfer(g::AbstractGraph, args...)
+function maximize_state_transfer(g::AbstractGraph, args...)
     if !is_simple(g)
         throw(ArgumentError("Graph must be undirected with no self-loops"))
     end
 
-    return max_state_transfer(adjacency_matrix(g), args...)
+    return maximize_state_transfer(adjacency_matrix(g), args...)
 end
 
-function max_state_transfer(
+function maximize_state_transfer(
     A::AbstractMatrix{<:Real},
     src::Tl,
     dst::Tl,
@@ -350,10 +350,8 @@ function _optimize_state_transfer_impl(input::_StateTransferProblemInput)
     end
 
     if method == :lipschitz_bb
-        lipschitz_constant = transfer_fidelity_deriv_bound(A, 1)
-        solver = LipschitzBranchAndBound(
-            epsilon, lipschitz_constant; target=target_infidelity
-        )
+        lipschitz_const = transfer_fidelity_deriv_bound(A, 1)
+        solver = LipschitzBranchAndBound(epsilon, lipschitz_const; target=target_infidelity)
     elseif method == :alpha_bb
         alpha = transfer_fidelity_deriv_bound(A, 2) / 2
         solver = AlphaBranchAndBound(epsilon, alpha; target=target_infidelity)

src/utils.jl

@@ -70,5 +70,23 @@ end
 [TODO: Proof sketch of bound, plus further relevant references?]
 """
 function transfer_fidelity_deriv_bound(A::Matrix{Float64}, order::Int)
-    return maximum(norm.eachcol(A^order))
+    return opnorm(A)^order # Equivalent to `opnorm(A^order)`, since `A` is symmetric
+end
+
+"""
+    mixing_uniformity_deriv_bound(A, order) -> Float64
+
+[TODO: Write here]
+
+# Arguments
+[TODO: Write here]
+
+# Returns
+[TODO: Write here]
+
+# Notes
+[TODO: Proof sketch of bound, plus further relevant references?]
+"""
+function mixing_uniformity_deriv_bound(A::Matrix{Float64}, order::Int)
+    return 2 * opnorm(A)^order # Equivalent to `2 * opnorm(A^order)`, since `A` is symmetric
 end