Add support to query the reduced cost of a variable

[mtanneau]

Querying the reduced cost of a variable is very common in linear programming (and the corresponding function is in all linear solvers' APIs), however there is no straightforward way of doing this through MOI.
Currently, to access the reduced cost of a variable x, one has to manually query the dual value of every possible bound constraint on x (if any), and aggregate the results.
That's very cumbersome, given that this feature is supported by all linear solvers.

Proposal:
Add a variable attribute to query the reduced cost of a variable, e.g. ReducedCost, or (more verbose) VariableReducedCost. For now, this would only apply to linear programming solvers (see list below).

The syntax would look like:

MOI.get(model, MOI.ReducedCost(), x)

Solvers wrappers to be updated:

• Clp.jl
• CPLEX.jl
• GLPK.jl
• Gurobi.jl
• MosekTools.jl
• SCIP.jl
• Xpress.jl

If it makes sense, support could be extended to non-linear solvers as well.

[mlubin]

This was discussed early in the design of MOI. I am not in favor of a ReducedCost attribute because it's redundant with ConstraintDual on the bound constraints. Another attribute means additional code on all solver interfaces and another set of tests to maintain. Furthermore solvers like Mosek, ECOS, SCS, OSQP, and Ipopt naturally return separate duals for the upper and lower bound constraints.

I would be in favor of having a helper function in MOI.Utilities that does the "query the dual value of every possible bound constraint on x (if any), and aggregate the results" because I agree that it's painful for to have to roll one's own version of this function.

[blegat]

The utility in MOI.Utilities could be get_fallback(::ModelLike, ::MOI.ReducedCost, ::MOI.VariableRef) allowing solvers to choose between the fallback or calling the solver API. That would require us to add it to MOI.Test though.

[mtanneau]

Having a universal fallback + a way to call the solver's API sounds like a good compromise. I also understand that adding an attribute and tests means more code maintenance.

One more thing I'd like to mention here (more for future readers of this issue).
As Miles pointed out, reduced costs only make sense for LP or when using a simplex-like algorithm (whenever there's some complementary slackness condition actually). I am not aware of the reduced cost concept being used in non-linear programming.
Therefore, if one were to use reduced costs, I would expect they do so in an LP context, for which they can use Clp, GLPK, or any commercial LP solver. If one really needs to use, say, SCS or ECOS and access reduced costs, then, IMO, it should be implemented at the solver level first, rather than MOI.
Expecting this to be supported by any and all solvers sounds (to me) like putting an unnecessary toll on MOI.

[mtanneau]

If the MOIU fallback + possibility to call the solver API directly (if one wishes to) is acceptable to all, then I can volunteer a PR for this. Most likely some time this summer.
It's probably best to wait for 0.9 to be released first?

[mlubin]

A ReducedCost attribute with a get_fallback implementation is reasonable. That's a more uniform way to handle the issue than having a one-off function in MOIU that won't be transparently accessible from JuMP.

It's probably best to wait for 0.9 to be released first?

The release wouldn't block the get_fallback work, only the implementation in the solvers, so up to you when to start.

[mtanneau]

@chriscoey and I discussed it a bit further and, rather than ReducedCost (which is only defined for LP + simplex), the name VariableDual might make more sense. (would it?)

• For LP, VariableDual would refer to the traditional reduced cost of said variable.
• For conic problems, VariableDual of a variable x ∈ K would be the dual variable s ∈ K* of that conic constraint. This does raise the question of whether VariableDual for a single scalar variable would make sense when x is a vector. @chriscoey @juan-pablo-vielma what are your thoughts on this?
• I don't know how much sense it makes for general NLP ; I guess the default fallback would be to query the dual variable associated to the bound(s) constraint(s) on x.

The release wouldn't block the get_fallback work, only the implementation in the solvers, so up to you when to start.

Time is not a pressing matter here (at least not for me). I'd rather think this through than rush into things.

[blegat]

This is related to constrained variables: #710.
For a constrained variable, it's dual is the dual of the constrained index returned by add_constrained_variables. Otherwise, it is unclear if the dual is part of several SingleVariable and/or VectorOfVariables constraints

[mtanneau]

This is related to constrained variables: #710.
For a constrained variable, it's dual is the dual of the constrained index returned by add_constrained_variables. Otherwise, it is unclear if the dual is part of several SingleVariable and/or VectorOfVariables constraints

I'm not sure I fully understand what you meant. I could not find add_constrained_variables in the doc (is that in a PR that hasn't been merged yet?).

Math-wise, do we agree that, for a conic (vector) variable x ∈ K, the corresponding VariableDual would be the dual (vector) s ∈ K* w.r.t that conic constraint? I have not seen any specific name for s in Conic optimization literature, which is why I was asking for some conic people's opinion.

Code-wise, I agree, if VariableDual is defined for a conic variable, then there should be a non-ambiguous way to access that conic constraint. Because vectors are involved, that's likely to be a bit trickier.

For example, if I have a conic variable (x1, x2, x3) ∈ SOC, and denote (s1, s2, s3) ∈ SOC* the corresponding dual vector, then

MOI.get(model, MOI.VariableDual(), x1)

would return s1. Obviously, s1 does not make sense by itself, so one would have to extract all components of s before doing anything meaningful.

[blegat]

is that in a PR that hasn't been merged yet?

Yes, see #759

Math-wise, do we agree that, for a conic (vector) variable x ∈ K, the corresponding VariableDual would be the dual (vector) s ∈ K* w.r.t that conic constraint?

Yes but currently (i.e. without add_constrained_variables), there are no conic variables. Variables are created free and some constraints are added on them. With add_constrained_variables, the variables are conic and their conic constraint is returned.

[mtanneau]

Then it's probably best to wait until #759 is merged. Happy to keep the discussion open though.

[mlubin]

I'm unsure I understand the motivation for defining another way to access the duals on x ∈ K constraints. You can already access duals on conic variables by using ConstraintDual on the VectorOfVariables-in-Cone constraint. #759 returns both a collection of variables and the corresponding constraint index that lets you access duals.

This also seems quite far away from the original issue of "Why can't I access reduced costs easily?".

[henriquebecker91]

May I contribute to this in some way? I am using the following excerpt in my code:

function reduced_cost(var) :: Float64
  rc = 0.0 
  has_upper_bound(var) && (rc += shadow_price(UpperBoundRef(var)))
  has_lower_bound(var) && (rc += shadow_price(LowerBoundRef(var)))
  is_fixed(var) && (rc += shadow_price(FixRef(var)))
  !has_upper_bound(var) && !has_lower_bound(var) && !is_fixed(var) &&
    @warn "CAUTION: reduce_cost was called over a variable with no bounds"
  rc  
end

But I would prefer to have a 'official' method to call instead. While I need to use reduced_costs on my research I have not a deep understanding of simplex and the maths related to it. So I would have a greater peace of mind if I could call a method already named reduced_cost (or taking an attribute with this name) and that was examined by people that understand more about the topic.

[odow]

Agreed that rather than at the MOI level, this might make sense at the JuMP level. We have shadow_price, for example.

[henriquebecker91]

What you would suggest to change in my method to qualify it as a contribution? With some direction I can make a PR.

[odow]

Probably just the following plus docs and tests. There's a question whether we want to use shadow_price or dual. I would probably vote for shadow_price since this function is likely going to be used by people first encountering JuMP.

function reduced_cost(x::VariableRef)::Float64
    if !has_duals(owner_model(x))
        error("... something descriptive ...")
    end
    if is_fixed(x)
        return shadow_price(FixRef(x))
    end
    rc = 0.0
    if has_upper_bound(x)
        rc += shadow_price(UpperBoundRef(x)))
    end
    if has_lower_bound(x)
        rc += shadow_price(LowerBoundRef(x))
    end
    return rc  
