Relax primal-dual system upon linesearch failure

examples/Project.toml

@@ -1,19 +1,13 @@
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-FastDifferentiation = "eb9bf01b-bf85-4b60-bf87-ee5de06c00be"
-FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 MixedComplementarityProblems = "6c9e26cb-9263-41b8-a6c6-f4ca104ccdcd"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
-Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 TrajectoryGamesBase = "ac1ac542-73eb-4349-ae1b-660ab3609574"
 TrajectoryGamesExamples = "ff3fa34c-8d8f-519c-b5bc-31760c52507a"
-Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [sources]
-MixedComplementarityProblems = {path = ".."}
+MixedComplementarityProblems = {path = ".."}

src/solver.jl

@@ -22,7 +22,8 @@ function solve(
     x₀ = zeros(mcp.unconstrained_dimension),
     y₀ = ones(mcp.constrained_dimension),
     tol = 1e-4,
-    max_inner_iters = 20
+    max_inner_iters = 20,
+    max_outer_iters = 50
 )
     x = x₀
     y = y₀
@@ -31,11 +32,12 @@ function solve(
     ϵ = 1.0
     kkt_error = Inf
     status = :solved
-    while kkt_error > tol && ϵ > tol
-        iters = 1
+    outer_iters = 1
+    while kkt_error > tol && ϵ > tol && outer_iters < max_outer_iters
+        inner_iters = 1
         status = :solved
 
-        while kkt_error > ϵ && iters < max_inner_iters
+        while kkt_error > ϵ && inner_iters < max_inner_iters
             # Compute the Newton step.
             # TODO! Can add some adaptive regularization.
             F = mcp.F(x, y, s; θ, ϵ)
@@ -63,10 +65,15 @@ function solve(
             y += α_y * δy
 
             kkt_error = maximum(abs.(F))
-            iters += 1
+            inner_iters += 1
         end
 
-        ϵ *= 1 - exp(-iters)
+        ϵ *= (status == :solved) ? 1 - exp(-inner_iters) : 1 + exp(-inner_iters)
+        outer_iters += 1
+    end
+
+    if outer_iters == max_outer_iters
+        status = :failed
     end
 
     (; status, x, y, s, kkt_error, ϵ)

test/runtests.jl

@@ -16,7 +16,7 @@ using FiniteDiff: FiniteDiff
     M = [2 1; 1 2]
     A = [1 0; 0 1]
     b = [1; 1]
-    θ = zeros(2)
+    θ = rand(2)
 
     G(x, y; θ) = M * x - A' * y - θ
     H(x, y; θ) = A * x - b
@@ -28,12 +28,12 @@ using FiniteDiff: FiniteDiff
     end
 
     function check_solution(sol)
-        @test all(abs.(G(sol.x, sol.y; θ)) .≤ 1e-3)
+        @test all(abs.(G(sol.x, sol.y; θ)) .≤ 2e-3)
         @test all(H(sol.x, sol.y; θ) .≥ 0)
         @test all(sol.y .≥ 0)
-        @test sum(sol.y .* H(sol.x, sol.y; θ)) ≤ 1e-3
-        @test all(sol.s .≤ 1e-3)
-        @test sol.kkt_error ≤ 1e-3
+        @test sum(sol.y .* H(sol.x, sol.y; θ)) ≤ 2e-3
+        @test all(sol.s .≤ 2e-3)
+        @test sol.kkt_error ≤ 2e-3
         @test sol.status == :solved
     end