[breaking] Fix a<=b to return a-b in Nonpositives instead of b-a in Nonnegatives

docs/src/changelog.md

@@ -22,6 +22,11 @@ changes.
       elementwise values. This does not affect the result of `relative_entropy`.
     * The function `constant` should be used instead of the type `Constant`
       (which now refers to exclusively real constants).
+    * The constraint `a <= b` now produces `a - b in Nonpositives()` instead of
+      `b - a in Nonnegatives()`. The primal solutions are equivalent, but **the
+      dual variable associated with such constraints is now reversed in sign**.
+      (Following the convention in MathOptInterface, the dual of `a <= b` is
+      always negative, regardless of optimization sense.) (#593)
  * The syntaxes `dot(*)`, `dot(/)` and `dot(^)` have been removed in favor of
    explicit broadcasting (`x .* y`, `x ./ y`, and `x .^ y`). These were (mild)
    type piracy. In addition, `vecdot(x,y)` has been removed. Call

src/MOI_wrapper.jl

@@ -310,19 +310,13 @@ function MOI.get(
     return MOI.get(model.context.model, attr, ci)
 end
 
-# See `constraints/constraints.jl`
-# `LessThan` constraints are reformulated as `rhs - lhs` unlike MOI while
-# `GreaterThan` constraints are reformulated as `lhs - rhs` like in MOI
-_flip_dual(x, ::Type{S}) where {S<:MOI.LessThan} = -x
-_flip_dual(x, ::Type{S}) where {S<:MOI.AbstractScalarSet} = x
-
 function MOI.get(
     model::Optimizer,
     attr::Union{MOI.ConstraintDual,MOI.ConstraintPrimal},
     ci::MOI.ConstraintIndex{MOI.ScalarNonlinearFunction,S},
 ) where {S<:MOI.AbstractScalarSet}
     ret = MOI.get(model.context.model, attr, model.constraint_map[ci.value])
-    return _flip_dual(ret[], S)
+    return ret[]
 end
 
 function MOI.get(model::Optimizer, I::Type{<:MOI.Index}, name::String)

src/atoms/exp_cone/LogSumExpAtom.jl

@@ -32,7 +32,7 @@ function new_conic_form!(context::Context, e::LogSumExpAtom)
     # log(sum(exp(x))) <= t  <=>  sum(exp(x)) <= exp(t) <=> sum(exp(x - t)) <= 1
     t = Variable()
     z = sum(exp(e.children[1] - t * ones(size(e.children[1]))))
-    add_constraint!(context, z <= 1)
+    add_constraint!(context, 1 >= z)
     return conic_form!(context, t)
 end
 

src/atoms/lp_cone/DotSortAtom.jl

@@ -50,7 +50,7 @@ function new_conic_form!(context::Context{T}, x::DotSortAtom) where {T}
         y = vec(y)
     end
     μ, ν, e = Variable(size(y)), Variable(size(y)), ones(T, size(y))
-    add_constraint!(context, y * x.w' <= e * ν' + μ * e')
+    add_constraint!(context, e * ν' + μ * e' >= y * x.w')
     return conic_form!(context, sum(μ) + sum(ν))
 end
 

src/atoms/lp_cone/MinAtom.jl

@@ -47,7 +47,7 @@ function new_conic_form!(context::Context, x::MinAtom)
     t = Variable(x.size)
     t_obj = conic_form!(context, t)
     for child in x.children
-        add_constraint!(context, t <= child)
+        add_constraint!(context, child >= t)
     end
     return t_obj
 end

src/atoms/lp_cone/MinimumAtom.jl

@@ -24,7 +24,7 @@ evaluate(x::MinimumAtom) = minimum(evaluate(x.children[1]))
 
 function new_conic_form!(context::Context, x::MinimumAtom)
     t = Variable()
-    add_constraint!(context, t <= x.children[1])
+    add_constraint!(context, x.children[1] >= t)
     return conic_form!(context, t)
 end
 

src/atoms/lp_cone/SumLargestAtom.jl

@@ -40,7 +40,7 @@ function new_conic_form!(context::Context, x::SumLargestAtom)
     # sum k largest given by the solution to
     # minimize sum(t) + k*q
     # subject to c <= t + q, t >= 0
-    add_constraint!(context, c <= t + q)
+    add_constraint!(context, t + q >= c)
     add_constraint!(context, t >= 0)
     return conic_form!(context, sum(t) + x.k * q)
 end

src/constraints/LessThanConstraint.jl

@@ -40,16 +40,16 @@ function vexity(c::LessThanConstraint)
 end
 
 function _add_constraint!(context::Context{T}, lt::LessThanConstraint) where {T}
-    f = conic_form!(context, lt.rhs - lt.lhs)
+    f = conic_form!(context, lt.lhs - lt.rhs)
     if f isa AbstractVector
         # a trivial constraint without variables like `5 <= 0`
-        if !all(f .>= -CONSTANT_CONSTRAINT_TOL[])
+        if !all(f .<= CONSTANT_CONSTRAINT_TOL[])
             @warn "Constant constraint is violated"
             context.detected_infeasible_during_formulation[] = true
         end
         return
     end
-    set = MOI.Nonnegatives(MOI.output_dimension(f))
+    set = MOI.Nonpositives(MOI.output_dimension(f))
     context.constr_to_moi_inds[lt] = MOI_add_constraint(context.model, f, set)
     return
 end
@@ -61,10 +61,6 @@ Base.:<=(lhs::AbstractExpr, rhs::Value) = <=(lhs, constant(rhs))
 Base.:<=(lhs::Value, rhs::AbstractExpr) = <=(constant(lhs), rhs)
 
 function populate_dual!(model::MOI.ModelLike, c::LessThanConstraint, indices)
-    # FIXME(odow): this dual is the 'wrong' sign, because Convex implements
-    # rhs - lhs in Nonnegatives for a LessThanConstraint, instead of
-    # lhs - rhs in Nonpositives.
-    # Should we fix this?
     ret = MOI.get(model, MOI.ConstraintDual(), indices)
     c.dual = output(reshape(ret, c.size))
     return

src/problem_depot/problems/affine.jl

@@ -762,7 +762,7 @@ end
 
     handle_problem!(p)
     if test
-        @test p.constraints[1].dual ≈ 1 atol = atol rtol = rtol
+        @test p.constraints[1].dual ≈ -1 atol = atol rtol = rtol
     end
 
     x = Variable()

src/problem_depot/problems/lp.jl

@@ -15,7 +15,7 @@
     if test
         @test p.optval ≈ 1 atol = atol rtol = rtol
         @test evaluate(abs(x)) ≈ 1 atol = atol rtol = rtol
-        @test p.constraints[1].dual ≈ 1 atol = atol rtol = rtol
+        @test p.constraints[1].dual ≈ -1 atol = atol rtol = rtol
     end
 
     x = Variable(2, 2)

test/test_constraints.jl

@@ -170,9 +170,7 @@ function test_LessThanConstraint_dual_minimize()
     c = 1 <= 3 * x
     p = minimize(2.0 * x + 1.0, [c])
     solve!(p, SCS.Optimizer; silent_solver = true)
-    # This dual is currently positive despite a <= constraint because Convex
-    # implements rhs - lhs in R_+ instead of lhs - rhs in R_-.
-    @test isapprox(c.dual, 2 / 3; atol = 1e-4)
+    @test isapprox(c.dual, -2 / 3; atol = 1e-4)
     return
 end
 
@@ -181,9 +179,7 @@ function test_LessThanConstraint_dual_maximize()
     c = 3 * x <= 1
     p = maximize(2.0 * x + 1.0, [c])
     solve!(p, SCS.Optimizer; silent_solver = true)
-    # This dual is currently positive despite a <= constraint because Convex
-    # implements rhs - lhs in R_+ instead of lhs - rhs in R_-.
-    @test isapprox(c.dual, 2 / 3; atol = 1e-4)
+    @test isapprox(c.dual, -2 / 3; atol = 1e-4)
     return
 end