end

[henriquebecker91]

Since the start of the week I have worked intermittently on the idea of a JuMP.reduced_cost method, and the more I think about the method the more I discover I am not sure of what reduced cost means and/or how it should be implemented. I need someone to bounce some ideas with, my difficulties follow:

0. What is the scope of the concept of reduced cost? It seems to be mainly discussed (and their values of some relevance) in the context of optimally solved LP models. The method for now only checks if duals are available, should it check if the model is a LP (and not a MILP or ILP), is this check even possible? (M)ILP optimally solved models sometimes have the duals available (so this method sometimes would work on them). The concept seems to be relevant inside the solving of the dual simplex (before it becomes optimal), but I think it is not this version of the concept that we mean to expose (even because I am not sure it exists if something like the barrier algorithm is used to solve the LP under the hood).
1. The most common "intuitive definition" of reduced cost I found is along the lines of "reduced cost value indicates how much the objective function coefficient on the corresponding variable must be improved before the value of the variable will be positive in the optimal solution." The problem with this definition (and many others I have seen) is that there seems to be an implicit assumption that variables cannot ever assume negative values and (if we use shadow_price to implement it) they always have a non-negativity trivial constraint, a lower bound constraint. This by itself raises several questions.
2. If the variable has no lower bound, it has a concept of reduced cost? Should the method throws if a variable with no lower bound is passed to it? The lower bound may not be explicit, there may be a constraint that acts as a lower bound, avoiding the variable goes below zero in the model, and so it seems like something that should have a reduced cost, but it has not a reduced cost because it is not "a lower bound constraint as recognized by JuMP". What if it has a lower bound, but it is like >= -5 and there is another constraint that forces it to be >= 0 and the variable has value zero because of it?
3. If the variable has a lower bound, but it is negative, we should consider it valid for the concept of reduced cost? The definition above says "how much the objective function coefficient on the corresponding variable must be improved before the value of the variable will be positive in the optimal solution", but only the shadow_price over a lower bound in the form x >= 0 gives this, a shadow_price over a constraint x >= -10 gives the "how much the objective function coefficient on the corresponding variable must be improved before the value of the variable will be higher than -10".
4. (Only relevant if we think that reduced_cost makes sense in the (M)ILP context.) Binary variables have implicit lower bounds on zero, which make them very compatible with the intuitive definition of reduced cost, however, in JuMP binary variables does not have "a lower bound constraint as recognized by JuMP" (i.e., has_lower_bound on a newly created binary variable returns false, and LowerBoundRef errors, unless explicit lower bounds are given). This makes the reduced_cost not work out-of-the-box with binary variables, what will probably be seen as confusing by the users. More than that, I did not even check if adding lower bounds to binary variables make these "theoretically redundant/unnecessary" lower bound constraints report the correct shadow price values. This is something I need to check.
5. If the variable has a lower bound, but it is positive, what we should report? If there is an optimal solution for a model with such lower bound, then the variable already has non-zero value, and as such "If the optimal value of a variable is positive (not zero), then the reduced cost is always zero.", however, computing the shadow price over such lower bound constraint will not obligatorily give zero, it can give: (a) zero if the value of the variable is above such lower bound plus one (i.e., incrementing the constant of the lower bound by one unity does not affect the objective function value); (b) absolutely any value, if the variable value is stuck at this positive lower bound, and the value of the variable coefficient in objective function needs to change (may by a large value) to make the variable "attractive enough" (not by force of the lower bound but by its value in the objective function alone) to have a higher positive value than it already has.
6. The code I and odow shared above is not correct. A very straightforward test follows:

