bounded LMO version for the matching polytope (#24)

This PR adds the simple BLMO functions for the matching polytope. The overhead cost of creating the perturbed direction doesn't seem to be high so this can stay the default, as opposed to creating a fully-fledged BLMO with bound management

Author: matbesancon

src/matchings.jl

@@ -61,6 +61,59 @@ function FrankWolfe.compute_extreme_point(
     return v
 end
 
+function Boscia.bounded_compute_extreme_point(lmo::MatchingLMO, direction, lb, ub, int_vars; kwargs...)
+    # any entry i fixed to zero -> use a positive direction, ensuring the edge is not taken
+    # any entry i fixed to one with neighbors (u, v) -> use positive direction for all other neighbors of u, of v
+    corrected_direction = copy(direction)
+    for (idx, edge) in enumerate(edges(lmo.original_graph))
+        if ub[idx] ≈ 0
+            @assert lb[idx] ≈ 0
+            corrected_direction[idx] = 1
+        elseif lb[idx] ≈ 1
+            @assert ub[idx] ≈ 1
+            (vtx1, vtx2) = Tuple(edge)
+            # negative cost ensures the edge is taken, since no neighbor will be in the matching
+            corrected_direction[idx] = -1
+            for (idx2, e2) in enumerate(edges(lmo.original_graph))
+                # check if e2 is adjacent to edge
+                # we want one of the two nodes to be the same
+                if xor(src(e2) in (vtx1, vtx2), dst(e2) in (vtx1, vtx2))
+                    # we should not have incompatible edges fixed to one
+                    @assert lb[idx2] ≈ 0 "incompatible edges $edge $e2"
+                    corrected_direction[idx2] = 1
+                end
+            end
+        end
+    end
+    v = FrankWolfe.compute_extreme_point(lmo, corrected_direction)
+    @debug begin
+        for idx in eachindex(direction)
+            if ub[idx] ≈ 0
+                @assert v[idx] ≈ 0
+            elseif lb[idx] ≈ 1
+                @assert v[idx] ≈ 1
+            end
+        end
+    end
+    return v
+end
+
+function Boscia.is_simple_linear_feasible(lmo::MatchingLMO, v)
+    vertex_matching = zeros(Graphs.nv(lmo.original_graph))
+    for (idx, edge) in enumerate(edges(lmo.original_graph))
+        if v[idx] ≈ 0
+            continue
+        end
+        vtx1, vtx2 = Tuple(edge)
+        vertex_matching[vtx1] += v[idx]
+        vertex_matching[vtx2] += v[idx]
+    end
+    # one vertex fractionally matched to multiple edges
+    if maximum(vertex_matching) >= 1 + 1e-4
+        return false
+    end
+    return true
+end
 
 """
 PerfectMatchingLMO{G}(g::Graphs)

test/runtests.jl

@@ -38,31 +38,81 @@ Random.seed!(StableRNG(42), 42)
 end
 
 @testset "Matching LMO" begin
-    N = 200
+    N = 50
     Random.seed!(9754)
     g = Graphs.complete_graph(N)
     iter = collect(Graphs.edges(g))
     M = length(iter)
     direction = randn(M)
     lmo = CO.MatchingLMO(g)
     v = FrankWolfe.compute_extreme_point(lmo, direction)
+    @test Boscia.is_simple_linear_feasible(lmo, v)
     adj_mat = spzeros(M, M)
-    for i in 1:M
-        adj_mat[src(iter[i]), dst(iter[i])] = direction[i]
+    for (i, edge) in enumerate(edges(g))
+        adj_mat[src(edge), dst(edge)] = direction[i]
     end
     match_result = GraphsMatching.maximum_weight_matching(g, HiGHS.Optimizer, -adj_mat)
     v_sol = spzeros(M)
     K = length(match_result.mate)
-    for i in 1:K
-        for j in 1:M
-            if (match_result.mate[i] == src(iter[j]) && dst(iter[j]) == i)
-                v_sol[j] = 1
+    for k in 1:K
+        for (i, edge) in enumerate(edges(g))
+            if (match_result.mate[k] == src(edge) && dst(edge) == k)
+                v_sol[i] = 1
             end
         end
     end
     @test v_sol == v
     v2 = FrankWolfe.compute_extreme_point(lmo, ones(M))
     @test norm(v2) == 0
+    @test v == Boscia.bounded_compute_extreme_point(lmo, direction, zeros(M), ones(M), 1:M)
+    @testset "Fix one entry to zero" begin
+        for one_idx in SparseArrays.nonzeroinds(v)
+            # upperbound one everywhere except one_idx fixed to zero
+            v_fixed1 = Boscia.bounded_compute_extreme_point(lmo, direction, zeros(M), (1:M) .!= one_idx, 1:M)
+            @test v_fixed1[one_idx] == 0
+            @test Boscia.is_simple_linear_feasible(lmo, v_fixed1)
+        end
+    end
+    @testset "Fix a single entry to one" begin
+        for idx in rand(1:M, 100)
+            # skip if entry already at one
+            if v[idx] == 1
+                continue
+            end
+            lb = (1:M) .== idx
+            ub = ones(M)
+            v_fixed2 = Boscia.bounded_compute_extreme_point(lmo, direction, lb, ub, 1:M)
+            @test v_fixed2[idx] == 1
+            @test Boscia.is_simple_linear_feasible(lmo, v_fixed2)
+        end
+    end
+    @testset "Fix two entries to one" begin
+        for (i1, e1) in enumerate(edges(g))
+            for (i2, e2) in enumerate(edges(g))
+                # lighter on computation
+                if i1 ÷ 2 + i2 ÷ 2 > 0
+                    continue
+                end
+                # non-adjacent edges
+                if isempty(intersect(Tuple(e1), Tuple(e2)))
+                    lb = zeros(M)
+                    ub = ones(M)
+                    lb[i1] = lb[i2] = 1
+                    v_fixed3 = Boscia.bounded_compute_extreme_point(lmo, direction, lb, ub, 1:M)
+                    @test v_fixed3[i1] == 1
+                    @test v_fixed3[i2] == 1
+                    @test Boscia.is_simple_linear_feasible(lmo, v_fixed3)
+                end
+            end
+        end
+    end
+    # non-matching vector
+    v_wrong = 1.0 * copy(v)
+    idx1 = findfirst(==(Edge(1, 2)), collect(edges(g)))
+    idx2 = findfirst(==(Edge(1, 3)), collect(edges(g)))
+    v_wrong[idx1] = 0.75
+    v_wrong[idx2] = 0.5
+    @test !Boscia.is_simple_linear_feasible(lmo, v_wrong)
 end