Approaching a Linear Program Correctly

Specific Domains Optimization (Mathematical)
, ,
1

Let’s see the Proposition 5.2.1 of the 2009 Bertsekas’s convex optimization theory:

(Linear Programming Duality Theorem)
(a) If any one of the 2 programs (primal and dual) has a finite optimal objective value, then the 2 programs’ optimal objective value coincide, moreover, attainment is ensured on both sides.

Actually, this theorem only concerns 1. the normal case (i.e. (a)), and 2. from unboundedness to infeasibility (i.e. (b) and (c), which I didn’t list here for conciseness).

Given this, in this post I’m addressing the missing abnormal circumstances (with step 1 and step 2 below).

Before starting, I mention this formula (★), which was proved in this post.
\sup \{y | y \ge 0\} = \inf \{ 0 | 0 \le x \le -1\} = + \infty

Let’s start here:
Suppose now we are dealing with an LP (■), with decision z, as follows
\sup \{ c^\top z | a^\top z \le b \}

The 1st step is:
Focus only on its feasibility system, namely \{z | a^\top z \le b\}. We ask “is it feasible”?
Well, if it is infeasible, then a^\top z \le b can be transformed (with variable substitution / elimination) to 0 \le x \le -1. The converse is also true.
If we have proved it is infeasible, then we quit this topic (no further steps needed). If we succeed in finding a feasible solution, then we continue with step 2.

The 2nd step is:
We ask “is it bounded”? Well, if (■) is unbounded, it should be equivalent to the leftmost “sup”-program in (★), whose dual feasibility system is infeasible (like 0 \le x \le -1). (Yes, this is where the name DUAL_INFEASIBLE is from). Note that we’ve done step 1 and in step 2 now, therefore DUAL_INFEASIBLE here really implies unbounded. (cf. Gurobi’s Optimization Status Codes, No. 5, which will politely suggests you to redo step 1 to avoid misinforming). If DUAL_INFEASIBLE is true, we quit this topic. If we succeed in finding a dual feasible solution, this indicates that (■) has a finite best dual bound, indicating that the (a) of [Bertsekas2009] above can be invoked. Then we continue with step 3.

The 3rd step:
Normality is ensured here, suggested by (a) of [Bertsekas2009] above.

2

How can we, users of Julia, help?

4 Likes
3

Hi @WalterMadelim, again it’s great to have people so actively posting in the forum :smile:

But I might ask that we convert your effort on tutorial-type math focused posts like this one and A naive explanation of infeasibility and unboundedness into improvements to the JuMP documentation?

Let’s try to keep the posts on the Julia forum focused on code problems with Julia and JuMP. For example, the bugs you’re finding in Gurobi.jl are excellent and very helpful!

6 Likes
4

If it contains a question, typically a “question” tag will be visible.
Otherwise typically it is about sharing experience to facilitate developing / referencing / etc., like this.

5 
help?> JuMP.DUAL_INFEASIBLE
  DUAL_INFEASIBLE::TerminationStatusCode

  An instance of the TerminationStatusCode enum.

  DUAL_INFEASIBLE: The algorithm concluded that no dual bound exists for the problem. If,
  additionally, a feasible (primal) solution is known to exist, this status typically implies that   
  the problem is unbounded, with some technical exceptions.

I wonder what is meant by saying “typically”, “some technical exceptions”?
I thought it should be deterministic, i.e.,
[ primal feasible + DUAL_INFEASIBLE ] directly implies [primal unbounded].

6

That is precisely why Oscar was suggesting contributions to the JuMP.jl documentation instead. While this forum is great for questions and interactive exchanges, the formal documentation is better suited for “documenting / referencing / etc.”

2 Likes
7

Okay, I’ll consider that when I’m confident.

1 Like
8

Oh, I seems to know why.

I here furnish an example (i.e., the data A, b, c), which renders the feasibility system on both the primal and dual side INFEASIBLE. (I have verified it using Gurobi and Mosek). The primal and dual LP formulation follows (2.12) and (2.13) in this page.

(A, b, c) = ([-15 10 0 11; -19 0 19 -17; 9 -16 15 -18], [-19, -13, 19], [-19, 0, 15, -12])

This is a perfect complement to the Theorem of [Bertsekas2009].
I’m happy on contributing this :slightly_smiling_face:.

9

The above 3-by-4 example was generated by Julia’s rand().
Here is an easier 2-by-3 example invented by my hand.
The runnable code is provided here

import JuMP, MosekTools, Ipopt
import Gurobi; const GRB_ENV = Gurobi.Env();
function optimise(m) return (JuMP.optimize!(m); (JuMP.termination_status(m), JuMP.primal_status(m))) end
function Model(solver, mode)
    solver == "Ipopt"  && (m = JuMP.Model(Ipopt.Optimizer))
    solver == "Mosek"  && (m = JuMP.Model(() -> MosekTools.Optimizer()))
    solver == "Gurobi" && (m = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)))
    occursin("m", mode) && JuMP.set_objective_sense(m, JuMP.MIN_SENSE)
    occursin("M", mode) && JuMP.set_objective_sense(m, JuMP.MAX_SENSE)
    occursin("s", mode) && JuMP.set_silent(m)
    m
end

A, b, c = [1 -1 0; 0 0 1], [0, -1], [0, -1, 1]

model = Model("Mosek", "s")
JuMP.@variable(model, x[eachindex(c)] >= 0)
JuMP.@constraint(model, A * x .== b)
t_primal = optimise(model)[1]

model = Model("Mosek", "s")
JuMP.@variable(model, y[eachindex(b)])
JuMP.@constraint(model, c .>= transpose(A) * y)
t_dual = optimise(model)[1]
println(t_primal, ", ", t_dual) # ✅ INFEASIBLE, INFEASIBLE
10