using JuMP, GLPK

@variable(model, 0 <= a <= 10) 
@objective(model, Max, a)
optimize!(model)
@show value(a)
@show reduced_cost(a)
@show shadow_price(LowerBoundRef(a))
@show shadow_price(UpperBoundRef(a))

gives

value(a) = 10.0
reduced_cost(a) = 1.0
shadow_price(LowerBoundRef(a)) = 0.0
shadow_price(UpperBoundRef(a)) = 1.0

What makes complete sense, given the definition of shadow_price but shows the reduced_cost code is wrong. "If the optimal value of a variable is positive (not zero), then the reduced cost is always zero." The method should only use the lower bound constraint, not both. But it is always the lower bound, or it depends on the optimization sense? If the variable can only assume negative values should it be the upper bound?

[mtanneau]

What is the scope of the concept of reduced cost?

The following should answer your points 0, 2, 3, and 5.

AFAIK, it is only used in LP -more specifically, for simplex-type algorithms.
If your LP writes

min    c'x
s.t.   A x = b
       l ⩽ x ⩽ u

then the dual writes

max    b'y + l's - u'z
s.t.   A' y + s - z = c
       s, z ⩾ 0

Note that s, z are the dual multipliers associated to the bound constraints l ⩽ x and x ⩽ u, respectively.
Mathematically, the "reduced cost" of x is c̄ = c - A'y = s-z.
99% of what you'll read in text books refers to the standard form where the only bounds are x ⩾ 0.

