Saturday, July 23, 2011

Farkas Certificates in CPLEX


This is a post about linear programs, specifically infeasible linear programs, and how to deal with them in CPLEX. In response to a question on a CPLEX forum, I played around a bit with Farkas certificates as implemented in CPLEX 12.2. I've never really worked with them, as I have a tendency to produce feasible models. The posted question, though, aroused my curiosity. It asked how to extend a Farkas certificate to a ray for the (unbounded) dual when the primal problem contains bounds on variables.

A Farkas certificate is a vector of constraint multipliers that proves a primal linear program to be infeasible. The name comes from Farkas's Lemma. Assume, without loss of generality, that we have a linear program of the form\[
\textrm{(P) }\begin{array}{lrcl}
\textrm{minimize} & 0\\
\textrm{subject to} & Ax & \ge & b\\
& x & \ge & l\\
& x & \le & u
\end{array}.
\]We are interested here only in the feasibility of (P), so we may as well assume the objective function (irrelevant to feasibility) is zero. Vectors $l$ and $u$ are lower and upper bounds for the variables $x$. To simplify the discussion, we will assume that $l$ and $u$ are finite (but not necessarily nonnegative). What will transpire extends easily to infinite bounds.

Assume that the lower and upper bounds are entered as bounds, rather than as constraints (so that the solver sees $Ax\ge b$ as the only set of constraints). A Farkas certificate for (P) will be a vector $y\ge0$ with the following property. Set $q'=y'A$ and define a vector $z$ of the same dimension as $x$ as follows:\[
z_{j}=\begin{cases}
l_{j} & \textrm{if }q_{j}<0\\
u_{j} & \textrm{if }q_{j}\ge0
\end{cases}.
\](Actually,the definition of $z_{j}$ is arbitrary and irrelevant if $q_{j}=0$; I chose $u_{j}$ just for concreteness.) The key property of the Farkas certificate is that it will satisfy\[
q'z=y'Az < y'b.
\]Now suppose that $x$ is any feasible solution, so that in particular $l\le x\le u$.  Then \[
q_{j}\ge0\implies z_{j}=u_{j}\ge x_{j}\implies q_{j}z_{j}\ge q_{j}x_{j}
\]and
\[
q_{j}<0\implies z_{j}=l_{j}\le x_{j}\implies q_{j}z_{j}\ge q_{j}x_{j}.
\]So $q'x\le q'z<y'b$. If $x$ is feasible in (P), though, we have $Ax\ge b$; given that $y\ge0$, it must follow that $q'x=y'Ax\ge y'b$, a contradiction. Thus the existence of a Farkas certificate for (P) tells us that there cannot be any feasible solutions $x$ to (P).

Now let's tie this to the dual (D) of (P), which is\[
\textrm{(D) }\begin{array}{lrcl}
\textrm{maximize} & y'b+v'l-w'u\\
\textrm{subject to} & y'A+v'-w' & = & 0\\
& y,v,w & \ge & 0
\end{array}
\]if we now include the bounds on $x$ as constraints (and treat $x$ as unrestricted in sign). The feasible region of (D) is a cone with vertex at the origin. Denoting by $h^{+}$ and $h^{-}$ the positive and negative parts respectively of a vector $h$ (i.e., $h_{j}^{+}=\max(h_{j,}0)$ and $h_{j}^{-}=-\min(h_{j},0)$), let $y$ be a Farkas certificate for (P), let $q'=y'A$ as before, and let $v=q^{-}$ and $w=q^{+}$. Then $(y',v',w')\ge0$ and $(y'A+v'-w')'=q+v-w=q+q^{-}-q^{+}=0$;
so $(y',v',w')$ is a feasible solution to (D). Moreover, based on our definition of $z$ above, $q'z=-(q^{-})'l+(q^{+})'u=-v'l+w'u$, and so $y'b+v'l-w'u=y'b-q'z>0$. By extending the Farkas certificate $y$ with $v=q^{-}$ and $w=q^{+}$, we obtain a feasible solution to the dual with a positive dual objective value. Since the dual is a cone, $(y',v',w')$ is a recession direction along which the dual objective is unbounded.

