Bug fix - LP direct solve (#537)

* Bug fix - LP direct solve * remove allocation with refactor * Add test for bug fix --------- Co-authored-by: Mathieu Besançon <[email protected]>
ZIB-IOL · Jan 11, 2025 · f51e293 · f51e293
1 parent 842f9e0
commit f51e293
Show file tree

Hide file tree

Showing 2 changed files with 42 additions and 8 deletions.
diff --git a/src/active_set_quadratic_direct_solve.jl b/src/active_set_quadratic_direct_solve.jl
@@ -284,18 +284,20 @@ function solve_quadratic_activeset_lp!(
     end
     sum_of_variables = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, λ), 0.0)
     MOI.add_constraint(o, sum_of_variables, MOI.EqualTo(1.0))
-    # Vᵗ A V λ == -Vᵗ b
-    for atom in as.atoms
+    # Wᵗ A V λ == -Wᵗ b
+    # V has columns vi
+    # W has columns vi - v1
+    for i in 2:nv
         lhs = MOI.ScalarAffineFunction{Float64}([], 0.0)
         Base.sizehint!(lhs.terms, nv)
         # replaces direct sum because of MOI and MutableArithmetic slow sums
         for j in 1:nv
             push!(
                 lhs.terms,
-                _compute_quadratic_constraint_term(atom, as.A, as.atoms[j], λ[j]),
+                _compute_quadratic_constraint_term(as.atoms[i], as.atoms[1], as.A, as.atoms[j], λ[j]),
             )
         end
-        rhs = -dot(atom, as.b)
+        rhs =  dot(as.atoms[1], as.b) - dot(as.atoms[i], as.b)
         MOI.add_constraint(o, lhs, MOI.EqualTo{Float64}(rhs))
     end
     MOI.set(o, MOI.ObjectiveFunction{typeof(sum_of_variables)}(), sum_of_variables)
@@ -373,12 +375,12 @@ function _compute_new_weights_wolfe_step(λ, ::Type{R}, old_weights, o::MOI.Abst
     return indices_to_remove, new_weights
 end
 
-function _compute_quadratic_constraint_term(atom1, A::AbstractMatrix, atom2, λ)
-    return MOI.ScalarAffineTerm(fast_dot(atom1, A, atom2), λ)
+function _compute_quadratic_constraint_term(atom1, atom0, A::AbstractMatrix, atom2, λ)
+    return MOI.ScalarAffineTerm(fast_dot(atom1, A, atom2) - fast_dot(atom0, A, atom2), λ)
 end
 
-function _compute_quadratic_constraint_term(atom1, A::Union{Identity,LinearAlgebra.UniformScaling}, atom2, λ)
-    return MOI.ScalarAffineTerm(A.λ * fast_dot(atom1, atom2), λ)
+function _compute_quadratic_constraint_term(atom1, atom0, A::Union{Identity,LinearAlgebra.UniformScaling}, atom2, λ)
+    return MOI.ScalarAffineTerm(A.λ * (fast_dot(atom1, atom2) - fast_dot(atom0, atom2)), λ)
 end
 
 struct LogScheduler{T}

diff --git a/test/as_quadratic_projection.jl b/test/as_quadratic_projection.jl
@@ -113,3 +113,35 @@ for traj in traj_data
     @test traj[end][2] ≤ 1e-8
     @test traj[end][4] ≤ 1e-7
 end
+
+lmo = FrankWolfe.LpNormLMO{Float64,5}(100.0)
+
+active_set_quadratic_manual_wolfe = FrankWolfe.ActiveSetQuadraticLinearSolve(
+    FrankWolfe.ActiveSet([(1.0, copy(x00))]),
+    2.0 * I, -2xp,
+    MOI.instantiate(MOI.OptimizerWithAttributes(HiGHS.Optimizer, MOI.Silent() => true)),
+    scheduler=FrankWolfe.LogScheduler(start_time=10, scaling_factor=1),
+    wolfe_step=true,
+)
+
+trajectoryBPCG_quadratic_manual = []
+x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient(
+    f,
+    grad!,
+    lmo,
+    active_set_quadratic_manual,
+    max_iteration=k,
+    callback=build_callback(trajectoryBPCG_quadratic_manual),
+);
+trajectoryBPCG_quadratic_manual_wolfe = []
+x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient(
+    f,
+    grad!,
+    lmo,
+    active_set_quadratic_manual_wolfe,
+    max_iteration=k,
+    callback=build_callback(trajectoryBPCG_quadratic_manual),
+);
+
+@test length(trajectoryBPCG_quadratic_manual) < 450
+@test length(trajectoryBPCG_quadratic_manual_wolfe) < 450