It then follows from complementary slackness that, if a bound constraint is not active, then the corresponding multiplier is zero. In the extreme case where one of the bounds in infinite, the corresponding multiplier is also zero.

Finally, for the duality-related conventions, you may want to check the JuMP docs and MOI docs on duality.

[odow]

Our definition of "reduced cost" is the dual on the bound constraints, which is s - z in @mtanneau's example.

To answer your questions:

0. Yes, mainly LPs. But our definition will "work" for any solver that returns duals of the variable bounds
1. Ignore that definition. Use the one above.
2. The reduced cost is on the bounds as recognized by JuMP. If x >= 0, and there is a linear constraint 1.0 * x >= 2, the reduced cost is 0, because the bound is not binding.
3. See our definition.
4. Solvers don't return dual variables for problems with binary variables. You can check if the solver found a dual solution with has_duals.
5. See our definition.
6. The code is correct.

[henriquebecker91]

I am sorry, but I think I will step out of this PR, at least for the moment. I discovered that:

• The definition of reduced cost the author of the paper I am reproducing is not the same in use by JuMP, so I will not be able to use this method.
• I am ignorant of the mathematical definitions and intricacies of the term "reduced cost" out of the "textbook" well-behaved cases. Seems clear that people with more experience in JuMP development are better qualified for polishing this method.
• The documentation I wrote until now was wrong taking into account the definition in use, and I also realized that the way the JuMP do the tests (using a MockOptimizer) does not match the way I was writing the tests. So nothing I had done since that excerpt of code I posted here is of use for the implementation of the JuMP.reduced_cost.

When I have more time, and if the situation of this issue continue to be open, I will try to give it some effort again.

[odow]

You're so close!

Ignore the "textbook" definitions.

Here is the correct code and documentation.

"""
    reduced_cost(x::VariableRef)

Return the reduced cost associated with variable `x`.

Equivalent to querying the shadow price of the active variable bound 
(if one exists and is active).

See also: [`shadow_price`](@ref).
"""
function reduced_cost(x::VariableRef)::Float64
    if !has_duals(owner_model(x))
        error("Unable to query reduced cost of $(x) because duals are not available.")
    end
    if is_fixed(x)
        return shadow_price(FixRef(x))
    end
    rc = 0.0
    if has_upper_bound(x)
        rc += shadow_price(UpperBoundRef(x)))
    end
    if has_lower_bound(x)
        rc += shadow_price(LowerBoundRef(x))
    end
    return rc  
end

[mlubin]

We should only be defining reduced cost if we're confident that it's compatible with standard textbook uses.

@chkwon, as a textbook author do you have thoughts on this?

[odow]

So this is the "textbook" value and sign of the reduced cost. For example, in @henriquebecker91's example, Excel (that bastion of a textbook implementation) also returns +1.

It also matches the implementation in Gurobi, which explicitly call it a reduced cost.

I more meant that @henriquebecker91 can ignore "textbook" definitions of the reduced cost that assume x >= 0 only.

[mlubin]

If we match up with Gurobi and Excel for both minimization and maximization and lower and upper bounds, then we should all set.

[odow]

For some reason, shadow_price on the bounds led to some negatives in the wrong places. I didn't check in enough detail. This works though:

using JuMP, Gurobi, Test

function reduced_cost(x::VariableRef)::Float64
    model = owner_model(x)
    if !has_duals(model)
        error("Unable to query reduced cost of $(x) because duals are not available.")
    end
    sense = objective_sense(model) == MOI.MIN_SENSE ? 1.0 : -1.0
    if is_fixed(x)
        return sense * dual(FixRef(x))
    end
    rc = 0.0
    if has_upper_bound(x)
        rc += dual(UpperBoundRef(x))
    end
    if has_lower_bound(x)
        rc += dual(LowerBoundRef(x))
    end
    return sense * rc  
