## Friday, July 30, 2010

### Infeasible LPs and Farkas Certificates

Yesterday I wrote about dealing with unbounded LPs in CPLEX. Today we'll look at infeasible LPs. CPLEX APIs provide a variety of mechanisms for dealing with infeasibility, including refineConflict()/getConflict() to isolate irreducible infeasible subsystems (IISes) and feasOpt() to find a “cheapest” relaxation of constraints and bounds that makes the model feasible. They also provide access to what is known as a Farkas certificate.

Thanks to Bo Jensen, who reminded me about Farkas certificates in a comment yesterday. He also provided a link to MOSEK's description of them, which saves me quite a bit of typing. Farkas certificates draw their name from Farkas's lemma, a well-known result in linear programming theory. While it is possible for both a primal linear program and its dual to be infeasible, that condition is rare in practice. The (much) more common case is that an infeasible LP has an unbounded dual. In essence, a Farkas certificate is a recession direction for the dual problem. It's existence proves the infeasibility of the primal.

For infeasible LPs, CPLEX provides the IloCplex.dualFarkas() method to get a Farkas certificate. (As always, I'll stick to the Java API, but all the APIs should provide some version of this.) The method has two arguments used for returning results (and not for providing inputs). One argument returns a list of primal constraints with non-zero dual values; the other vector holds the corresponding dual values.

Consider the following primal LP:\begin{eqnarray*} \mathrm{minimize} & u+v-2w\\ \mathrm{s.t.} & u-2v-w & \ge3\\ & -2u+v-w & \ge2\\ & u,v,w & \ge0\end{eqnarray*}whose dual just happens, by sheer coincidence, to be the LP from yesterday's post. Dual variables $x$ and $y$ correspond to the first and second constraints respectively. The dual problem, shown in Figure 1, is obviously unbounded, so the primal is infeasible. We can easily deduce infeasibility directly from the formulation: add the two constraints together to get $-u-v-2w\ge5$, which cannot be satisfied with nonnegative values of the variables.

Continuing yesterday's Java code (without repeating it), I'll start by constructing the new primal problem:
IloCplex primal = new IloCplex();
IloNumVar u = primal.numVar(0, Double.MAX_VALUE, "u");
IloNumVar v = primal.numVar(0, Double.MAX_VALUE, "v");
IloNumVar w = primal.numVar(0, Double.MAX_VALUE, "w");
IloNumVar[] vars = new IloNumVar[] {u, v, w};
primal.scalProd(new double[] {1.0, 1.0, -2.0}, vars));
IloConstraint c1 =
primal.scalProd(new double[] {1.0, -2.0, -1.0}, vars),
3.0);
IloConstraint c2 =
primal.scalProd(new double[] {-2.0, 1.0, -1.0}, vars),
2.0);
Note that I assign the constraints to program variables as I add them to the model, so that I can refer back to them later.

Since I cannot be sure which constraints dualFarkas will return, nor in what order, it will be handy to construct a map associating each constraint with the corresponding dual variable.
HashMap<IloConstraint, IloNumVar> map =
new HashMap<IloConstraint, IloNumVar>();
map.put(c1, x);
map.put(c2, y);
The dual variables x and y (instances of IloNumVar) were defined in yesterday's portion of the code.

Now I turn off the presolver, solve the model, and verify that it is infeasible. (As was the case yesterday with an unbounded primal, fail to turn off the presolver will cause CPLEX to throw an exception when I try to access the Farkas certificate.)  I also need to use the dual simplex algorithm; any other algorithm (primal simplex, barrier, ...) will cause CPLEX to throw an exception. There is no explicit statement in the code since CPLEX's default choice for algorithm on this problem is dual simplex, but in general it should probably be specified just in case CPLEX prefers some other algorithm.
primal.setParam(IloCplex.BooleanParam.PreInd, false);
primal.solve();
System.out.println("CPLEX status: "
+ primal.getCplexStatus());
The output here is CPLEX status: Infeasible.

The next step is to access the Farkas certificate:
IloRange[] farkasConstraints = new IloRange[2];
double[] farkasValues = new double[2];
primal.dualFarkas(farkasConstraints, farkasValues);
System.out.print("Farkas certificate:");
for (int i = 0; i < 2; i++) {
System.out.print("\t" + map.get(farkasConstraints[i])
+ " -> " + farkasValues[i]);
}
System.out.println("");
which results in Farkas certificate: x -> 0.6666666666666666 y -> 0.3333333333333333. So the Farkas certificate proves to be the direction vector of the extreme ray anchored at $A$ in Figure 1. Note the use of the previously constructed map to match the constraints returned in farkasConstraints with the corresponding dual values. Although both constraints have nonzero direction components in this example, in general not all constraints need be returned, nor can I be sure that those returned will be in the order I added them to the model.

