refactoring problem to avoid needing to construct separate mcps, by p…

…arameterizing M, A, b

benchmark/path.jl

@@ -1,16 +1,35 @@
 """ Generate a random large (convex) quadratic problem of the form
-                               min_x 0.5 xᵀ M x - θᵀ x
+                               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,
-    sparsity_rate = 0.9,
-)
+function generate_test_problem(; num_primals, num_inequalities)
+    G(x, y; θ) =
+        let
+            (; M, A, ϕ) = unpack_parameters(θ; num_primals, num_inequalities)
+            M * x - ϕ - A' * y
+        end
+
+    H(x, y; θ) =
+        let
+            (; A, b) = unpack_parameters(θ; num_primals, num_inequalities)
+            A * x - b
+        end
+
+    K(z, θ) =
+        let
+            x = z[1:num_primals]
+            y = z[(num_primals + 1):end]
+
+            [G(x, y; θ); H(x, y; θ)]
+        end
+
+    (; G, H, K)
+end
+
+"Generate a random parameter vector Θ corresponding to a convex QP."
+function generate_random_parameter(rng; num_primals, num_inequalities, sparsity_rate)
     bernoulli = Distributions.Bernoulli(1 - sparsity_rate)
 
     M = let
@@ -24,85 +43,92 @@ function generate_test_problem(
         randn(rng, num_inequalities, num_primals) .*
         rand(rng, bernoulli, num_inequalities, num_primals)
     b = randn(rng, num_inequalities)
+    ϕ = randn(rng, num_primals)
 
-    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]
+    [reshape(M, length(M)); reshape(A, length(A)); b; ϕ]
+end
 
-        [G(x, y; θ); H(x, y; θ)]
-    end
+"Unpack a parameter vector θ into the components of a convex QP."
+function unpack_parameters(θ; num_primals, num_inequalities)
+    M = reshape(θ[1:(num_primals^2)], num_primals, num_primals)
+    A = reshape(
+        θ[(num_primals^2 + 1):(num_primals^2 + num_inequalities * num_primals)],
+        num_inequalities,
+        num_primals,
+    )
 
-    (; G, H, K)
+    b =
+        θ[(num_primals^2 + num_inequalities * num_primals + 1):(num_primals^2 + num_inequalities * (num_primals + 1))]
+    ϕ = θ[(num_primals^2 + num_inequalities * (num_primals + 1) + 1):end]
+
+    (; M, A, b, ϕ)
 end
 
 "Benchmark interior point solver against PATH on a bunch of random QPs."
 function benchmark(;
-    num_problems = 10,
-    num_samples_per_problem = 100,
+    num_samples = 1000,
     num_primals = 100,
     num_inequalities = 100,
     sparsity_rate = 0.9,
 )
     rng = Random.MersenneTwister(1)
 
-    # Generate random problems and parameters.
-    problems = @showprogress desc = "Generating test problems..." map(1:num_problems) do _
-        generate_test_problem(rng; num_primals, num_inequalities, sparsity_rate)
-    end
+    # Generate problem and random parameters.
+    @info "Generating random problems..."
+    problem = generate_test_problem(; num_primals, num_inequalities)
 
-    θs = map(1:num_samples_per_problem) do _
-        randn(rng, num_primals)
+    θs = map(1:num_samples) do _
+        generate_random_parameter(rng; num_primals, num_inequalities, sparsity_rate)
     end
 
     # Generate corresponding MCPs.
-    ip_mcps = @showprogress desc = "Generating 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
-
-    path_mcps = @showprogress desc = "Generating 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
-
-    ip = @showprogress desc = "Solving IP MCPs..." map(ip_mcps) do mcp
-        # Make sure everything is compiled in a dry run.
-        MixedComplementarityProblems.solve(
+    @info "Generating IP MCP..."
+    parameter_dimension = length(first(θs))
+    ip_mcp = MixedComplementarityProblems.PrimalDualMCP(
+        problem.G,
+        problem.H;
+        unconstrained_dimension = num_primals,
+        constrained_dimension = num_inequalities,
+        parameter_dimension,
+    )
+
+    @info "Generating PATH MCP..."
+    lower_bounds = [fill(-Inf, num_primals); fill(0, num_inequalities)]
+    upper_bounds = fill(Inf, num_primals + num_inequalities)
+    path_mcp = ParametricMCPs.ParametricMCP(
+        problem.K,
+        lower_bounds,
+        upper_bounds,
+        parameter_dimension,
+    )
+
+    # Warm up the solvers.
+    @info "Warming up IP solver..."
+    MixedComplementarityProblems.solve(
+        MixedComplementarityProblems.InteriorPoint(),
+        ip_mcp,
+        first(θs),
+    )
+
+    @info "Warming up PATH solver..."
+    ParametricMCPs.solve(path_mcp, first(θs))
+
+    # Solve and time.
+    ip = @showprogress desc = "Solving IP MCPs..." map(θs) do θ
+        elapsed_time = @elapsed sol = MixedComplementarityProblems.solve(
             MixedComplementarityProblems.InteriorPoint(),
-            mcp,
-            zeros(num_primals),
+            ip_mcp,
+            θ,
         )
 
-        # Solve and time.
-        map(θs) do θ
-            elapsed_time = @elapsed sol = MixedComplementarityProblems.solve(
-                MixedComplementarityProblems.InteriorPoint(),
-                mcp,
-                θ,
-            )
-
-            (; elapsed_time, success = sol.status == :solved)
-        end
+        (; elapsed_time, success = sol.status == :solved)
     end
 
-    path = @showprogress desc = "Solving PATH MCPs..." map(path_mcps) do mcp
-        # Make sure everything is compiled in a dry run.
-        ParametricMCPs.solve(mcp, zeros(num_primals))
-
+    path = @showprogress desc = "Solving PATH MCPs..." map(θs) do θ
         # Solve and time.
-        map(θs) do θ
-            elapsed_time = @elapsed sol = ParametricMCPs.solve(mcp, θ)
+        elapsed_time = @elapsed sol = ParametricMCPs.solve(path_mcp, θ)
 
-            (; elapsed_time, success = sol.status == PATHSolver.MCP_Solved)
-        end
+        (; elapsed_time, success = sol.status == PATHSolver.MCP_Solved)
     end
 
     (; ip, path)
@@ -119,25 +145,15 @@ end
 
 "Estimate mean and standard deviation of runtimes for all problems."
 function runtime_stats(solver_data)
-    stats = map(solver_data) do problem_data
-        filtered_times = map(
-            datum -> datum.elapsed_time,
-            filter(datum -> datum.success, problem_data),
-        )
-        μ = Statistics.mean(filtered_times)
-        σ = Statistics.stdm(filtered_times, μ)
+    filtered_times =
+        map(datum -> datum.elapsed_time, filter(datum -> datum.success, solver_data))
+    μ = Statistics.mean(filtered_times)
+    σ = Statistics.stdm(filtered_times, μ)
 
-        (; μ, σ)
-    end
-
-    μ = map(datum -> datum.μ, stats)
-    σ = map(datum -> datum.σ, stats)
-    (; μ, σ, mean_μ = Statistics.mean(μ), mean_σ = Statistics.mean(σ))
+    (; μ, σ)
 end
 
 "Compute fraction of problems solved."
 function fraction_solved(solver_data)
-    Statistics.mean(solver_data) do problem_data
-        Statistics.mean(datum -> datum.success, problem_data)
-    end
+    Statistics.mean(datum -> datum.success, solver_data)
 end