From f51e2937466f9192fb48dc98e111391ce7176528 Mon Sep 17 00:00:00 2001
From: "Jannis H." <44263801+JannisHal@users.noreply.github.com>
Date: Sat, 11 Jan 2025 10:16:44 +0100
Subject: [PATCH] Bug fix - LP direct solve (#537)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

* Bug fix - LP direct solve

* remove allocation with refactor

* Add test for bug fix

---------

Co-authored-by: Mathieu Besançon <mathieu.besancon@gmail.com>
---
 src/active_set_quadratic_direct_solve.jl | 18 +++++++------
 test/as_quadratic_projection.jl          | 32 ++++++++++++++++++++++++
 2 files changed, 42 insertions(+), 8 deletions(-)

diff --git a/src/active_set_quadratic_direct_solve.jl b/src/active_set_quadratic_direct_solve.jl
index 894c7c63b..a9fecdff6 100644
--- 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
index 2d9933430..8815765ae 100644
--- 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
\ No newline at end of file