To find the vertex at which it is anchored, I can access the last feasible dual solution obtained before CPLEX stopped.
double[] q = primal.getDuals(farkasConstraints);
System.out.print("Dual vertex:");
for (int i = 0; i < 2; i++) {
System.out.print("\t" + map.get(farkasConstraints[i])
+ " -> " + q[i]);
}
System.out.println("");
The output is as follows: Dual vertex: x -> 1.6666666666666665 y -> 0.33333333333333326. This is vertex $A$ in Figure 1. So, at least in this example, the Farkas certificate produces the same extreme ray that I got yesterday by solving today's dual (yesterday's primal) directly.

I still have one concern. The definition of a Farkas certificate only guarantees a recession direction for the dual; it does not need to be an extreme dual direction to qualify as a Farkas certificate. I suspect, but do not know with certainty, that the mechanism by which solvers select a Farkas certificate will consistently lead to an extreme direction.

1. The properties of the farkas ray returned, depends on the solver used. If you use an interior point optimizer (which allows you access to the interior point solution before basis identification), then the ray will not be an extreme ray if multiple solution exists. If the optimizer used is simplex, then the ray is an extreme ray. All the information needed is right there in the simplex method when it is either infeasible or unbounded, you just have to get it back through the postsolve. I think it only gets more confusing for the user if they have to switch off presolve and use certain algorithms etc. The optimizer should return farkas rays pr. default and report the quality of them, then there can be no discussion if the problem is infeasible or not, it is a proof. Farkas rays can be used in many cases such as DW, where you practically can implement phase one as phase two.

2. If you're using primal simplex on an unbounded problem or dual simplex on an infeasible problem, the Farkas certificate found should be an extreme direction. What if you're using dual simplex on an unbounded problem or primal simplex (phase I) on an infeasible one? Do you still have an obvious Farkas certificate in the last iteration, or do you need to do additional computation? If it takes more work to get the certificate in those cases, then someone solving a large problem who only cares that the model is unbounded/infeasible (not looking for a proof) might not want to spend the extra CPU cycles.

That said, I agree with you that forcing the user to jump through hoops to get the certificate is not a good idea. It seems particularly problematic to me that in CPLEX the user has to turn off presolving. Assuming the user is not anticipating an infeasible/unbounded result, the user will want presolving on for the gain in computational efficiency; so CPLEX is forcing her/him to solve the problem a second time to get the certificate (which, again, is a PITA if the model is large).

Do you have a blog? An entry about the use of Farkas certificates in DW would be make for interesting reading.

3. If you are using dual simplex and the problem is dual infeasible, then you still get a primal farkas ray for free, so thats no excuse for being sloppy :-)
I have a empty blog, I never found time to write and I think trafic is too limited. Actually I think OR blogs should be gathered in one place such as or-exchange, which would make it a lot more active. Blogs is meant to have a lot of comments and share ideas, but most blogs even Dr. Trick's have few or none comments, which is a shame.
The information you have written about infeasibility/unboundedness is very useful for many people, the right topics for blogs is everything :-) I love blogs with a practical view.

4. I should add, if you use some other phase one strategy than minimizing sum of infeasibilities (such as penalty heuristics), then you don't have a certificate, but on the other hand you can't either for sure declare the problem infeasible. Most commercial codes use some variant of minimizing sum of infeasibilities.

5. I have answered in blog

http://erlingdandersen.blogspot.com/2010/08/exploiting-farkas-certificate-in.html

how to exploit the Farkas certificate in DW decomposition.

6. @Erling: Excellent! I think the post is very informative, and I hope that people working with D-W also find it. It does seem to be one of those useful things that somehow slips through the cracks.

7. Hi Paul.Thanks for the post.Please do you know if I can get the extreme ray/certificate of infeasibility information from CPLEX API for MATLAB?If yes, what code do I use?I have been looking for that information for sometime now. Thanks.
Emmanuel

1. Emmanuel: According to this FAQ, it is not possible from MATLAB.

8. Dear Prof. Rubin,

A friend of mine used dualFarkas method without any concern on the ordering of returned ray values. (he assumes the returned constraints are in order) and he finds the correct result. Is this a coincidence? Should I be careful about mapping? Thanks in advance.

1. I'm not sure. The corresponding call in the C API does *not* pass a vector of constraints, which implies that the returned vector of dual variables corresponds to the constraints in the order they appear in the model. So that would suggest your friend's assumption was sound. On the other hand, the Java API *does* pass a vector of constraints, which suggests that your friend just got lucky (else why bother with that argument?). The C++ API also has that argument.

The higher level APIs (C++, Java, ...) don't always behave the same way that the C API does, so my guess is that CPLEX does not guarantee a particular order for the certificate obtained via dualFarkas(). It's just a guess, but I would be careful about the ordering if I were doing it.

9. Dear Prof. Rubin,

thank you for your post but I have one question, what if I'm getting a feasible solution returned by the solver for the master problem and I need to get the dual prices of the constraints in order to pass them to the pricing problem !? the dualFarkas() method requires an infeasible solution .. what could be the way to get the dual prices in this case ?

Thank you very much in advance,

1. Invoke the getDual or getDuals method on the master problem. You can find a call to getDuals in my sample code above.

If this is your first time commenting on the blog, please read the Ground Rules for Comments. In particular, if you want to ask an operations research-related question not relevant to this post, consider asking it on OR-Exchange.