first cut of benchmark loop

benchmark/Project.toml

@@ -2,6 +2,9 @@
 MixedComplementarityProblems = "6c9e26cb-9263-41b8-a6c6-f4ca104ccdcd"
 PATHSolver = "f5f7c340-0bb3-5c69-969a-41884d311d1b"
 ParametricMCPs = "9b992ff8-05bb-4ea1-b9d2-5ef72d82f7ad"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [sources]
 MixedComplementarityProblems = {path = ".."}

benchmark/SolverBenchmarks.jl

@@ -0,0 +1,13 @@
+"Module for benchmarking different solvers against one another."
+module SolverBenchmarks
+
+using MixedComplementarityProblems: MixedComplementarityProblems
+using ParametricMCPs: ParametricMCPs
+using Random: Random
+using Statistics: Statistics
+using PATHSolver: PATHSolver
+using ProgressMeter: @showprogress
+
+include("path.jl")
+
+end # module SolverBenchmarks

benchmark/path.jl

@@ -0,0 +1,101 @@
+""" Generate a random large (convex) quadratic problem of the form
+                               min_x 0.5 xᵀ M x - θᵀ x
+                               s.t.  Ax - b ≥ 0.
+using Base: parameter_upper_bound
+
+NOTE: the problem may not be feasible!
+"""
+function generate_test_problem(
+    rng = Random.MersenneTwister(1);
+    num_primals = 1000,
+    num_inequalities = 1000,
+)
+    M = let
+        P = randn(rng, num_primals, num_primals)
+        P' * P
+    end
+
+    A = randn(rng, num_inequalities, num_primals)
+    b = randn(rng, num_inequalities)
+
+    G(x, y; θ) = M * x - θ - A' * y
+    H(x, y; θ) = A * x - b
+    K(z, θ) = begin
+        x = z[1:size(M, 1)]
+        y = z[(size(M, 1) + 1):end]
+
+        [G(x, y; θ); H(x, y; θ)]
+    end
+
+    (; G, H, K)
+end
+
+"Benchmark interior point solver against PATH on a bunch of random QPs."
+function benchmark(;
+    num_problems = 100,
+    num_samples_per_problem = 10,
+    num_primals = 1000,
+    num_inequalities = 1000,
+)
+    rng = Random.MersenneTwister(1)
+
+    # Generate random problems and parameters.
+    @showprogress desc = "Generating test problems..." problems = map(1:num_problems) do _
+        problem = generate_test_problem(rng; num_primals, num_inequalities)
+    end
+
+    θs = map(1:num_samples_per_problem) do _
+        randn(rng, num_primals)
+    end
+
+    # Generate corresponding MCPs.
+    @showprogress desc = "Generating IP MCPs... " ip_mcps = map(problems) do p
+        MixedComplementarityProblems.PrimalDualMCP(
+            p.G,
+            p.H;
+            unconstrained_dimension = num_primals,
+            constrained_dimension = num_inequalities,
+            parameter_dimension = num_primals,
+        )
+    end
+
+    @showprogress desc = "Generating PATH MCPs..." path_mcps = map(problems) do p
+        lower_bounds = [fill(-Inf, num_primals); fill(0, num_inequalities)]
+        upper_bounds = fill(Inf, num_primals + num_inequalities)
+        ParametricMCPs.ParametricMCP(p.K, lower_bounds, upper_bounds, num_primals)
+    end
+
+    @showprogress desc = "Solving IP MCPs..." ip_times = map(ip_mcps) do mcp
+        # Make sure everything is compiled in a dry run.
+        MixedComplementarityProblems.solve(
+            MixedComplementarityProblems.InteriorPoint(),
+            mcp,
+            zeros(num_primals),
+        )
+
+        # Solve and time.
+        map(θs) do θ
+            elapsed_time = @elapsed sol = MixedComplementarityProblems.solve(
+                MixedComplementarityProblems.InteriorPoint(),
+                mcp,
+                θ,
+            )
+
+            (; elapsed_time, sol.status)
+        end
+    end
+
+    @showprogress desc = "Solving PATH MCPs..." path_times = map(path_mcps) do mcp
+        # Make sure everything is compiled in a dry run.
+        ParametricMCPs.solve(mcp, zeros(num_primals))
+
+        # Solve and time.
+        map(θs) do θ
+            elapsed_time = @elapsed sol = ParametricMCPs.solve(mcp, θ)
+
+            (; elapsed_time, sol.status)
+        end
+    end
+
+    (; ip_times, path_times)
+end

test/runtests.jl

@@ -18,7 +18,7 @@ using FiniteDiff: FiniteDiff
     b = [1; 1]
     θ = rand(2)
 
-    G(x, y; θ) = M * x - A' * y - θ
+    G(x, y; θ) = M * x - θ - A' * y
     H(x, y; θ) = A * x - b
     K(z; θ) = begin
         x = z[1:size(M, 1)]