diff --git a/benchmark/path.jl b/benchmark/path.jl index 33d06c5..ba63434 100644 --- a/benchmark/path.jl +++ 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