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.addMinimize(
  primal.scalProd(new double[] {1.0, 1.0, -2.0}, vars));
IloConstraint c1 =
  primal.addGe(
    primal.scalProd(new double[] {1.0, -2.0, -1.0}, vars),
    3.0);
IloConstraint c2 =
  primal.addGe(
    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.