From 5adbb2affb6263effcd0706de0a35d8cd0ca2a50 Mon Sep 17 00:00:00 2001 From: Deborah Hendrych Date: Mon, 23 Mar 2026 22:20:56 +0100 Subject: [PATCH 1/8] Extend to Boscia functionality. --- src/spanning_tree.jl | 224 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 224 insertions(+) diff --git a/src/spanning_tree.jl b/src/spanning_tree.jl index 49dc269..6e857fe 100644 --- a/src/spanning_tree.jl +++ b/src/spanning_tree.jl @@ -9,6 +9,36 @@ struct SpanningTreeLMO{G} <: FrankWolfe.LinearMinimizationOracle graph::G end +""" +Union-find find with path compression. +Used to detect cycles and connectivity in fixing checks. +""" +function uf_find!(parent::Vector{Int}, x::Int) + while parent[x] != x + parent[x] = parent[parent[x]] + x = parent[x] + end + return x +end + +""" +Union-find union. +Returns `false` when `a` and `b` are already connected (cycle). +""" +function uf_union!(parent::Vector{Int}, a::Int, b::Int) + ra = uf_find!(parent, a) + rb = uf_find!(parent, b) + if ra == rb + return false + end + parent[rb] = ra + return true +end + +""" +Minimum spanning tree LMO using Kruskal on the weighted graph. +Returns an incidence vector over `edges(g)`. +""" function FrankWolfe.compute_extreme_point( lmo::SpanningTreeLMO, direction::M; @@ -34,3 +64,197 @@ function FrankWolfe.compute_extreme_point( end return v end + +""" +Bound-aware LMO for spanning trees. +Contracts forced edges, removes forbidden edges, and runs Kruskal +on the reduced graph, then lifts the solution to the original graph. +""" +function Boscia.bounded_compute_extreme_point( + lmo::SpanningTreeLMO, + direction, + lb, + ub, + int_vars; + kwargs..., +) + N = length(direction) + edges_iter = collect(Graphs.edges(lmo.graph)) + @assert length(edges_iter) == N + check_spanning_tree_fixings(lmo.graph, edges_iter, lb, ub) + + # Contract all forced edges into components via union-find. + # Connected nodes will have the same parent node at the end. + parent = collect(1:Graphs.nv(lmo.graph)) + for (i, edge) in enumerate(edges_iter) + if lb[i] ≈ 1 + @assert ub[i] ≈ 1 + uf_union!(parent, src(edge), dst(edge)) + end + end + # Count the number of connected components + # which will be contracted to super nodes. + comp_id = Dict{Int,Int}() + comp = Vector{Int}(undef, Graphs.nv(lmo.graph)) + k = 0 + for vtx in 1:Graphs.nv(lmo.graph) + root = uf_find!(parent, vtx) + if !haskey(comp_id, root) + k += 1 + comp_id[root] = k + end + comp[vtx] = comp_id[root] + end + + # Initialize solution with all forced edges. + v = spzeros(N) + for i in 1:N + if lb[i] ≈ 1 + v[i] = 1 + end + end + + # Build reduced graph on components using the cheapest allowed + # edge for each component pair. (Note that after contracted the graph + # we might have multiple edges between the same super nodes.) + # Then run Kruskal on the reduced graph. + if k > 1 + edge_choice = Dict{Tuple{Int,Int},Tuple{eltype(direction),Int}}() + for (i, edge) in enumerate(edges_iter) + # Edge is forbidden. + if ub[i] ≈ 0 + continue + end + c1 = comp[src(edge)] + c2 = comp[dst(edge)] + # source and destination nodes lie in the same component, + # so the edge can be ignored. + if c1 == c2 + continue + end + # order independent key + if c1 > c2 + c1, c2 = c2, c1 + end + key = (c1, c2) + w = direction[i] + # If multiple edges connect the same super nodes + # choose the cheapest one. + if !haskey(edge_choice, key) || w < edge_choice[key][1] + edge_choice[key] = (w, i) + end + end + # Build reduced graph and weight matrix. + reduced_graph = SimpleGraph(k) + for key in keys(edge_choice) + add_edge!(reduced_graph, key[1], key[2]) + end + distmx = spzeros(k, k) + for (key, (w, _)) in edge_choice + distmx[key[1], key[2]] = w + distmx[key[2], key[1]] = w + end + reduced_span = Graphs.kruskal_mst(reduced_graph, distmx) + for edge in reduced_span + c1 = src(edge) + c2 = dst(edge) + # order independent key + if c1 > c2 + c1, c2 = c2, c1 + end + idx = edge_choice[(c1, c2)][2] + v[idx] = 1 + end + end + # Optional sanity check. + @debug begin + for i in 1:N + if ub[i] ≈ 0 + @assert v[i] ≈ 0 + elseif lb[i] ≈ 1 + @assert v[i] ≈ 1 + end + end + end + return v +end + +""" +Lightweight feasibility check for a candidate `v`. +Enforces total edge count and singleton-cut constraints. +""" +function Boscia.is_simple_linear_feasible(lmo::SpanningTreeLMO, v) + n = Graphs.nv(lmo.graph) + if n == 0 + return true + end + total = sum(v) + if abs(total - (n - 1)) > 1e-4 + return false + end + # detect cycles among edges that are (almost) fully selected + parent = collect(1:n) + for (idx, edge) in enumerate(edges(lmo.graph)) + if v[idx] < 1 - 1e-4 + continue + end + if !(uf_union!(parent, src(edge), dst(edge))) + return false + end + end + # singleton cut constraints: each vertex must have at least one incident edge + degrees = zeros(eltype(v), n) + for (idx, edge) in enumerate(edges(lmo.graph)) + if v[idx] ≈ 0 + continue + end + degrees[src(edge)] += v[idx] + degrees[dst(edge)] += v[idx] + end + if minimum(degrees) < 1 - 1e-4 + return false + end + return true +end + +""" +Feasibility of bounds alone for spanning trees. +Returns `Boscia.OPTIMAL` if some spanning tree can satisfy the bounds. +""" +function Boscia.check_feasibility( + lmo::SpanningTreeLMO, + lb, + ub, + int_vars, + n +) + edges_iter = collect(Graphs.edges(lmo.graph)) + if n <= 1 + return Boscia.OPTIMAL + end + # The forced edges (lb=ub=1) must be acyclic. + parent = collect(1:n) + for (i, edge) in enumerate(edges_iter) + if lb[i] ≈ 1 + if !uf_union!(parent, src(edge), dst(edge)) + @debug "Forced edges form a cycle" + return Boscia.INFEASIBLE + end + end + end + # The graph must stay connected after removing forbidden edges. + parent = collect(1:n) + for (i, edge) in enumerate(edges_iter) + if !(ub[i] ≈ 0) + uf_union!(parent, src(edge), dst(edge)) + end + end + root = uf_find!(parent, 1) + for vtx in 2:n + if uf_find!(parent, vtx) != root + @debug "Forbidden edges disconnect graph" + return Boscia.INFEASIBLE + end + end + return Boscia.OPTIMAL +end From 2b03ca45c64964fc81df91cdc86f5f710124daed Mon Sep 17 00:00:00 2001 From: Deborah Hendrych Date: Mon, 23 Mar 2026 22:28:44 +0100 Subject: [PATCH 2/8] Test. --- src/spanning_tree.jl | 1 - test/runtests.jl | 51 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 51 insertions(+), 1 deletion(-) diff --git a/src/spanning_tree.jl b/src/spanning_tree.jl index 6e857fe..204209b 100644 --- a/src/spanning_tree.jl +++ b/src/spanning_tree.jl @@ -81,7 +81,6 @@ function Boscia.bounded_compute_extreme_point( N = length(direction) edges_iter = collect(Graphs.edges(lmo.graph)) @assert length(edges_iter) == N - check_spanning_tree_fixings(lmo.graph, edges_iter, lb, ub) # Contract all forced edges into components via union-find. # Connected nodes will have the same parent node at the end. diff --git a/test/runtests.jl b/test/runtests.jl index ecaf42a..f9802e1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -145,6 +145,57 @@ end v = FrankWolfe.compute_extreme_point(lmo, direction) @test dot(v, direction) < -4e-7 end + + @testset "Bounded and feasibility checks" begin + g4 = Graphs.complete_graph(4) + lmo4 = CO.SpanningTreeLMO(g4) + iter4 = collect(Graphs.edges(g4)) + M4 = length(iter4) + + idx12 = findfirst(==(Edge(1, 2)), iter4) + idx23 = findfirst(==(Edge(2, 3)), iter4) + idx13 = findfirst(==(Edge(1, 3)), iter4) + idx14 = findfirst(==(Edge(1, 4)), iter4) + idx24 = findfirst(==(Edge(2, 4)), iter4) + idx34 = findfirst(==(Edge(3, 4)), iter4) + + # force edges (1,2) and (2,3); best connecting edge should be (3,4) + direction = ones(M4) + direction[idx14] = 5.0 + direction[idx24] = 2.0 + direction[idx34] = -1.0 + lb = zeros(M4) + ub = ones(M4) + lb[idx12] = 1.0 + lb[idx23] = 1.0 + v = Boscia.bounded_compute_extreme_point(lmo4, direction, lb, ub, 1:M4) + @test v[idx12] == 1 + @test v[idx23] == 1 + @test v[idx34] == 1 + + # forbid the cheapest edge (1,4); must choose the next best connection + direction[idx14] = -10.0 + ub[idx14] = 0.0 + v2 = Boscia.bounded_compute_extreme_point(lmo4, direction, lb, ub, 1:M4) + @test v2[idx14] == 0 + @test v2[idx24] == 1 || v2[idx34] == 1 + + # feasibility: cycle in forced edges + lb_cycle = zeros(M4) + ub_cycle = ones(M4) + lb_cycle[idx12] = 1.0 + lb_cycle[idx23] = 1.0 + lb_cycle[idx13] = 1.0 + @test Boscia.check_feasibility(lmo4, lb_cycle, ub_cycle, 1:M4, M4) == Boscia.INFEASIBLE + + # feasibility: disconnect node 4 + lb_disc = zeros(M4) + ub_disc = ones(M4) + ub_disc[idx14] = 0.0 + ub_disc[idx24] = 0.0 + ub_disc[idx34] = 0.0 + @test Boscia.check_feasibility(lmo4, lb_disc, ub_disc, 1:M4, M4) == Boscia.INFEASIBLE + end end @testset "Shortest path" begin From 9855a58fc634da2645d72a8e73d6c0fb595852f3 Mon Sep 17 00:00:00 2001 From: Deborah Hendrych Date: Tue, 24 Mar 2026 17:09:31 +0100 Subject: [PATCH 3/8] Spanning trees should also be connected. --- src/spanning_tree.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/spanning_tree.jl b/src/spanning_tree.jl index 204209b..87e88ab 100644 --- a/src/spanning_tree.jl +++ b/src/spanning_tree.jl @@ -213,6 +213,21 @@ function Boscia.is_simple_linear_feasible(lmo::SpanningTreeLMO, v) if minimum(degrees) < 1 - 1e-4 return false end + # ensure support is connected (prevents disjoint forests passing) + parent = collect(1:n) + for (idx, edge) in enumerate(edges(lmo.graph)) + if v[idx] <= 1e-4 + continue + end + uf_union!(parent, src(edge), dst(edge)) + end + # All nodes should have a common root if the graph is connected. + root = uf_find!(parent, 1) + for vtx in 2:n + if uf_find!(parent, vtx) != root + return false + end + end return true end From 8f32b0dbdd98af34088843a09d6747e44b317104 Mon Sep 17 00:00:00 2001 From: Deborah Hendrych Date: Tue, 24 Mar 2026 19:28:05 +0100 Subject: [PATCH 4/8] Functions for DICG. --- src/spanning_tree.jl | 99 ++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 77 ++++++++++++++++++++++++++++++++++ 2 files changed, 176 insertions(+) diff --git a/src/spanning_tree.jl b/src/spanning_tree.jl index 87e88ab..bee3313 100644 --- a/src/spanning_tree.jl +++ b/src/spanning_tree.jl @@ -272,3 +272,102 @@ function Boscia.check_feasibility( end return Boscia.OPTIMAL end + +""" +Spanning trees support decomposition-invariant oracles since in-faces +are defined by fixing edges to 0/1. +""" +function Boscia.is_decomposition_invariant_oracle_simple(::SpanningTreeLMO) + return true +end + +function FrankWolfe.is_decomposition_invariant_oracle(::SpanningTreeLMO) + return true +end + +""" +In-face bounded LMO: adds fixings implied by `x` (0/1 entries) +on top of existing bounds, then solves the bounded LMO. +""" +function Boscia.bounded_compute_inface_extreme_point( + lmo::SpanningTreeLMO, + direction, + x, + lb, + ub, + int_vars; + kwargs..., +) + N = length(direction) + lb2 = copy(lb) + ub2 = copy(ub) + for i in 1:N + if x[i] ≤ eps() + ub2[i] = 0.0 + elseif x[i] ≥ 1 - eps() + lb2[i] = 1.0 + end + end + return Boscia.bounded_compute_extreme_point(lmo, direction, lb2, ub2, int_vars; kwargs...) +end + +""" +Check whether `a` lies on the minimal face of `x` and respects bounds. +""" +function Boscia.is_simple_inface_feasible( + lmo::SpanningTreeLMO, + a, + x, + lb, + ub, + int_vars, +) + for i in eachindex(a) + if x[i] ≤ eps() && a[i] > 1e-6 + return false + elseif x[i] ≥ 1 - eps() && a[i] < 1 - 1e-6 + return false + end + end + return true +end + +""" +Maximum DICG step size respecting bounds and in-face fixings from `x`. +""" +function Boscia.bounded_dicg_maximum_step( + lmo::SpanningTreeLMO, + x, + direction, + lb, + ub, + int_vars; + kwargs..., +) + T = promote_type(eltype(x), eltype(direction)) + gamma_max = one(T) + for i in eachindex(x) + if direction[i] == 0 + continue + end + lower = lb[i] + upper = ub[i] + if x[i] ≤ eps() + upper = min(upper, zero(T)) + elseif x[i] ≥ 1 - eps() + lower = max(lower, one(T)) + end + if direction[i] > 0 + if x[i] ≥ upper - 1e-12 + return zero(gamma_max) + end + gamma_max = min(gamma_max, (upper - x[i]) / direction[i]) + else + if x[i] ≤ lower + 1e-12 + return zero(gamma_max) + end + gamma_max = min(gamma_max, (lower - x[i]) / direction[i]) + end + end + return gamma_max +end diff --git a/test/runtests.jl b/test/runtests.jl index f9802e1..743da1e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -127,6 +127,7 @@ end Random.seed!(1645) g = Graphs.complete_graph(N) lmo = CO.SpanningTreeLMO(g) + @test FrankWolfe.is_decomposition_invariant_oracle(lmo) iter = collect(Graphs.edges(g)) M = length(iter) @testset "Basic tree properties" begin @@ -149,6 +150,7 @@ end @testset "Bounded and feasibility checks" begin g4 = Graphs.complete_graph(4) lmo4 = CO.SpanningTreeLMO(g4) + @test Boscia.is_decomposition_invariant_oracle_simple(lmo4) iter4 = collect(Graphs.edges(g4)) M4 = length(iter4) @@ -196,6 +198,80 @@ end ub_disc[idx34] = 0.0 @test Boscia.check_feasibility(lmo4, lb_disc, ub_disc, 1:M4, M4) == Boscia.INFEASIBLE end + + @testset "bounded_dicg_maximum_step" begin + g4 = Graphs.complete_graph(4) + lmo4 = CO.SpanningTreeLMO(g4) + iter4 = collect(Graphs.edges(g4)) + M4 = length(iter4) + + idx12 = findfirst(==(Edge(1, 2)), iter4) + idx23 = findfirst(==(Edge(2, 3)), iter4) + + lb = zeros(M4) + ub = ones(M4) + + # If x is fixed at 0 and direction wants to increase, step is zero. + x = zeros(M4) + direction = zeros(M4) + direction[idx12] = 1.0 + γ0 = Boscia.bounded_dicg_maximum_step(lmo4, x, direction, lb, ub, 1:M4) + @test γ0 == 0.0 + + # Otherwise, gamma is limited by the tightest bound among active directions. + x = zeros(M4) + direction = zeros(M4) + x[idx12] = 0.2 + x[idx23] = 0.8 + direction[idx12] = 1.0 + direction[idx23] = -2.0 + ub[idx12] = 0.9 + lb[idx23] = 0.1 + γ = Boscia.bounded_dicg_maximum_step(lmo4, x, direction, lb, ub, 1:M4) + @test isapprox(γ, 0.35, atol=1e-12, rtol=0.0) + end + + @testset "Decomposition-invariant in-face" begin + g4 = Graphs.complete_graph(4) + lmo4 = CO.SpanningTreeLMO(g4) + iter4 = collect(Graphs.edges(g4)) + M4 = length(iter4) + + idx12 = findfirst(==(Edge(1, 2)), iter4) + idx23 = findfirst(==(Edge(2, 3)), iter4) + idx13 = findfirst(==(Edge(1, 3)), iter4) + idx14 = findfirst(==(Edge(1, 4)), iter4) + idx24 = findfirst(==(Edge(2, 4)), iter4) + idx34 = findfirst(==(Edge(3, 4)), iter4) + + direction = ones(M4) + direction[idx14] = 1.0 + direction[idx24] = 2.0 + direction[idx34] = -5.0 + + x = fill(0.3, M4) + x[idx12] = 1.0 + x[idx23] = 1.0 + x[idx13] = 0.0 + + lb = zeros(M4) + ub = ones(M4) + v_if = Boscia.bounded_compute_inface_extreme_point(lmo4, direction, x, lb, ub, 1:M4) + + @test v_if[idx12] == 1.0 + @test v_if[idx23] == 1.0 + @test v_if[idx13] == 0.0 + @test v_if[idx34] == 1.0 + @test Boscia.is_simple_inface_feasible(lmo4, v_if, x, lb, ub, 1:M4) + + v_bad = copy(v_if) + v_bad[idx13] = 1.0 + @test !Boscia.is_simple_inface_feasible(lmo4, v_bad, x, lb, ub, 1:M4) + + v_bad2 = copy(v_if) + v_bad2[idx12] = 0.0 + @test !Boscia.is_simple_inface_feasible(lmo4, v_bad2, x, lb, ub, 1:M4) + end end @testset "Shortest path" begin @@ -238,6 +314,7 @@ end n = 4 d = randn(rng, n, n) lmo = CO.BirkhoffLMO(n) + @test FrankWolfe.is_decomposition_invariant_oracle(lmo) x = ones(n, n) ./ n # test without fixings v_if = CO.compute_inface_extreme_point(lmo, d, x) From 374df87ac7ee6832bb7cc5bf84c785418d18fa33 Mon Sep 17 00:00:00 2001 From: Deborah Hendrych Date: Tue, 24 Mar 2026 20:29:11 +0100 Subject: [PATCH 5/8] minor. --- src/spanning_tree.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spanning_tree.jl b/src/spanning_tree.jl index bee3313..0ca8d5d 100644 --- a/src/spanning_tree.jl +++ b/src/spanning_tree.jl @@ -80,7 +80,7 @@ function Boscia.bounded_compute_extreme_point( ) N = length(direction) edges_iter = collect(Graphs.edges(lmo.graph)) - @assert length(edges_iter) == N + @assert length(edges_iter) == N "length(edges_iter) = $(length(edges_iter)) != N = $N" # Contract all forced edges into components via union-find. # Connected nodes will have the same parent node at the end. From 714b7624c110c1691535a0391ba4015fe7e96b3d Mon Sep 17 00:00:00 2001 From: Deborah Hendrych Date: Tue, 24 Mar 2026 21:29:02 +0100 Subject: [PATCH 6/8] Minor change. --- src/spanning_tree.jl | 10 ++++++---- test/runtests.jl | 4 ++-- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/spanning_tree.jl b/src/spanning_tree.jl index 0ca8d5d..421dd17 100644 --- a/src/spanning_tree.jl +++ b/src/spanning_tree.jl @@ -243,11 +243,13 @@ function Boscia.check_feasibility( n ) edges_iter = collect(Graphs.edges(lmo.graph)) - if n <= 1 + n_local = Graphs.nv(lmo.graph) + @debug "n_local = $n_local n = $n" + if n_local <= 1 return Boscia.OPTIMAL end # The forced edges (lb=ub=1) must be acyclic. - parent = collect(1:n) + parent = collect(1:n_local) for (i, edge) in enumerate(edges_iter) if lb[i] ≈ 1 if !uf_union!(parent, src(edge), dst(edge)) @@ -257,14 +259,14 @@ function Boscia.check_feasibility( end end # The graph must stay connected after removing forbidden edges. - parent = collect(1:n) + parent = collect(1:n_local) for (i, edge) in enumerate(edges_iter) if !(ub[i] ≈ 0) uf_union!(parent, src(edge), dst(edge)) end end root = uf_find!(parent, 1) - for vtx in 2:n + for vtx in 2:n_local if uf_find!(parent, vtx) != root @debug "Forbidden edges disconnect graph" return Boscia.INFEASIBLE diff --git a/test/runtests.jl b/test/runtests.jl index 743da1e..f259608 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -188,7 +188,7 @@ end lb_cycle[idx12] = 1.0 lb_cycle[idx23] = 1.0 lb_cycle[idx13] = 1.0 - @test Boscia.check_feasibility(lmo4, lb_cycle, ub_cycle, 1:M4, M4) == Boscia.INFEASIBLE + @test Boscia.check_feasibility(lmo4, lb_cycle, ub_cycle, 1:M4, nv(g4)) == Boscia.INFEASIBLE # feasibility: disconnect node 4 lb_disc = zeros(M4) @@ -196,7 +196,7 @@ end ub_disc[idx14] = 0.0 ub_disc[idx24] = 0.0 ub_disc[idx34] = 0.0 - @test Boscia.check_feasibility(lmo4, lb_disc, ub_disc, 1:M4, M4) == Boscia.INFEASIBLE + @test Boscia.check_feasibility(lmo4, lb_disc, ub_disc, 1:M4, nv(g4)) == Boscia.INFEASIBLE end @testset "bounded_dicg_maximum_step" begin From 8134749413c5768ccc08dcc30a78bc1bdd8f8c69 Mon Sep 17 00:00:00 2001 From: Deborah Hendrych Date: Wed, 25 Mar 2026 14:17:13 +0100 Subject: [PATCH 7/8] Convention fix. --- src/spanning_tree.jl | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/spanning_tree.jl b/src/spanning_tree.jl index 421dd17..fdfba4c 100644 --- a/src/spanning_tree.jl +++ b/src/spanning_tree.jl @@ -354,21 +354,22 @@ function Boscia.bounded_dicg_maximum_step( end lower = lb[i] upper = ub[i] - if x[i] ≤ eps() - upper = min(upper, zero(T)) - elseif x[i] ≥ 1 - eps() - lower = max(lower, one(T)) - end + # FrankWolfe updates as: x_new = x - gamma * direction. + # We need gamma >= 0 such that lower <= x_new <= upper. if direction[i] > 0 - if x[i] ≥ upper - 1e-12 + # x_new decreases with gamma; only the lower bound can become active: + # lower <= x[i] - gamma*dir => gamma <= (x[i] - lower)/dir + if x[i] ≤ lower + 1e-12 return zero(gamma_max) end - gamma_max = min(gamma_max, (upper - x[i]) / direction[i]) + gamma_max = min(gamma_max, (x[i] - lower) / direction[i]) else - if x[i] ≤ lower + 1e-12 + # direction[i] < 0: x_new increases with gamma; only the upper bound can become active: + # x[i] - gamma*dir <= upper => gamma <= (upper - x[i]) / (-dir) + if x[i] ≥ upper - 1e-12 return zero(gamma_max) end - gamma_max = min(gamma_max, (lower - x[i]) / direction[i]) + gamma_max = min(gamma_max, (upper - x[i]) / (-direction[i])) end end return gamma_max From eed6005ba638e0225541c540c3dc14891dd58cc1 Mon Sep 17 00:00:00 2001 From: Deborah Hendrych Date: Sun, 29 Mar 2026 13:01:55 +0200 Subject: [PATCH 8/8] Some corrections. --- src/spanning_tree.jl | 28 ++++++++++++++++++++-------- 1 file changed, 20 insertions(+), 8 deletions(-) diff --git a/src/spanning_tree.jl b/src/spanning_tree.jl index fdfba4c..42e8874 100644 --- a/src/spanning_tree.jl +++ b/src/spanning_tree.jl @@ -303,10 +303,13 @@ function Boscia.bounded_compute_inface_extreme_point( N = length(direction) lb2 = copy(lb) ub2 = copy(ub) + # In-face fixings: iterates are floating; using machine `eps()` is too strict and + # can prevent DICG from identifying the correct minimal face. + tol = 1e-8 for i in 1:N - if x[i] ≤ eps() + if x[i] ≤ tol ub2[i] = 0.0 - elseif x[i] ≥ 1 - eps() + elseif x[i] ≥ 1 - tol lb2[i] = 1.0 end end @@ -324,10 +327,11 @@ function Boscia.is_simple_inface_feasible( ub, int_vars, ) + tol = 1e-8 for i in eachindex(a) - if x[i] ≤ eps() && a[i] > 1e-6 + if x[i] ≤ tol && a[i] > 1e-6 return false - elseif x[i] ≥ 1 - eps() && a[i] < 1 - 1e-6 + elseif x[i] ≥ 1 - tol && a[i] < 1 - 1e-6 return false end end @@ -348,6 +352,10 @@ function Boscia.bounded_dicg_maximum_step( ) T = promote_type(eltype(x), eltype(direction)) gamma_max = one(T) + # Safety margin: keep x - gamma*direction strictly inside bounds + # to prevent FrankWolfe's domain oracle from repeatedly rejecting points + # due to floating-point noise during line search. + tol = T(1e-12) for i in eachindex(x) if direction[i] == 0 continue @@ -359,17 +367,21 @@ function Boscia.bounded_dicg_maximum_step( if direction[i] > 0 # x_new decreases with gamma; only the lower bound can become active: # lower <= x[i] - gamma*dir => gamma <= (x[i] - lower)/dir - if x[i] ≤ lower + 1e-12 + if x[i] ≤ lower + tol return zero(gamma_max) end - gamma_max = min(gamma_max, (x[i] - lower) / direction[i]) + num = x[i] - lower - tol + num <= 0 && return zero(gamma_max) + gamma_max = min(gamma_max, num / direction[i]) else # direction[i] < 0: x_new increases with gamma; only the upper bound can become active: # x[i] - gamma*dir <= upper => gamma <= (upper - x[i]) / (-dir) - if x[i] ≥ upper - 1e-12 + if x[i] ≥ upper - tol return zero(gamma_max) end - gamma_max = min(gamma_max, (upper - x[i]) / (-direction[i])) + num = upper - x[i] - tol + num <= 0 && return zero(gamma_max) + gamma_max = min(gamma_max, num / (-direction[i])) end end return gamma_max