Let's go in the other direction. Suppose that (D) is unbounded (making (P) infeasible), and let $(y',v',w')$ be a feasible dual solution with objective value $y'b+v'l-w'u>0$. Let $q'=y'A=w'-v'$ and define $z$ as before. Then $l\le z\le u$, so\[
q'z=w'z-v'z\le w'u-v'l<y'b
\](the last step since the dual objective is positive at $(y',v',w')$); thus $y$ must be a Farkas certificate for (P).

To answer the original question (finally!), suppose we are solving an infeasible primal problem in CPLEX. CPLEX provides a mechanism (the IloCplex.dualFarkas method in the Java API) to retrieve a Farkas certificate $y$. To extend that to the corresponding dual multipliers $v$, $w$ for the lower and upper bounds respectively, compute $y'A$ and assign the positive parts to $w$ and the negative parts to $v$. It would be nice if there were a more direct method (one that would avoid the multiplication $y'A$), but I do not know of one.

One last note: the point of extending the Farkas certificate to a full dual solution is usually to obtain an extreme ray of the dual. As shown above, any solution to (D) with positive objective value provides a Farkas certificate, not just extreme rays. So we rely on the mechanism by which the solver generates the certificate to pick one corresponding to an extreme ray.

9 comments:

  1. Hi!

    I think you should use (A | I | I) as constraint matrix of (D) instead of (A | I | -I) for the following reason: If the bounds of (P) were constraints, the constraint matrix of (P) would be (A' | I' | I')' (or the upper bound-constraints should be -x >= -u). w should then be a non-positive variable ( w<=0 ).

    The answer to the original Problem is then: Set v = w = 0. Compute q' = -y'A. For i = 1..n: If q[i] >= 0 then v[i] = q[i] else w[i] = q[i].

    Best regards,
    Thomas

    ReplyDelete
  2. Thomas,

    Making the third block of the constraint matrix of (D) I, making w <= 0, and adding w'u in the dual objective is equivalent to what I did (third block -I, w >= 0, subtract w'u in the objective). I guess it's a matter of taste which way you write it; I lean toward all nonnegative variables where possible. Among other things, I think it makes the connection between q, v and w a bit cleaner.

    I think, though, that your q may have the opposite signs of what CPLEX (and other solvers?) produces.

    ReplyDelete
  3. Hi!

    Thank you for your reply.

    Of course, this is equivalent. We are working on benders cuts and I spent a lot of time in finding a mistake in our code until I saw that I got confused with the signs. Maybe others have similar problems.

    For me, q was just some variable. Maybe I should write: Set v = w = 0. Compute q' = y'A. For i = 1..n: If q[i] <= 0 then v[i] = -q[i] else w[i] = -q[i].

    By the way: Thank you for this page. It helped us a lot!

    Best regards,
    Thomas

    ReplyDelete
  4. You are most welcome - and thank you for the comments!

    Cheers,
    Paul

    ReplyDelete
  5. This comment has been removed by the author.

    ReplyDelete
  6. Hi Paul.Thanks for the post. I have two questions:
    1) In CPLEX, if the variables bounds are explicitly entered as constraints, rather than bounds, then can the dualFarkas return the complete extreme ray? As for the method you stated in this post, I wonder how to calculate y′A in CPLEX?
    2) I found the following descriptions about dualFarkas(), at the CPLEX Optimizer (https://www.ibm.com/developerworks/community/forums/html/topic?id=77777777-0000-0000-0000-000014433103). Does this mean we can obtain the extreme ray by combining the results returned by dualFarkas() and getReducedcost()?

    CPXdualfarkas() is somewhat tricky to understand. The method provides a set of dual multipliers y', such that the aggregated row y'Ax >= y'b is a valid inequality (the signs of y' match the inequality sense of the constraints) that proves infeasibility of the system. Namely, this aggregated row a'x >= b' cannot be made feasible even if you substitute the "most feasible" bound of x_j for each variable x_j. In other words, max{a'x | l <= x <= u} < b'. In this sense, y' is not a dual ray (because y'A does not need to be 0), but you can extend it with suitable reduced cost values (i.e., duals for the bounds) to a dual ray.

    For example, consider the following system:

    1x1 + 3x2 - 2x3 >= 5
    1x1 - 4x2 + 4x3 >= 6

    with bounds 0 <= x1 <= 3, 0 <= x2 <= 5, 0 <= x3 <= 1.

    Then, the dual vector y'=(1,2) yields the aggregated row y'Ax >= y'b that reads

    3x1 - 5x2 + 6x3 >= 17.

    This shows the infeasibility of the model, because substituting the "most feasible" bounds for x yields

    3*3 - 5*0 + 6*1 = 15 < 17.

    Now, if you extend the dual vector y' with reduced costs r' = (-3,5,-6), you get a dual ray:

    1x1 + 3x2 - 2x3 >= 5 (y'_1 = 1)
    1x1 - 4x2 + 4x3 >= 6 (y'_2 = 2)
    x1 <= 3 (r'_1 = -3)
    x2 >= 0 (r'_2 = 5)
    x3 <= 1 (r'_3 = -6)

    Summing this up with the weights of the dual ray yields

    0x1 + 0x2 + 0x3 >= 2,

    which proves the infeasibility of the model.

    ReplyDelete
    Replies
    1. For the first part of (1), I think that putting the bounds (both upper and lower) in as explicit constraints should work, with one caveat. When you declare the variables, give them lower (upper) bounds strictly less (greater) than the actual bounds (the ones you enter as constraints). For instance, if you are going to put $0 \le x \le 5$ in the constraints, give $x$ a lower bound of -1 and an upper bound of 6. That's to ensure that the bounds are nonbinding and that the actual dual multiplier shows up as the dual of the constraint and not as the reduced cost of a binding bound.

      For the second part of (1) (calculating y'A), I would typically do this in my code. If you solve the dual explicitly, rather than solving the primal and then call dualFarkas(), then you can use the getValue() method to get the value of the constraint expressions (which would be y'A).

      For (2), that's how I interpret the post you cited. I have not actually tested the method, though. Also, I think you have to be careful about the signs of the reduced costs, since they relate to which bound (lower or upper) is binding.

      Delete
    2. Thank you for your reply. It helps me a lot!
      I still have some more questions on the calculation of y'A.
      If we solve the dual explicitly, I think we do not have to call dualFarkas( )since getRay( )can give us the extreme ray. If we solve the primal and then call dualFarkas(), I don't think calculating y'A is a easy job. I don't know there is some tools in CPLEX to get the coefficients vector for a specific variable in all constraints, which can be used to obtain a column of y'A. The only way in my mind is to get the expression in each constraints and then go through it with an iterator to construct the coefficients vector for each variable. Do you have any idea?

      As for (2), I tried it in CPLEX, with the example I cited, but combining the results returned by dualFarkas() and getReducedcosts() doesn't work. In fact, the constraint multiplier (y') CPLEX returned is (1, 0.75), I think the reduced costs is relevant to the objective, as the example did not give one, I tried two objectives: 0 and x1 + x2 + x3, neither of them works. I'm trying to figure this out, do you have any idea?


      Best regards

      Delete
    3. Regarding the computation of $y^\prime A$, I think that if you build the model using an instance of IloLPMatrix (rather than row-wise using IloRange instances), you can then use IloLPMatrix.getIndex() to map variables to column indices and IloLPMatrix.getCols() to get $A$ column-wise. I have to confess, though, that I never use the matrix modeling approach, so I could be wrong here.

      For your second question, I'm afraid I don't have any ideas off the top of my head, and right now don't have time to give it any further study.

      Delete

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.