From ecf8db0b4201fb000724b6df23178fb8dce69ae0 Mon Sep 17 00:00:00 2001 From: dhendryc <92737336+dhendryc@users.noreply.github.com> Date: Wed, 7 Jun 2023 11:18:08 +0200 Subject: [PATCH] Support self concordant functions (#130) Added options for user to pick FW variant. Currently blended pairwise conditional gradient, away frank wolfe, vanilla frank wolfe and blended conditional gradient are supported. The user is also able to add their own variant of Frank-Wolfe. --- Project.toml | 2 +- examples/portfolio.jl | 2 +- src/Boscia.jl | 1 + src/callbacks.jl | 58 +++++----- src/custom_bonobo.jl | 7 ++ src/frank_wolfe_variants.jl | 217 ++++++++++++++++++++++++++++++++++++ src/heuristics.jl | 26 +++-- src/interface.jl | 122 +++++++++++++------- src/node.jl | 74 ++++++++++-- test/interface_test.jl | 105 +++++++++++++++++ test/mean_risk.jl | 2 +- test/runtests.jl | 1 + 12 files changed, 524 insertions(+), 93 deletions(-) create mode 100644 src/frank_wolfe_variants.jl diff --git a/Project.toml b/Project.toml index 48761e9db..d3b4bf897 100644 --- a/Project.toml +++ b/Project.toml @@ -22,7 +22,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Bonobo = "0.1.3" DataStructures = "0.18" Distributions = "0.25" -FrankWolfe = "0.2.11" +FrankWolfe = "0.2.24" HiGHS = "1" MathOptInterface = "1" MathOptSetDistances = "0.2" diff --git a/examples/portfolio.jl b/examples/portfolio.jl index 3e09b587b..f5b57741f 100644 --- a/examples/portfolio.jl +++ b/examples/portfolio.jl @@ -53,7 +53,7 @@ const Mi = (Ai + Ai') / 2 return storage end - x, _, result = Boscia.solve(f, grad!, lmo, verbose=true) + x, _, result = Boscia.solve(f, grad!, lmo, verbose=true, time_limit=600) @test dot(ai, x) <= bi + 1e-6 @test f(x) <= f(result[:raw_solution]) + 1e-6 end diff --git a/src/Boscia.jl b/src/Boscia.jl index 08ba19c26..e73ee7b51 100644 --- a/src/Boscia.jl +++ b/src/Boscia.jl @@ -14,6 +14,7 @@ import MathOptSetDistances as MOD include("time_tracking_lmo.jl") include("bounds.jl") +include("frank_wolfe_variants.jl") include("node.jl") include("custom_bonobo.jl") include("callbacks.jl") diff --git a/src/callbacks.jl b/src/callbacks.jl index 629a81f54..03b1d297b 100644 --- a/src/callbacks.jl +++ b/src/callbacks.jl @@ -26,23 +26,25 @@ function build_FW_callback(tree, min_number_lower, check_rounding_value::Bool, f @assert is_linear_feasible(tree.root.problem.lmo, state.v) end # @assert is_linear_feasible(tree.root.problem.lmo, state.x) - push!(fw_iterations, state.t) - if ncalls != state.lmo.ncalls - ncalls = state.lmo.ncalls - (best_v, best_val) = - find_best_solution(tree.root.problem.f, tree.root.problem.lmo.lmo.o, vars) - if best_val < tree.incumbent - tree.root.updated_incumbent[] = true - node = tree.nodes[tree.root.current_node_id[]] - sol = FrankWolfeSolution(best_val, best_v, node, :SCIP) - push!(tree.solutions, sol) - if tree.incumbent_solution === nothing || - sol.objective < tree.incumbent_solution.objective - tree.incumbent_solution = sol + + if state.lmo !== nothing # can happen with using Blended Conditional Gradient + if ncalls != state.lmo.ncalls + ncalls = state.lmo.ncalls + (best_v, best_val) = + find_best_solution(tree.root.problem.f, tree.root.problem.lmo.lmo.o, vars, tree.root.options[:domain_oracle]) + if best_val < tree.incumbent + tree.root.updated_incumbent[] = true + node = tree.nodes[tree.root.current_node_id[]] + sol = FrankWolfeSolution(best_val, best_v, node, :SCIP) + push!(tree.solutions, sol) + if tree.incumbent_solution === nothing || + sol.objective < tree.incumbent_solution.objective + tree.incumbent_solution = sol + end + tree.incumbent = best_val + Bonobo.bound!(tree, node.id) end - tree.incumbent = best_val - Bonobo.bound!(tree, node.id) end end @@ -50,19 +52,21 @@ function build_FW_callback(tree, min_number_lower, check_rounding_value::Bool, f return false end - val = tree.root.problem.f(state.v) - if val < tree.incumbent - tree.root.updated_incumbent[] = true - #TODO: update solution without adding node - node = tree.nodes[tree.root.current_node_id[]] - sol = FrankWolfeSolution(val, copy(state.v), node, :vertex) - push!(tree.solutions, sol) - if tree.incumbent_solution === nothing || - sol.objective < tree.incumbent_solution.objective - tree.incumbent_solution = sol + if tree.root.options[:domain_oracle](state.v) + val = tree.root.problem.f(state.v) + if val < tree.incumbent + tree.root.updated_incumbent[] = true + #TODO: update solution without adding node + node = tree.nodes[tree.root.current_node_id[]] + sol = FrankWolfeSolution(val, copy(state.v), node, :vertex) + push!(tree.solutions, sol) + if tree.incumbent_solution === nothing || + sol.objective < tree.incumbent_solution.objective + tree.incumbent_solution = sol + end + tree.incumbent = val + Bonobo.bound!(tree, node.id) end - tree.incumbent = val - Bonobo.bound!(tree, node.id) end node = tree.nodes[tree.root.current_node_id[]] diff --git a/src/custom_bonobo.jl b/src/custom_bonobo.jl index 837378bd8..0fcf91775 100644 --- a/src/custom_bonobo.jl +++ b/src/custom_bonobo.jl @@ -106,5 +106,12 @@ function Bonobo.get_solution( tree::Bonobo.BnBTree{N,R,V,S}; result=1, ) where {N,R,V,S<:FrankWolfeSolution{N,V}} + if isempty(tree.solutions) + @warn "There is no solution in the tree. This behaviour can happen if you have supplied + \na custom domain oracle. In that case, try to increase the time limit. If you have not specified a + \ndomain oracle, please report!" + @assert tree.root.problem.solving_stage == TIME_LIMIT_REACHED + return nothing + end return tree.solutions[result].solution end diff --git a/src/frank_wolfe_variants.jl b/src/frank_wolfe_variants.jl new file mode 100644 index 000000000..0c18a6318 --- /dev/null +++ b/src/frank_wolfe_variants.jl @@ -0,0 +1,217 @@ + +""" +Frank-Wolfe variant used to compute the problems at node level. +A `FrankWolfeVariant` must implement +``` +solve_frank_wolfe(fw::FrankWolfeVariant, f, grad!, lmo, active_set, line_search, epsilon, max_iteration, + added_dropped_vertices, use_extra_vertex_storage, callback, lazy, timeout, verbose, workspace)) +``` +It may also implement `build_frank_wolfe_workspace(x)` which creates a +workspace structure that is passed as last argument to `solve_frank_wolfe`. +""" +abstract type FrankWolfeVariant end + +# default printing for FrankWolfeVariant is just showing the type +Base.print(io::IO, fw::FrankWolfeVariant) = print(io, split(string(typeof(fw)), ".")[end]) + +""" + solve_frank_wolfe(fw::FrankWolfeVariant, f, grad!, lmo, active_set, line_search, epsilon, max_iteration, + added_dropped_vertices, use_extra_vertex_storage, callback, lazy, timeout, verbose, workspace) + +Returns the optimal solution x to the node problem, its primal and dual gap and the active set. +""" +function solve_frank_wolfe end + +build_frank_wolfe_workspace(::FrankWolfeVariant, x) = nothing + + +""" + Away-Frank-Wolfe + +In every iteration, it computes the worst performing vertex, called away vertex, in the active set with regard to the gradient. +If enough local progress can be made, weight is shifted from the away vertex to all other vertices. + +In case lazification is activated, the FW vertex is only computed if not enough local progress can be guaranteed. +""" +struct AwayFrankWolfe <: FrankWolfeVariant end + +function solve_frank_wolfe( + frank_wolfe_variant::AwayFrankWolfe, + f, + grad!, + lmo, + active_set; + line_search::FrankWolfe.LineSearchMethod=FrankWolfe.Adaptive(), + epsilon=1e-7, + max_iteration=10000, + add_dropped_vertices=false, + use_extra_vertex_storage=false, + extra_vertex_storage=nothing, + callback=nothing, + lazy=false, + timeout=Inf, + verbose=false, + workspace=nothing +) + x, _, primal, dual_gap, _, active_set = FrankWolfe.away_frank_wolfe( + f, + grad!, + lmo, + active_set, + epsilon=epsilon, + max_iteration=max_iteration, + line_search=line_search, + callback=callback, + lazy=lazy, + timeout=timeout, + add_dropped_vertices=add_dropped_vertices, + use_extra_vertex_storage=use_extra_vertex_storage, + extra_vertex_storage=extra_vertex_storage, + verbose=verbose, + ) + + return x, primal, dual_gap, active_set +end + +Base.print(io::IO, ::AwayFrankWolfe) = print(io, "Away-Frank-Wolfe") + + +""" + Blended Conditional Gradient +""" +struct Blended <: FrankWolfeVariant end + +function solve_frank_wolfe( + frank_wolfe_variant::Blended, + f, + grad!, + lmo, + active_set; + line_search::FrankWolfe.LineSearchMethod=FrankWolfe.Adaptive(), + epsilon=1e-7, + max_iteration=10000, + add_dropped_vertices=false, + use_extra_vertex_storage=false, + extra_vertex_storage=nothing, + callback=nothing, + lazy=false, + timeout=Inf, + verbose=false, + workspace=nothing +) + x,_, primal, dual_gap,_, active_set = blended_conditional_gradient( + f, + grad!, + lmo, + active_set, + line_search=line_search, + epsilon=epsilon, + max_iteration=max_iteration, + add_dropped_vertices=add_dropped_vertices, + use_extra_vertex_storage=use_extra_vertex_storage, + extra_vertex_storage=extra_vertex_storage, + callback=callback, + timeout=timeout, + verbose=verbose, + ) + + return x, primal, dual_gap, active_set +end + +Base.print(io::IO, ::Blended) = print(io, "Blended Conditional Gradient") + +""" + Blended Pairwise Conditional Gradient +""" +struct BPCG <: FrankWolfeVariant end + +function solve_frank_wolfe( + frank_wolfe_variant::BPCG, + f, + grad!, + lmo, + active_set; + line_search::FrankWolfe.LineSearchMethod=FrankWolfe.Adaptive(), + epsilon=1e-7, + max_iteration=10000, + add_dropped_vertices=false, + use_extra_vertex_storage=false, + extra_vertex_storage=nothing, + callback=nothing, + lazy=false, + timeout=Inf, + verbose=false, + workspace=nothing +) + x, _, primal, dual_gap, _, active_set = FrankWolfe.blended_pairwise_conditional_gradient( + f, + grad!, + lmo, + active_set, + line_search=line_search, + epsilon=epsilon, + max_iteration=max_iteration, + add_dropped_vertices=add_dropped_vertices, + use_extra_vertex_storage=use_extra_vertex_storage, + extra_vertex_storage=extra_vertex_storage, + callback=callback, + lazy=lazy, + timeout=timeout, + verbose=verbose + ) + + return x, primal, dual_gap, active_set +end + +Base.print(io::IO, ::BPCG) = print(io, "Blended Pairwise Conditional Gradient") + +""" + Vanilla-Frank-Wolfe + +The standard variant of Frank-Wolfe. In each iteration, the vertex v minimizing ∇f * (x-v) is computed. + +Lazification cannot be used in this setting. +""" +struct VanillaFrankWolfe <: FrankWolfeVariant end + +function solve_frank_wolfe( + frank_wolfe_variant::VanillaFrankWolfe, + f, + grad!, + lmo, + active_set; + line_search::FrankWolfe.LineSearchMethod=FrankWolfe.Adaptive(), + epsilon=1e-7, + max_iteration=10000, + add_dropped_vertices=false, + use_extra_vertex_storage=false, + extra_vertex_storage=nothing, + callback=nothing, + lazy=false, + timeout=Inf, + verbose=false, + workspace=nothing +) + # If the flag away_steps is set to false, away_frank_wolfe performs Vanilla. + # Observe that the lazy flag is only observed if away_steps is set to true, so it can neglected. + x, _, primal, dual_gap, _, active_set = FrankWolfe.away_frank_wolfe( + f, + grad!, + lmo, + active_set, + away_steps=false, + epsilon=epsilon, + max_iteration=max_iteration, + line_search=line_search, + callback=callback, + timeout=timeout, + add_dropped_vertices=add_dropped_vertices, + use_extra_vertex_storage=use_extra_vertex_storage, + extra_vertex_storage=extra_vertex_storage, + verbose=verbose, + ) + + return x, primal, dual_gap, active_set +end + +Base.print(io::IO, ::VanillaFrankWolfe) = print(io, "Vanilla-Frank-Wolfe") diff --git a/src/heuristics.jl b/src/heuristics.jl index f9e412b34..ebbcd47a4 100644 --- a/src/heuristics.jl +++ b/src/heuristics.jl @@ -3,34 +3,38 @@ Finds the best solution in the SCIP solution storage, based on the objective function `f`. Returns the solution vector and the corresponding best value. """ -function find_best_solution(f::Function, o::SCIP.Optimizer, vars::Vector{MOI.VariableIndex}) +function find_best_solution(f::Function, o::SCIP.Optimizer, vars::Vector{MOI.VariableIndex}, domain_oracle) sols_vec = unsafe_wrap(Vector{Ptr{Cvoid}}, SCIP.LibSCIP.SCIPgetSols(o), SCIP.LibSCIP.SCIPgetNSols(o)) best_val = Inf best_v = nothing for sol in sols_vec v = SCIP.sol_values(o, vars, sol) - val = f(v) - if val < best_val - best_val = val - best_v = v + if domain_oracle(v) + val = f(v) + if val < best_val + best_val = val + best_v = v + end end end - @assert isfinite(best_val) + #@assert isfinite(best_val) -> not necessarily the case if the domain oracle is not the default. return (best_v, best_val) end -function find_best_solution(f::Function, o::MOI.AbstractOptimizer, vars::Vector{MOI.VariableIndex}) +function find_best_solution(f::Function, o::MOI.AbstractOptimizer, vars::Vector{MOI.VariableIndex}, domain_oracle) nsols = MOI.get(o, MOI.ResultCount()) @assert nsols > 0 best_val = Inf best_v = nothing for sol_idx in 1:nsols xv = [MOI.get(o, MOI.VariablePrimal(sol_idx), xi) for xi in vars] - val = f(xv) - if val < best_val - best_val = val - best_v = xv + if domain_oracle(xv) + val = f(xv) + if val < best_val + best_val = val + best_v = xv + end end end return (best_v, best_val) diff --git a/src/interface.jl b/src/interface.jl index a30d1a39e..07bab846e 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -1,34 +1,44 @@ """ solve -f - objective function oracle. -g - oracle for the gradient of the objective. -lmo - a MIP solver instance (SCIP) encoding the feasible region. -traverse_strategy - encodes how to choose the next node for evaluation. +f - objective function oracle. +g - oracle for the gradient of the objective. +lmo - a MIP solver instance (SCIP) encoding the feasible region. +traverse_strategy - encodes how to choose the next node for evaluation. By default the node with the best lower bound is picked. -branching_strategy - by default we branch on the entry which is the farthest +branching_strategy - by default we branch on the entry which is the farthest away from being an integer. -fw_epsilon - the tolerance for FrankWolfe in the root node. -verbose - if true, a log and solution statistics are printed. -dual_gap - if this absolute dual gap is reached, the algorithm stops. -rel_dual_gap - if this relative dual gap is reached, the algorithm stops. -time_limit - algorithm will stop if the time limit is reached. Depending on the problem +variant - variant of FrankWolfe to be used to solve the node problem. + Options: FW -- Vanilla FrankWolfe + AFW -- Away FrankWolfe + BPCG -- Blended Pairwise Conditional Gradient +line_search - specifies the Line Search method used in the FrankWolfe variant. + Default is the Adaptive Line Search. For other types, check the FrankWolfe.jl package. +fw_epsilon - the tolerance for FrankWolfe in the root node. +verbose - if true, a log and solution statistics are printed. +dual_gap - if this absolute dual gap is reached, the algorithm stops. +rel_dual_gap - if this relative dual gap is reached, the algorithm stops. +time_limit - algorithm will stop if the time limit is reached. Depending on the problem it is possible that no feasible solution has been found yet. -print_iter - encodes after how many proccessed nodes the current node and solution status +print_iter - encodes after how many proccessed nodes the current node and solution status is printed. Will always print if a new integral solution has been found. -dual_gap_decay_factor - the FrankWolfe tolerance at a given level i in the tree is given by +dual_gap_decay_factor - the FrankWolfe tolerance at a given level i in the tree is given by fw_epsilon * dual_gap_decay_factor^i until we reach the min_node_fw_epsilon. -max_fw_iter - maximum number of iterations in a FrankWolfe run. -min_number_lower - If not Inf, evaluation of a node is stopped if at least min_number_lower nodes have a better +max_fw_iter - maximum number of iterations in a FrankWolfe run. +min_number_lower - If not Inf, evaluation of a node is stopped if at least min_number_lower nodes have a better lower bound. -min_node_fw_epsilon - smallest fw epsilon possible, see dual_gap_decay_factor. -min_fw_iterations - the minimum number of FrankWolfe iterations in the node evaluation. -max_iteration_post - maximum number of iterations in a FrankWolfe run during postsolve -dual_tightening - whether to use dual tightening techniques (make sure your function is convex!) +min_node_fw_epsilon - smallest fw epsilon possible, see dual_gap_decay_factor. +min_fw_iterations - the minimum number of FrankWolfe iterations in the node evaluation. +max_iteration_post - maximum number of iterations in a FrankWolfe run during postsolve +dual_tightening - whether to use dual tightening techniques (make sure your function is convex!) global_dual_tightening - dual tightening maintained globally valid (when new solutions are found) -bnb_callback - an optional callback called at every node of the tree, for example for heuristics -strong_convexity - strong convexity of the function, used for tighter dual bound at every node -start_solution - initial solution to start with an incumbent +bnb_callback - an optional callback called at every node of the tree, for example for heuristics +strong_convexity - strong convexity of the function, used for tighter dual bound at every node +domain_oracle - For a point x: returns true if x is in the domain of f, else false. Per default is true. + In case of the non trivial domain oracle, the starting point has to be feasible for f. Also, depending + on the Line Search method, you might have to provide the domain oracle to it, too. +start_solution - initial solution to start with an incumbent +fw_verbose - if true, FrankWolfe logs are printed """ function solve( f, @@ -36,6 +46,9 @@ function solve( lmo; traverse_strategy=Bonobo.BestFirstSearch(), branching_strategy=Bonobo.MOST_INFEASIBLE(), + variant::FrankWolfeVariant=BPCG(), + line_search::FrankWolfe.LineSearchMethod=FrankWolfe.Adaptive(), + active_set::Union{Nothing, FrankWolfe.ActiveSet} = nothing, fw_epsilon=1e-2, verbose=false, dual_gap=1e-6, @@ -53,7 +66,9 @@ function solve( global_dual_tightening=true, bnb_callback=nothing, strong_convexity=0.0, + domain_oracle= x->true, start_solution=nothing, + fw_verbose = false, kwargs..., ) if verbose @@ -61,6 +76,8 @@ function solve( println("Parameter settings.") println("\t Tree traversal strategy: ", _value_to_print(traverse_strategy)) println("\t Branching strategy: ", _value_to_print(branching_strategy)) + println("\t FrankWolfe variant: $(variant)") + println("\t Line Search Method: $(line_search)") @printf("\t Absolute dual gap tolerance: %e\n", dual_gap) @printf("\t Relative dual gap tolerance: %e\n", rel_dual_gap) @printf("\t Frank-Wolfe subproblem tolerance: %e\n", fw_epsilon) @@ -136,10 +153,20 @@ function solve( global_bounds[idx, :lessthan] = MOI.LessThan(1.0) end - direction = collect(1.0:n) - v = compute_extreme_point(lmo, direction) - v[integer_variables] = round.(v[integer_variables]) - active_set = FrankWolfe.ActiveSet([(1.0, v)]) + v = [] + if active_set === nothing + direction = collect(1.0:n) + v = compute_extreme_point(lmo, direction) + v[integer_variables] = round.(v[integer_variables]) + active_set = FrankWolfe.ActiveSet([(1.0, v)]) + vertex_storage = FrankWolfe.DeletedVertexStorage(typeof(v)[], 1) + else + @assert FrankWolfe.active_set_validate(active_set) + for a in active_set.atoms + @assert is_linear_feasible(lmo.o, a) + end + v = active_set.atoms[1] + end vertex_storage = FrankWolfe.DeletedVertexStorage(typeof(v)[], 1) m = Boscia.SimpleOptimizationProblem(f, grad!, n, integer_variables, time_lmo, global_bounds) @@ -183,6 +210,11 @@ function solve( :dual_tightening => dual_tightening, :global_dual_tightening => global_dual_tightening, :strong_convexity => strong_convexity, + :variant => variant, + :lineSearch => line_search, + :domain_oracle => domain_oracle, + :usePostsolve => use_postsolve, + :fwVerbose => fw_verbose, ), ), branch_strategy=branching_strategy, @@ -269,22 +301,24 @@ function solve( Bonobo.optimize!(tree; callback=bnb_callback) - x = postsolve(tree, result, time_ref, verbose, use_postsolve, max_iteration_post) + x = postsolve(tree, result, time_ref, verbose, max_iteration_post) # Check solution and polish x_polished = x - if !is_linear_feasible(tree.root.problem.lmo, x) - error("Reported solution not linear feasbile!") - end - if !is_integer_feasible(tree.root.problem.integer_variables, x, atol=1e-16, rtol=1e-16) - @info "Polish solution" - for i in tree.root.problem.integer_variables - x_polished[i] = round(x_polished[i]) + if x !== nothing + if !is_linear_feasible(tree.root.problem.lmo, x) + error("Reported solution not linear feasbile!") end - if !is_linear_feasible(tree.root.problem.lmo, x_polished) - @warn "Polished solution not linear feasible" - else - x = x_polished + if !is_integer_feasible(tree.root.problem.integer_variables, x, atol=1e-16, rtol=1e-16) && x !== nothing + @info "Polish solution" + for i in tree.root.problem.integer_variables + x_polished[i] = round(x_polished[i]) + end + if !is_linear_feasible(tree.root.problem.lmo, x_polished) + @warn "Polished solution not linear feasible" + else + x = x_polished + end end end println() # cleaner output @@ -492,7 +526,13 @@ function build_bnb_callback( if Bonobo.terminated(tree) Bonobo.sort_solutions!(tree.solutions, tree.sense) x = Bonobo.get_solution(tree) - primal_value = tree.root.problem.f(x) + # x can be nothing if the user supplied a custom domain oracle and the time limit is reached + if x === nothing + @assert tree.root.problem.solving_stage == TIME_LIMIT_REACHED + end + primal_value = x !== nothing ? tree.root.problem.f(x) : Inf + # deactivate postsolve if there is no solution + tree.root.options[:usePostsolve] = x === nothing ? false : tree.root.options[:usePostsolve] # TODO: here we need to calculate the actual state @@ -533,9 +573,9 @@ Runs the post solve both for a cleaner solutiona and to optimize for the continuous variables if present. Prints solution statistics if verbose is true. """ -function postsolve(tree, result, time_ref, verbose, use_postsolve, max_iteration_post) +function postsolve(tree, result, time_ref, verbose, max_iteration_post) x = Bonobo.get_solution(tree) - primal = tree.incumbent_solution.objective + primal = x !== nothing ? tree.incumbent_solution.objective : Inf status_string = "FIX ME" # should report "feasible", "optimal", "infeasible", "gap tolerance met" if isempty(tree.nodes) @@ -549,7 +589,7 @@ function postsolve(tree, result, time_ref, verbose, use_postsolve, max_iteration end only_integer_vars = tree.root.problem.nvars == length(tree.root.problem.integer_variables) - if use_postsolve && !only_integer_vars + if tree.root.options[:usePostsolve] && !only_integer_vars # Build solution lmo fix_bounds = IntegerBounds() for i in tree.root.problem.integer_variables diff --git a/src/node.jl b/src/node.jl index f066bda3f..817200a92 100644 --- a/src/node.jl +++ b/src/node.jl @@ -72,11 +72,11 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN end new_bound_left = lower_bound_base + μ/2 * (x[vidx] - floor(x[vidx]))^2 new_bound_right = lower_bound_base + μ/2 * (ceil(x[vidx]) - x[vidx])^2 - if new_bound_left >= tree.incumbent + if new_bound_left > tree.incumbent @debug "prune left, from $(node.lb) -> $new_bound_left, ub $(tree.incumbent), lb $(node.lb)" prune_left = true end - if new_bound_right >= tree.incumbent + if new_bound_right > tree.incumbent @debug "prune right, from $(node.lb) -> $new_bound_right, ub $(tree.incumbent), lb $(node.lb)" prune_right = true end @@ -157,11 +157,21 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN local_potential_tightenings=0, dual_gap=NaN, ) - nodes = if !prune_left && !prune_right + + # in case of non trivial domain oracle: Only split if the iterate is still domain feasible + x_left = FrankWolfe.compute_active_set_iterate!(active_set_left) + x_right = FrankWolfe.compute_active_set_iterate!(active_set_right) + domain_oracle = tree.root.options[:domain_oracle] + + nodes = if !prune_left && !prune_right #&& domain_oracle(x_left) && domain_oracle(x_right) [node_info_left, node_info_right] - elseif prune_left + elseif prune_left [node_info_right] - else + elseif prune_right + [node_info_left] + elseif domain_oracle(x_right) + [node_info_right] + elseif domain_oracle(x_left) [node_info_left] end return nodes @@ -262,21 +272,63 @@ function Bonobo.evaluate_node!(tree::Bonobo.BnBTree, node::FrankWolfeNode) end end - # call blended_pairwise_conditional_gradient - x, _, primal, dual_gap, _, active_set = FrankWolfe.blended_pairwise_conditional_gradient( + # x = zeros(tree.root.problem.nvars) + #dual_gap = primal = 0 + #active_set = FrankWolfe.ActiveSet([(1.0, x)]) + x=FrankWolfe.compute_active_set_iterate!(node.active_set) + domain_oracle = tree.root.options[:domain_oracle] + + x, primal, dual_gap, active_set = solve_frank_wolfe( + tree.root.options[:variant], tree.root.problem.f, tree.root.problem.g, tree.root.problem.lmo, - node.active_set, + node.active_set; epsilon=node.fw_dual_gap_limit, max_iteration=tree.root.options[:max_fw_iter], - line_search=FrankWolfe.Adaptive(verbose=false), + line_search=tree.root.options[:lineSearch], + lazy=true, + timeout=tree.root.options[:time_limit], add_dropped_vertices=true, use_extra_vertex_storage=true, extra_vertex_storage=node.discarded_vertices, callback=tree.root.options[:callback], - lazy=true, - ) + verbose=tree.root.options[:fwVerbose], + )#= + if tree.root.options[:variant] == BPCG + # call blended_pairwise_conditional_gradient + x, _, primal, dual_gap, _, active_set = FrankWolfe.blended_pairwise_conditional_gradient( + tree.root.problem.f, + tree.root.problem.g, + tree.root.problem.lmo, + node.active_set, + epsilon=node.fw_dual_gap_limit, + max_iteration=tree.root.options[:max_fw_iter], + line_search=tree.root.options[:lineSearch], + add_dropped_vertices=true, + use_extra_vertex_storage=true, + extra_vertex_storage=node.discarded_vertices, + callback=tree.root.options[:callback], + lazy=true, + timeout=tree.root.options[:time_limit], + verbose=tree.root.options[:fwVerbose], + ) + elseif tree.root.options[:variant] == AFW + # call away_frank_wolfe + x, _, primal, dual_gap, _, active_set = FrankWolfe.away_frank_wolfe( + tree.root.problem.f, + tree.root.problem.g, + tree.root.problem.lmo, + node.active_set, + epsilon=node.fw_dual_gap_limit, + max_iteration=tree.root.options[:max_fw_iter], + line_search=tree.root.options[:lineSearch], + callback=tree.root.options[:callback], + lazy=true, + timeout=tree.root.options[:time_limit], + verbose=tree.root.options[:fwVerbose], + ) + end=# node.fw_time = Dates.now() - time_ref node.dual_gap = dual_gap diff --git a/test/interface_test.jl b/test/interface_test.jl index 867fa2ad0..d3b5039f0 100644 --- a/test/interface_test.jl +++ b/test/interface_test.jl @@ -2,6 +2,7 @@ using LinearAlgebra using Distributions import Random using SCIP +using HiGHS import MathOptInterface const MOI = MathOptInterface import Boscia @@ -170,3 +171,107 @@ Ns = 0.1 @test sum(x2[p+1:2p]) <= k @test f(x2) == f(x) end + +n = 20 +diffi = Random.rand(Bool, n) * 0.6 .+ 0.3 + +@testset "Different FW variants" begin + + function build_model() + o = SCIP.Optimizer() + MOI.set(o, MOI.Silent(), true) + MOI.empty!(o) + x = MOI.add_variables(o, n) + for xi in x + MOI.add_constraint(o, xi, MOI.GreaterThan(0.0)) + MOI.add_constraint(o, xi, MOI.LessThan(1.0)) + MOI.add_constraint(o, xi, MOI.ZeroOne()) # or MOI.Integer() + end + lmo = FrankWolfe.MathOptLMO(o) + return lmo + end + + function f(x) + return 0.5 * sum((x[i] - diffi[i])^2 for i in eachindex(x)) + end + function grad!(storage, x) + @. storage = x - diffi + end + + direction = rand(n) + + lmo = build_model() + v = FrankWolfe.compute_extreme_point(lmo, direction) + x_afw, _, result_afw = Boscia.solve(f, grad!, lmo, verbose=false, variant=Boscia.AwayFrankWolfe()) + + lmo = build_model() + v = FrankWolfe.compute_extreme_point(lmo, direction) + x_blended, _, result_blended = Boscia.solve(f, grad!, lmo, verbose=false, variant=Boscia.Blended()) + + lmo = build_model() + v = FrankWolfe.compute_extreme_point(lmo, direction) + x_bpcg, _, result_bpcg = Boscia.solve(f, grad!, lmo, verbose=false, variant=Boscia.BPCG()) + + lmo = build_model() + v = FrankWolfe.compute_extreme_point(lmo, direction) + x_vfw, _, result_vfw = Boscia.solve(f, grad!, lmo, verbose=false, variant=Boscia.VanillaFrankWolfe()) + + @test isapprox(f(x_afw), f(result_afw[:raw_solution]), atol=1e-6, rtol=1e-3) + @test isapprox(f(x_blended), f(result_blended[:raw_solution]), atol=1e-6, rtol=1e-3) + @test isapprox(f(x_bpcg), f(result_bpcg[:raw_solution]), atol=1e-6, rtol=1e-3) + @test isapprox(f(x_vfw), f(result_vfw[:raw_solution]), atol=1e-6, rtol=1e-3) + + @test sum(isapprox.(x_afw, x_blended, atol=1e-6, rtol=1e-3)) == n + @test sum(isapprox.(x_blended, x_bpcg, atol=1e-6, rtol=1e-3)) == n + @test sum(isapprox.(x_bpcg, x_vfw, atol=1e-6, rtol=1e-3)) == n + @test sum(isapprox.(x_vfw, x_afw, atol=1e-6, rtol=1e-3)) == n +end + +@testset "Different line search types" begin + + function build_model() + o = SCIP.Optimizer() + MOI.set(o, MOI.Silent(), true) + MOI.empty!(o) + x = MOI.add_variables(o, n) + for xi in x + MOI.add_constraint(o, xi, MOI.GreaterThan(0.0)) + MOI.add_constraint(o, xi, MOI.LessThan(1.0)) + MOI.add_constraint(o, xi, MOI.ZeroOne()) # or MOI.Integer() + end + lmo = FrankWolfe.MathOptLMO(o) + return lmo + end + + function f(x) + return 0.5 * sum((x[i] - diffi[i])^2 for i in eachindex(x)) + end + function grad!(storage, x) + @. storage = x - diffi + end + + direction = rand(n) + + lmo = build_model() + v = FrankWolfe.compute_extreme_point(lmo, direction) + line_search = FrankWolfe.Adaptive() + x_adaptive, _, result_adaptive = Boscia.solve(f, grad!, lmo, verbose=false, line_search=line_search) + + lmo = build_model() + v = FrankWolfe.compute_extreme_point(lmo, direction) + line_search = FrankWolfe.MonotonicStepSize() + x_monotonic, _, result_monotonic = Boscia.solve(f, grad!, lmo, verbose=false, line_search=line_search, time_limit=120) + + lmo = build_model() + v = FrankWolfe.compute_extreme_point(lmo, direction) + line_search = FrankWolfe.Agnostic() + x_agnostic, _, result_agnostic = Boscia.solve(f, grad!, lmo, verbose=false, line_search=line_search, time_limit=120) + + @test isapprox(f(x_adaptive), f(result_adaptive[:raw_solution]), atol=1e-6, rtol=1e-3) + @test isapprox(f(x_monotonic), f(result_monotonic[:raw_solution]), atol=1e-6, rtol=1e-3) + @test isapprox(f(x_agnostic), f(result_agnostic[:raw_solution]), atol=1e-6, rtol=1e-3) + + @test sum(isapprox.(x_adaptive, x_monotonic, atol=1e-6, rtol=1e-3)) == n + @test sum(isapprox.(x_agnostic, x_monotonic, atol=1e-6, rtol=1e-3)) == n + @test sum(isapprox.(x_adaptive, x_agnostic, atol=1e-6, rtol=1e-3)) == n +end \ No newline at end of file diff --git a/test/mean_risk.jl b/test/mean_risk.jl index 178e653c9..c51eb53b6 100644 --- a/test/mean_risk.jl +++ b/test/mean_risk.jl @@ -61,6 +61,6 @@ const M1 = (A1 + A1') / 2 x,_,result = Boscia.solve(f, grad!, lmo, verbose= true) - @test dot(a,x) <= b + 1e-6 + @test dot(a,x) <= b + 1e-4 @test f(x) <= f(result[:raw_solution]) + 1e-6 end diff --git a/test/runtests.jl b/test/runtests.jl index d3e0966f7..df0d92e79 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -131,6 +131,7 @@ end @test isapprox(x, round.(diff1), atol=1e-5, rtol=1e-5) end + for file in readdir(joinpath(@__DIR__, "../examples/"), join=true) if endswith(file, "jl") include(file)