end

@testset "Reduced cost" begin
    @testset "Error" begin
        model = direct_model(Gurobi.Optimizer())
        set_silent(model)
        @variable(model, x >= 1)
        @objective(model, Min, x)
        @test_throws ErrorException reduced_cost(x)
    end
    @testset "Free variable" begin
        model = direct_model(Gurobi.Optimizer())
        set_silent(model)
        @variable(model, x == 1.0)
        @variable(model, y)
        @objective(model, Min, x)
        optimize!(model)
        @test reduced_cost(y) == 0.0 ==
            Gurobi.get_dblattrelement(backend(model).inner, "RC", 2)
    end
    @testset "Fixed Variables" begin
        model = direct_model(Gurobi.Optimizer())
        set_silent(model)
        @variable(model, x == 10)
        @objective(model, Min, x)
        optimize!(model)
        @test reduced_cost(x) == 1.0 ==
            Gurobi.get_dblattrelement(backend(model).inner, "RC", 1)

        @objective(model, Max, x)
        optimize!(model)
        @test reduced_cost(x) == 1.0 ==
            Gurobi.get_dblattrelement(backend(model).inner, "RC", 1)

        @objective(model, Min, -x)
        optimize!(model)
        @test reduced_cost(x) == -1.0 ==
            Gurobi.get_dblattrelement(backend(model).inner, "RC", 1)

        @objective(model, Max, -x)
        optimize!(model)
        @test reduced_cost(x) == -1.0 ==
            Gurobi.get_dblattrelement(backend(model).inner, "RC", 1)
    end
    @testset "Bounds" begin
        model = direct_model(Gurobi.Optimizer())
        set_silent(model)
        @variable(model, 0 <= x <= 10)
        @objective(model, Min, x)
        optimize!(model)
        @test reduced_cost(x) == 1.0 ==
            Gurobi.get_dblattrelement(backend(model).inner, "RC", 1)

        @objective(model, Max, x)
        optimize!(model)
        @test reduced_cost(x) == 1.0 ==
            Gurobi.get_dblattrelement(backend(model).inner, "RC", 1)

        @objective(model, Min, -x)
        optimize!(model)
        @test reduced_cost(x) == -1.0 ==
            Gurobi.get_dblattrelement(backend(model).inner, "RC", 1)

        @objective(model, Max, -x)
        optimize!(model)
        @test reduced_cost(x) == -1.0 ==
            Gurobi.get_dblattrelement(backend(model).inner, "RC", 1)
    end
end

[henriquebecker91]

This is exactly the kind of thing I would not catch, but as you already provided the implementation and documentation I will try to whip up some tests, and make a PR. I am stuck at the fact that trying to do the same kind of test that the shadow_price has does not exactly work (the shadow_price tests do not need to access/"refer to" any variables), but it should not be too hard to overcome this problem. I am not sure how much I will work this weekend, but in the worst case I will try to have something by Monday.

[chkwon]

I think reduced_cost() will be a good addition. Especially for students to check homework answers :) In classrooms, it will be rare to consider "non-textbook" cases. Most of the time, it will be for nonnegative variables for minimization problems.

Reading this thread, I learned many new things about the meaning of the reduced costs. And I also have a new doubt.

Is it called 'reduced cost' for maximization problems also? What does it mean by 'cost'? For max c'x is it still c-z? Or, z-c? Or, after changing c to -c as in minimization and then -c-z? I think it can be quite confusing for maximization problems.

As far as I know, almost all graduate level textbooks are written for minimization problems. So no confusion. Interestingly, many undergraduate level textbooks are written for maximization problems. For example, Hillier & Lieberman is based on maximization and avoids the term 'reduced cost' at all. I see some books, written for maximization problems, just call c-z reduced cost at the same time. But it's not the best name, I think. Some use 'opportunity cost' or 'indicator'.

I'm not suggesting any new name for reduced_cost(), since it is called reduced cost. But I wonder if there is a consensus for the maximization case.

[odow]

Closed by jump-dev/JuMP.jl#2205