Add shortcut

Author: sebastiendesignolle

src/bell_frank_wolfe.jl

@@ -39,6 +39,7 @@ function bell_frank_wolfe(
         v0 = one(T),
         epsilon = 10Base.rtoldefault(T),
         verbose = 0,
+        verbose_init = verbose > 0,
         shr2 = NaN,
         mode::Int = 0,
         nb::Int = 10^2,
@@ -54,6 +55,7 @@ function bell_frank_wolfe(
         reset_dots_A::Bool = false, # warm start (technical)
         reset_dots_b::Bool = true, # warm start (technical)
         lazy::Bool = true, # default in FW package is false
+        shortcut = 0, # early termination criterion (outside only)
         max_iteration::Int = 10^9, # default in FW package is 10^4
         nb_increment_interval::Int = 10^4,
         callback_interval::Int = verbose > 0 ? 10^4 : typemax(Int),
@@ -66,7 +68,7 @@ function bell_frank_wolfe(
         kwargs...,
     ) where {T <: Number, N}
     Random.seed!(seed)
-    LMO, DS, m, o, sym, deflate, inflate = _bfw_init(p, v0, prob, marg, o, sym, deflate, inflate, verbose)
+    LMO, DS, m, o, sym, deflate, inflate = _bfw_init(p, v0, prob, marg, o, sym, deflate, inflate, verbose_init)
     if verbose > 0
         println("Visibility: ", v0)
     end
@@ -111,16 +113,14 @@ function bell_frank_wolfe(
             println("Active set initialised")
         end
     end
-    if verbose > 0
-        println()
-    end
     callback = build_callback(
         rp,
         v0,
         ro,
         shr2 ^ (prob ? (N ÷ 2) / 2 : N / 2),
         verbose,
         epsilon,
+        shortcut,
         nb_increment_interval,
         callback_interval,
         hyperplane_interval,
@@ -149,8 +149,7 @@ function bell_frank_wolfe(
     primal = res.primal
     dual_gap = res.dual_gap
     as = res.active_set
-    if verbose ≥ 2
-        println()
+    if verbose == 2
         @printf("Primal: %.2e\n", primal)
         @printf("FW gap: %.2e\n", dual_gap)
         @printf("#Atoms: %d\n", length(as))
@@ -196,14 +195,14 @@ function bell_frank_wolfe(
     if verbose > 0
         if verbose ≥ 2 && mode_last ≥ 0
             @printf("FW gap: %.2e\n", dual_gap) # recomputed FW gap (usually with a more reliable heuristic)
-            println()
         end
         if primal > dual_gap
             @printf("v_c ≤ %f\n", β)
         elseif !isnan(shr2)
             ν = 1 / (1 + norm(vp - as.x, 2))
             @printf("v_c ≥ %f (%f)\n", shr2^(N / 2) * ν * v0, shr2^(N / 2) * v0)
         end
+        println()
     end
     if save
         serialize(file * ".dat", ActiveSetStorage(as))

src/callback.jl

@@ -5,6 +5,7 @@ function build_callback(
         shr2,
         verbose,
         epsilon,
+        shortcut,
         nb_increment_interval,
         callback_interval,
         hyperplane_interval,
@@ -37,34 +38,14 @@ function build_callback(
     verbose_hyperplane = (verbose || save) && hyperplane_interval != typemax(Int)
     verbose_bound = verbose && bound_interval != typemax(Int)
     if verbose
-        @printf(
-            stdout,
-            "%s    %s    %s    %s    %s   %s    %s\n",
-            lpad("Iteration", 12),
-            lpad("Primal", 10),
-            lpad("Dual gap", 10),
-            lpad("Time (sec)", 10),
-            lpad("#It/sec", 10),
-            lpad("#Atoms", 7),
-            lpad("#LMO", 7)
-        )
+        _print_headers()
     end
     function callback(state, active_set, args...)
         if mod(state.t, nb_increment_interval) == 0
             state.lmo.lmo.nb += 1
         end
         if verbose && mod(state.t, callback_interval) == 0
-            @printf(
-                stdout,
-                "%s    %.4e    %.4e    %.4e    %.4e   %s    %s\n",
-                lpad(state.t, 12),
-                state.primal,
-                state.dual_gap,
-                state.time,
-                state.t / state.time,
-                lpad(length(active_set), 7),
-                lpad(state.lmo.lmo.cnt, 7)
-            )
+            _print_callback(state.t, state, active_set)
         end
         if verbose_hyperplane && mod(state.t, hyperplane_interval) == 0
             a = -state.gradient # v*p+(1-v)*o-active_set.x
@@ -83,10 +64,47 @@ function build_callback(
         if save && mod(state.t, save_interval) == 0
             serialize(file * "_tmp.dat", ActiveSetStorage(active_set))
         end
-        # if state.dual_gap < state.primal / 2
-            # return false
-        # end
-        return state.primal > epsilon
+        if shortcut > 0 && state.dual_gap < state.primal / shortcut
+            if verbose
+                _print_callback("Shortcut", state, active_set)
+            end
+            return false
+        end
+        if state.primal ≤ epsilon || state.dual_gap ≤ epsilon
+            if verbose
+                _print_callback("Last", state, active_set)
+            end
+            return false
+        end
+        return true
     end
     return callback
 end
+
+function _print_headers()
+    @printf(
+            stdout,
+            "%s    %s    %s    %s    %s   %s    %s\n",
+            lpad("Iteration", 12),
+            lpad("Primal", 10),
+            lpad("Dual gap", 10),
+            lpad("Time (sec)", 10),
+            lpad("#It/sec", 10),
+            lpad("#Atoms", 7),
+            lpad("#LMO", 7)
+           )
+end
+
+function _print_callback(it, state, active_set)
+    @printf(
+            stdout,
+            "%s    %.4e    %.4e    %.4e    %.4e   %s    %s\n",
+            lpad(it, 12),
+            state.primal,
+            state.dual_gap,
+            state.time,
+            state.t / state.time,
+            lpad(length(active_set), 7),
+            lpad(state.lmo.lmo.cnt, 7)
+           )
+end

src/nonlocality_threshold.jl

@@ -25,17 +25,18 @@ function nonlocality_threshold(
         deflate = identity,
         inflate = identity,
         verbose = 0,
+        shortcut = 100,
         kwargs...,
     ) where {T <: Number, N}
-    _, _, _, o, sym, deflate, inflate = _bfw_init(p, 0, prob, marg, o, sym, deflate, inflate, verbose)
+    _, _, _, o, sym, deflate, inflate = _bfw_init(p, 0, prob, marg, o, sym, deflate, inflate, verbose > 0)
     lower_bound = zero(T)
     upper_bound = one(T)
     local_model = nothing
     bell_inequality = nothing
     while upper_bound - lower_bound > 10.0^(-precision)
-        res = bell_frank_wolfe(p; v0, epsilon, marg, o, sym, deflate, inflate, kwargs...)
+        res = bell_frank_wolfe(p; v0, epsilon, marg, o, sym, deflate, inflate, verbose, verbose_init = false, shortcut, kwargs...)
         x, ds, primal, dual_gap, as, M, β = res
-        if primal > 10epsilon && dual_gap > 10epsilon
+        if dual_gap * shortcut ≥ primal && primal > 10epsilon && dual_gap > 10epsilon
             @warn "Please increase nb or max_iteration"
         end
         if dual_gap < primal

src/utils.jl

@@ -1151,7 +1151,7 @@ function _bfw_init(p::Array{T, N}, v0, prob, marg, o, sym, deflate, inflate, ver
             sym = false
         end
     end
-    if verbose > 1
+    if verbose
         println("   #Inputs: ", all(diff(m) .== 0) ? m[end] - (marg && !prob) : m .- (marg && !prob))
         println(" Symmetric: ", sym)
         println(" Dimension: ", length(deflate(p)))