We have had a number of discussions over the years about the extent to which we want to make the JuMP documentation about JuMP, or whether we want to extend it to be more of a textbook for optimization more generally (e.g., Suggestions for documentation improvements · Issue #2348 · jump-dev/JuMP.jl · GitHub).

We’ve tended to stick fairly close to JuMP. So we’d regard the knowledge that an LP can be both primal and dual infeasible as theory knowledge that is more appropriate somewhere else. It’s up for debate though. This information could be useful in the JuMP documentation.

11

If one inspects the No.5 termination_status_code UNBOUNDED of Gurobi’s doc.

Model was proven to be unbounded. Important note: an unbounded status indicates the presence of an unbounded ray that allows the objective to improve without limit. It says nothing about whether the model has a feasible solution. If you require information on feasibility, you should set the objective to zero and reoptimize.

He might feel a bit weird. Since:
From a rigorous point of view, if we are handling a model.
The first and foremost thing we should make sure is “is it feasible?”
The secondary thing is then asking “is it unbounded?”
Therefore, strictly speaking, claiming [model is “unbounded”] before assuring [ model has a feasible solution ] is meaningless.
Therefore, the “unbounded ray” which Gurobi claims to exists is also rendered meaningless. (Indeed, only the ray who starts at a feasible point and extends to infinity without hitting the infeasible region is meaningful and desirable in practice).

By comparison, JuMP’s explanation for DUAL_INFEASIBLE is less disquieting:
Its second statement is flawless: If, additionally, a feasible (primal) solution is known to exist, this status typically implies that the problem is unbounded, with some technical exceptions. where the exceptions can be exemplified by my 2-by-3 case above.

But I can’t fully agree its first statement The algorithm concluded that no dual bound exists for the problem. , which I think should be substituted by The dual program is proved to be infeasible. (just as its name verbatim, my humble personal opinion, I don’t know how solvers perform behind the scenes).

Things only become unambiguous when model is an LP, within which:

  1. DUAL_INFEASIBLE + PRIMAL_FEASIBLE implies PRIMAL_UNBOUNDED
  2. The situation DUAL_INFEASIBLE + PRIMAL_INFEASIBLE does exist
12

The case that solvers report some equivalent to unbounded when they really mean “I found an unbounded ray but not necessarily a feasible point” is precisely why JuMP has the status DUAL_INFEASIBLE instead of UNBOUNDED.

The doc strings are here: MathOptInterface.jl/src/attributes.jl at 8f2f2539e00fd5039f1eccc015b48c69cc925613 · jump-dev/MathOptInterface.jl · GitHub

Feel free to open a PR with suggested improvements. I’ve opened on here: [docs] improve the docstrings of the various enum values by odow · Pull Request #2698 · jump-dev/MathOptInterface.jl · GitHub

1 Like
13

Yes, DUAL_INFEASIBLE alone does not implies PRIMAL_UNBOUNDED, even in the simplest LP setting.
It appears I’m clearer now. I’ll later update

14

In this post I expound the essential implication of the termination status DUAL_INFEASIBLE.

We need this Linear Farkas’ Lemma (Proposition 5.1.2, [Bertsekas 2009])

(LFL) The following 2 condition are equivalent

  1. +\infty = \sup_y \{ b^\top y | A^\top y \le 0 \}
  2. \{ x | x \ge 0, A x= b \} = \emptyset

Now suppose our primal problem is
\sup_y \{ b^\top y | A^\top y \le c \}
We routinely write its dual problem as (in accordance with this)
\inf_x \{ c^\top x | x \ge 0, A x= b \}

In view of (LFL), we conclude that DUAL_INFEASIBLE is essentially saying that

  • (:spade_suit:) the modified primal problem \sup_y \{ b^\top y | A^\top y \le 0 \} is UNBOUNDED.

We can see that the modified primal problem admits of a trivial solution y = 0.
But we cannot claim that the primal problem is feasible, unless we set its objective b^\top y to 0 and then hand it to a solver.

If a feasible solution of the primal problem is already assured, then at this stage can we immediately conclude that PRIMAL_UNBOUNDED (due to (:spade_suit:)).

15 
The algorithm concluded that no dual bound exists for the problem.
 
If the problem is a conic optimization problem (thus also a linear program), this status means the dual
problem is infeasible.

If a primal feasible solution exists, this status typically implies that the
problem is unbounded, with some technical exceptions (for example, if the
problem is a conic optimization problem in which strong duality does not
hold).

To check if the primal is unbounded, set the objective sense to
[`FEASIBILITY_SENSE`](@ref) and re-solve the problem. If a primal feasible
point exists, the original problem is unbounded. If a primal feasible point
does not exist, the original problem is both primal and dual infeasible.

Is this the revised docstring? It’s probably not correct either, is it?

The last (4th) segment, shouldn’t it be “To check if the primal is feasible, set …”.
If a primal feasible point exists, the original problem is unbounded. this sentence is also not correct (It is correct only if the primal is an LP).
These can be verified just based on the 2-by-2 SDP example you have listed in the Github.

Additionally, I think in the current status, the info is a bit dense. The 2nd and 3rd segment is not problematic, though. But I think the 1st statement can be removed now. I think the 1st statement is a bit aggressive, whose essential meaning can be embodied in my revised one below.
A nontrivial bound on the optimal objective value cannot be obtained through further exploring the dual problem.
in which trivial bound means: -Inf for a inf-program, +Inf for a sup-program

16

How about: [docs] improve the docstring for DUAL_INFEASIBLE by odow · Pull Request #2701 · jump-dev/MathOptInterface.jl · GitHub?

(Let’s keep future comments to the GitHub PR.)

1 Like