From 8f9cd837b4f81c74e17b31ddeca7e72377e8fa53 Mon Sep 17 00:00:00 2001 From: odow Date: Fri, 2 Oct 2020 09:16:41 +1300 Subject: [PATCH 1/3] Implement QCQP duals Co-authored-by: metab0t --- src/MOI/MOI_wrapper.jl | 88 +++++++++++++++++++++++++--- test/MathOptInterface/MOI_wrapper.jl | 57 +++++++----------- 2 files changed, 99 insertions(+), 46 deletions(-) diff --git a/src/MOI/MOI_wrapper.jl b/src/MOI/MOI_wrapper.jl index 135d73f6..67e71db0 100644 --- a/src/MOI/MOI_wrapper.jl +++ b/src/MOI/MOI_wrapper.jl @@ -2787,15 +2787,61 @@ function MOI.get( return _dual_multiplier(model) * p[] end -# function MOI.get( -# model::Optimizer, attr::MOI.ConstraintDual, -# c::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, <:Any} -# ) -# _throw_if_optimize_in_progress(model, attr) -# MOI.check_result_index_bounds(model, attr) -# pi = model.cached_solution.quadratic_dual[_info(model, c).row] -# return _dual_multiplier(model) * pi -# end +function MOI.get( + model::Optimizer, + attr::MOI.ConstraintDual, + c::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, <:Any}, +) + # For more information on QCP duals, see + # https://www.ibm.com/support/knowledgecenter/SSSA5P_12.10.0/ilog.odms.cplex.help/CPLEX/UsrMan/topics/cont_optim/qcp/17_QCP_duals.html + _throw_if_optimize_in_progress(model, attr) + MOI.check_result_index_bounds(model, attr) + # The derivative of a quadratic f(x) = x^TQx + a^Tx + b <= 0 is + # ∇f(x) = Q^Tx + Qx + a + # The dual is undefined if x is at the point of the cone. This can only be + # checked to numeric tolerances. We use `cone_top_tol`. + cone_top, cone_top_tol = true, 1e-6 + x = zeros(length(model.variable_info)) + ret = CPXgetx(model.env, model.lp, x, 0, length(x) - 1) + _check_ret(model, ret) + ∇f = zeros(length(x)) + a_i, a_v, qrow, qcol, qval = _CPXgetqconstr(model, c) + for (i, j, v) in zip(qrow, qcol, qval) + ∇f[i + 1] += v * x[j + 1] + ∇f[j + 1] += v * x[i + 1] + if abs(x[i + 1]) > cone_top_tol || abs(x[j + 1]) > cone_top_tol + cone_top = false + end + end + for (i, v) in zip(a_i, a_v) + ∇f[i + 1] += v + if abs(x[i + 1]) > cone_top_tol + cone_top = false + end + end + # TODO(odow): if at top of cone (x = 0) dual multiplier is ill-formed. + if cone_top + return NaN + end + qind = Cint(_info(model, c).row - 1) + nz_p, surplus_p = Ref{Cint}(), Ref{Cint}() + CPXgetqconstrdslack( + model.env, model.lp, qind, nz_p, C_NULL, C_NULL, 0, surplus_p + ) + ind = Vector{Cint}(undef, -surplus_p[]) + val = Vector{Cdouble}(undef, -surplus_p[]) + ret = CPXgetqconstrdslack( + model.env, model.lp, qind, nz_p, ind, val, -surplus_p[], surplus_p + ) + _check_ret(model, ret) + ∇f_max, ∇f_i = findmax(abs.(∇f)) + for (i, v) in zip(ind, val) + if i + 1 == ∇f_i + return _dual_multiplier(model) * (∇f_max > cone_top_tol ? v / ∇f[∇f_i] : 0.0) + end + end + return 0.0 +end function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue) _throw_if_optimize_in_progress(model, attr) @@ -3464,3 +3510,27 @@ function MOI.set( return end +function MOI.get( + model::Optimizer, + ::MOI.ConstraintDual, + c::MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.SecondOrderCone}, +) + f = MOI.get(model, MOI.ConstraintFunction(), c) + qind = Cint(_info(model, c).row - 1) + surplus_p = Ref{Cint}() + CPXgetqconstrdslack( + model.env, model.lp, qind, C_NULL, C_NULL, C_NULL, 0, surplus_p + ) + ind = Vector{Cint}(undef, -surplus_p[]) + val = Vector{Cdouble}(undef, -surplus_p[]) + ret = CPXgetqconstrdslack( + model.env, model.lp, qind, C_NULL, ind, val, -surplus_p[], surplus_p + ) + _check_ret(model, ret) + slack = zeros(length(model.variable_info)) + for (i, v) in zip(ind, val) + slack[i + 1] += v + end + indices = [_info(model, v).column for v in f.variables] + return _dual_multiplier(model) * slack[indices] +end diff --git a/test/MathOptInterface/MOI_wrapper.jl b/test/MathOptInterface/MOI_wrapper.jl index b15416ba..bbdf3fbc 100644 --- a/test/MathOptInterface/MOI_wrapper.jl +++ b/test/MathOptInterface/MOI_wrapper.jl @@ -12,14 +12,10 @@ const CONFIG = MOIT.TestConfig() const OPTIMIZER = CPLEX.Optimizer() MOI.set(OPTIMIZER, MOI.Silent(), true) -const BRIDGED_OPTIMIZER = MOI.Bridges.full_bridge_optimizer(OPTIMIZER, Float64) +MOI.set(OPTIMIZER, MOI.RawParameter("CPX_PARAM_REDUCE"), 0) +MOI.set(OPTIMIZER, MOI.RawParameter("CPX_PARAM_PRELINEAR"), 0) -const CERTIFICATE_OPTIMIZER = CPLEX.Optimizer() -MOI.set(CERTIFICATE_OPTIMIZER, MOI.Silent(), true) -MOI.set(CERTIFICATE_OPTIMIZER, MOI.RawParameter("CPX_PARAM_REDUCE"), 0) -MOI.set(CERTIFICATE_OPTIMIZER, MOI.RawParameter("CPX_PARAM_PRELINEAR"), 0) -const BRIDGED_CERTIFICATE_OPTIMIZER = - MOI.Bridges.full_bridge_optimizer(CERTIFICATE_OPTIMIZER, Float64) +const BRIDGED_OPTIMIZER = MOI.Bridges.full_bridge_optimizer(OPTIMIZER, Float64) function test_basic_constraint_tests() MOIT.basic_constraint_tests(BRIDGED_OPTIMIZER, CONFIG; exclude = [ @@ -62,12 +58,7 @@ function test_unittest() # TODO(odow): bug! We can't delete a vector of variables if one is in # a second order cone. "delete_soc_variables", - - # CPLEX returns INFEASIBLE_OR_UNBOUNDED without extra parameters. - # See below for the test. - "solve_unbounded_model", ]) - MOIT.solve_unbounded_model(CERTIFICATE_OPTIMIZER, CONFIG) end function test_modificationtest() @@ -76,18 +67,11 @@ end function test_contlineartest() MOIT.contlineartest(BRIDGED_OPTIMIZER, CONFIG, [ - # These tests require extra parameters to be set. - "linear8a", "linear8b", "linear8c", - # TODO(odow): This test requests the infeasibility certificate of a # variable bound. "linear12" ]) - MOIT.linear8atest(CERTIFICATE_OPTIMIZER, CONFIG) - MOIT.linear8btest(CERTIFICATE_OPTIMIZER, CONFIG) - MOIT.linear8ctest(CERTIFICATE_OPTIMIZER, CONFIG) - MOIT.linear12test(OPTIMIZER, MOIT.TestConfig(infeas_certificates=false)) end @@ -101,35 +85,34 @@ end function test_contquadratictest() MOIT.contquadratictest( - BRIDGED_CERTIFICATE_OPTIMIZER, - # TODO(odow): duals for quadratic problems. - MOIT.TestConfig(duals = false, atol = 1e-3, rtol = 1e-3), + BRIDGED_OPTIMIZER, + MOIT.TestConfig(atol = 1e-3, rtol = 1e-3), ["ncqcp"], # CPLEX doesn't support non-convex problems ) end function test_contconic() - MOIT.lintest(BRIDGED_OPTIMIZER, CONFIG, [ - # These tests require extra parameters to be set. - "lin3", "lin4" - ]) - - MOIT.lin3test(BRIDGED_CERTIFICATE_OPTIMIZER, CONFIG) - MOIT.lin4test(BRIDGED_CERTIFICATE_OPTIMIZER, CONFIG) - - # TODO(odow): duals for SOC constraints. - soc_config = MOIT.TestConfig(duals = false, atol=5e-3) + MOIT.lintest(BRIDGED_OPTIMIZER, CONFIG) - MOIT.soctest(BRIDGED_OPTIMIZER, soc_config, [ - "soc3" - ]) + soc_config = MOIT.TestConfig(atol=5e-3) + # TODO(odow): investigate why infeasibility certificates not generated for + # SOC. + MOIT.soctest(BRIDGED_OPTIMIZER, soc_config, ["soc3"]) MOIT.soc3test( BRIDGED_OPTIMIZER, - MOIT.TestConfig(duals = false, infeas_certificates = false, atol = 1e-3) + MOIT.TestConfig(atol = 1e-3, infeas_certificates = false) ) - MOIT.rsoctest(BRIDGED_OPTIMIZER, soc_config) + MOIT.rsoctest(BRIDGED_OPTIMIZER, soc_config, ["rotatedsoc2"]) + MOIT.rotatedsoc2test( + BRIDGED_OPTIMIZER, + # Need for `duals = false` fixed by https://github.com/jump-dev/MathOptInterface.jl/pull/1171 + # Remove in a future MOI 0.9.18+ release. + MOIT.TestConfig( + atol = 1e-3, infeas_certificates = false, duals = false + ), + ) MOIT.geomeantest(BRIDGED_OPTIMIZER, soc_config) end From 1d1dda1a6fccf9d49c8236ce0867c91622ba525e Mon Sep 17 00:00:00 2001 From: odow Date: Sun, 4 Oct 2020 20:01:16 +1300 Subject: [PATCH 2/3] Explain parameter choice --- test/MathOptInterface/MOI_wrapper.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/MathOptInterface/MOI_wrapper.jl b/test/MathOptInterface/MOI_wrapper.jl index bbdf3fbc..3f8c5635 100644 --- a/test/MathOptInterface/MOI_wrapper.jl +++ b/test/MathOptInterface/MOI_wrapper.jl @@ -12,8 +12,9 @@ const CONFIG = MOIT.TestConfig() const OPTIMIZER = CPLEX.Optimizer() MOI.set(OPTIMIZER, MOI.Silent(), true) +# Turn off presolve reductions so CPLEX will generate infeasibility +# certificates. MOI.set(OPTIMIZER, MOI.RawParameter("CPX_PARAM_REDUCE"), 0) -MOI.set(OPTIMIZER, MOI.RawParameter("CPX_PARAM_PRELINEAR"), 0) const BRIDGED_OPTIMIZER = MOI.Bridges.full_bridge_optimizer(OPTIMIZER, Float64) From fb469929bacd120a7f6c037b2ed419f268e51802 Mon Sep 17 00:00:00 2001 From: odow Date: Mon, 5 Oct 2020 09:17:36 +1300 Subject: [PATCH 3/3] Simplify dual calculations --- src/MOI/MOI_wrapper.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/MOI/MOI_wrapper.jl b/src/MOI/MOI_wrapper.jl index 67e71db0..9ad4a2c0 100644 --- a/src/MOI/MOI_wrapper.jl +++ b/src/MOI/MOI_wrapper.jl @@ -2835,9 +2835,11 @@ function MOI.get( ) _check_ret(model, ret) ∇f_max, ∇f_i = findmax(abs.(∇f)) - for (i, v) in zip(ind, val) - if i + 1 == ∇f_i - return _dual_multiplier(model) * (∇f_max > cone_top_tol ? v / ∇f[∇f_i] : 0.0) + if ∇f_max > cone_top_tol + for (i, v) in zip(ind, val) + if i + 1 == ∇f_i + return _dual_multiplier(model) * v / ∇f[∇f_i] + end end end return 0.0 @@ -3531,6 +3533,6 @@ function MOI.get( for (i, v) in zip(ind, val) slack[i + 1] += v end - indices = [_info(model, v).column for v in f.variables] - return _dual_multiplier(model) * slack[indices] + z = _dual_multiplier(model) + return [z * slack[_info(model, v).column] for v in f.variables] end