test runs but getting numerical issues

Project.toml

@@ -7,6 +7,8 @@ version = "0.1.0"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 FastDifferentiation = "eb9bf01b-bf85-4b60-bf87-ee5de06c00be"
+Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 TrajectoryGamesBase = "ac1ac542-73eb-4349-ae1b-660ab3609574"
@@ -15,6 +17,8 @@ TrajectoryGamesBase = "ac1ac542-73eb-4349-ae1b-660ab3609574"
 BlockArrays = "0.16.43"
 DataStructures = "0.18.20"
 FastDifferentiation = "0.4.2"
+Infiltrator = "1.8.3"
+LinearAlgebra = "1.11.0"
 SparseArrays = "1.11.0"
 Symbolics = "6.19.0"
 TrajectoryGamesBase = "0.3.10"

src/MCPSolver.jl

@@ -3,6 +3,9 @@ module MCPSolver
 using SparseArrays: SparseArrays
 using FastDifferentiation: FastDifferentiation as FD
 using Symbolics: Symbolics
+using LinearAlgebra: I, norm
+
+using Infiltrator
 
 include("SymbolicUtils.jl")
 include("sparse_utils.jl")

src/mcp.jl

@@ -51,7 +51,7 @@ function PrimalDualMCP(
 ) where {T<:Union{FD.Node,Symbolics.Num}}
     # Create symbolic slack variable `s` and parameter `ϵ`.
     s_symbolic = SymbolicUtils.make_variables(backend, :s, length(y_symbolic))
-    ϵ_symbolic = SymbolicUtils.make_variables(backend, :ϵ, 1)
+    ϵ_symbolic = only(SymbolicUtils.make_variables(backend, :ϵ, 1))
     z_symbolic = [x_symbolic; y_symbolic; s_symbolic]
 
     F_symbolic = [

src/solver.jl

@@ -1,6 +1,8 @@
 abstract type SolverType end
 struct InteriorPoint <: SolverType end
 
+using Infiltrator: @infiltrate
+
 """ Basic interior point solver, based on Nocedal & Wright, ch. 19.
 Computes step directions `δz` by solving the relaxed primal-dual system, i.e.
                          ∇F(z; ϵ) δz = -F(z; ϵ).
@@ -19,39 +21,49 @@ function solve(
     ::InteriorPoint,
     mcp::PrimalDualMCP;
     x₀ = zeros(mcp.unconstrained_dimension),
-    y₀ = zeros(mcp.constrained_dimension),
+    y₀ = ones(mcp.constrained_dimension),
     tol = 1e-4
 )
     x = x₀
     y = y₀
     s = ones(length(y))
 
-    ϵ = 1.0
+    ϵ = 10.0
     kkt_error = Inf
-    while kkt_error > tol
+    while kkt_error > tol && ϵ > tol
         iters = 1
         while kkt_error > ϵ
             # Compute the Newton step.
-            F = mcp.F(x, y, s, ϵ)
-            δz = -mcp.∇F(x, y, s, ϵ) \ F
+            F = mcp.F(x, y, s; ϵ)
+            δz = -(mcp.∇F(x, y, s; ϵ) + ϵ*I) \ F
 
             # Fraction to the boundary linesearch.
             δx = @view δz[1:mcp.unconstrained_dimension]
             δy = @view δz[(mcp.unconstrained_dimension + 1):(mcp.unconstrained_dimension + mcp.constrained_dimension)]
             δs = @view δz[(mcp.unconstrained_dimension + mcp.constrained_dimension + 1):end]
 
-            α_s = fraction_to_the_boundary_linesearch(s, δs)
-            α_y = fraction_to_the_boundary_linesearch(y, δy)
+            α_s = fraction_to_the_boundary_linesearch(s, δs; tol)
+            α_y = fraction_to_the_boundary_linesearch(y, δy; tol)
+
+            if isnan(α_s) || isnan(α_y)
+                @warn "Linesearch failed. Exiting prematurely."
+                break
+            end
 
             # Update variables accordingly.
             x += α_s * δx
             s += α_s * δs
             y += α_y * δy
 
-            kkt_error = norm(F)
+            kkt_error = maximum(abs.(F))
             iters += 1
+            @info iters, s
         end
 
+        @info x
+        @info y
+        @info s
+        @info ϵ
         ϵ *= 1 - exp(-iters)
     end
 
@@ -61,12 +73,11 @@ end
 """Helper function to compute the step size `α` which solves:
                    α* = max(α ∈ [0, 1] : v + α δ ≥ (1 - τ) v).
 """
-function fraction_to_the_boundary_linesearch(v, δ; τ = 0.995, decay = 0.5)
+function fraction_to_the_boundary_linesearch(v, δ; τ = 0.995, decay = 0.5, tol = 1e-4)
     α = 1.0
-    while any(v + α * δ .≥ (1 - τ) * v)
-        if α < 1e-2
-            @warn "Fraction to the boundary linesarch failed."
-            break
+    while any(v + α * δ .< (1 - τ) * v)
+        if α < tol
+            return NaN
         end
 
         α *= decay

test/runtests.jl

@@ -1,25 +1,28 @@
-using Test: @testset, @test
-using MCPSolver
-
-
 """ Test for the following QP:
                            min_x 0.5 xᵀ M x
                            s.t.  Ax - b ≥ 0.
 Taking `y ≥ 0` as a Lagrange multiplier, we obtain the KKT conditions:
                              G(x, y) = Mx - Aᵀy = 0
                              0 ≤ y ⟂ H(x, y) = Ax - b ≥ 0.
 """
+
+using Test: @testset, @test
+using MCPSolver
+
+
 @testset "QPTestProblem" begin
-    M = [2, 1; 1, 2]
-    A = [1, 0.3]
-    b = 1.1
+    M = [2 1; 1 2]
+    A = [1.0 0.0; 0.0 1.0]
+    b = [1.0; 1.0]
 
     G(x, y) = M * x - A' * y
     H(x, y) = A * x - b
 
-    mcp = MCPSolver.to_symbolic_mcp(G, H, 2, 1)
+    mcp = MCPSolver.to_symbolic_mcp(G, H, size(M, 1), length(b))
     sol = MCPSolver.solve(MCPSolver.InteriorPoint(), mcp)
 
+    println(sol)
+
     @test all(abs.(G(sol.x, sol.y)) .≤ 1e-3)
     @test all(H(sol.x, sol.y) .≥ 0)
     @test all(sol.y .≥ 0)