Linear programming feasibility problem with strict positivity constraints

There is a system of linear constraints ${\bf Ax} \leq {\bf b}$ . I wish to find a strictly positive vector ${\bf x} > 0$ that satisfies these constraints. That means, $x_i > 0$ is required for every component $x_i$ of ${\bf x}$. How can I use an LP solver to find such a strictly positive vector ${\bf x}$ (or confirm that no ${\bf x}$ exists)? I cannot simply introduce another system of constraints $x_i > 0$, because equality must always be permitted in an LP—but I can use the LP solver several times, with changing objective functions. I think I should use the slack variable method, but I don't know how.

migrated from stackoverflow.comJan 24 '12 at 7:15

This question came from our site for professional and enthusiast programmers.

You can circumvent the problem of choosing a small $$\epsilon>0$$ by being a bit more ambitious: Try to find $$\mathbf{x}$$ such that $$\mathbf{Ax}\leq \mathbf{b}$$ and that the smallest entry in $$\mathbf{x}$$ is largest possible.

To that end, introduce a new variable $$\mathbf{y} = \begin{bmatrix} \mathbf{x}\\ \epsilon\end{bmatrix}\in\mathbb{R}^{n+1}$$ (if $$\mathbf{x}$$ was in $$\mathbb{R}^n$$) and solve the following problem by an LP-solver $$\max_y [0 \dots 0\ 1]\cdot \mathbf{y}\quad \text{s.t.}\quad [A\ \mathbf{0}]\mathbf{y}\leq\mathbf{b}\quad \text{and}\quad \mathbf{0}\leq \begin{bmatrix}1 & 0 & \cdots & 0 & -1 \\ 0 & 1 & \cdots & 0 & -1 \\ & \ddots & &\vdots \\ 0 & \cdots & & 1 & -1 \end{bmatrix}\mathbf{y}.$$

This is a reformulation of the following problem: $$\max \epsilon \quad\text{s.t}\quad \mathbf{Ax\leq b}\quad\text{and}\quad \mathbf{x}\geq\epsilon\mathbf{1}.$$

• well done, this is equivalent to a trick a coauthor and I just used in a recent paper, and definitely superior to the approach I suggested. – Aron Ahmadia Jan 24 '12 at 12:26
• Agreed. Well played, sir. – Geoff Oxberry Jan 24 '12 at 13:00
• The reformulated problem may have an unbounded objective in cases where the answer to the original problem is trivial. For example, if the system of constraints is just $x \ge -1$. That is fine as long as you check for feasible, optimal or unbounded in the return status of your lp solver, or explicitly bound the $\epsilon$. – David Nehme Feb 17 '12 at 3:50
• @DavidNehme: One can add the constraint $y_{n+1}\le 1$ to get a bounded objective. – Arnold Neumaier Jun 22 '12 at 19:51

For an LP feasibility problem, I wouldn't use standard simplex. Standard primal (or dual) simplex algorithms will only visit the vertices of the feasible set of the primal (or dual) problems.

Let the feasible set of the problem you actually want to solve be $F = \{\mathbf{x}: \mathbf{Ax} \leq \mathbf{b}, \mathbf{x} > \mathbf{0}\}$, and suppose instead you were to solve the problem ($F_{\varepsilon}$):

\begin{alignat}{1} & \min_{\mathbf{x}} \quad 0 \\ \textrm{s.t.} & \quad \mathbf{Ax} \leq \mathbf{b} \\ & \quad \mathbf{x} \geq \varepsilon \cdot \mathbf{1}. \end{alignat}

The closest approximant of the problem you want to solve is $F_{0}$, which admits slightly too many points. The problem is that the boundary of the positive orthant (i.e., the set $B = \{\mathbf{x}: \mathbf{x} \geq 0, \exists{i}: x_{i} = 0\}$ could make up part of the boundary of the feasible set of $F_{0}$. We'd like to exclude those points. One way of doing that is to do what Aron suggested, which is to set $\varepsilon$ to some small positive value, and then use any standard LP algorithm. This strategy is a good one, and will probably work in a wide variety of situations. However, it will fail if $F_{\varepsilon}$ is infeasible. We know that $F_{0} \subset F \subset F_{\varepsilon}$ for all $\varepsilon > 0$ (to abuse notation and refer to a feasible set by its corresponding problem), and it's possible that even if you pick small positive values of $\varepsilon$, the LP solver will indicate that your LP is infeasible.

For an LP solver, I'd use any interior point algorithm for LPs that starts with a feasible point and stays feasible, which is another way to exclude points in $B$. You needn't have to supply a feasible point to these algorithms; standard solvers will do it for you. Methods like affine scaling, potential reduction, and barrier methods set up auxiliary LPs that will find feasible solutions, and the iterates for these algorithms traverse the interior of the feasible region. You only need to locate one point in your feasible region, so as long as the auxiliary problems used by the LP solvers locate a feasible point for your problem, and that feasible point is strictly positive, you should be all right. If solving $F_{\varepsilon}$ fails for small positive values of $\varepsilon$, you might still be able to use these methods to locate a strictly positive feasible point within $F_{0}$.

Don't use simplex, though, because it will only explore the vertices of $F_{\varepsilon}$, which is exactly what you want to avoid doing.

Feasibility problems are a slightly trickier game than general linear problems, which you have noted. If you are solving approximately (by using a floating-point representation of the system of equations and constraints), it is legitimate to require $x_i >= \epsilon$, where $\epsilon$ is some very small numerical value, big enough to assure that $x_i$ actually lives in $\Re_+$, but small enough that a solution on the boundary is not considered.

You might have to adjust $\epsilon$, and your solution will be qualified to "within a factor of $\epsilon$", but this is sufficient for many situations.

The answer given by aeismail is to be read carefully, regard the lp

$\max( x_1 + x_2)$

s.t.

$x_1 + x_2 \le 1$

$x_1,x_2 \ge 0$

It has solutions $(1,0)$ and $(0,1)$ as well as others (degenerated). The generality of the question implys that these cases need to be treated as well.

Since you are able to choose you objective function, you could try to modify it iteratively. E.g. Start with all coefficients for all variables equal to one, check wether you gain an approprate solution. If one variable is zero, rise it's coefficient and start again...

Though I can not give a mathematical prove that this works (or a well defined procedure how to modify the objective function). I hope this helps :)

• However, if you have a large number of degenerate solutions, how would you deal with this numerically? Won't pretty much any numerical solver throw up a warning (or worse) about solving this problem? – aeismail Jan 24 '12 at 9:42
• No, they won't; they'll just return the first optimal solution encountered. The way that you would continue to generate solutions is to add cutting planes (or other constraints) that exclude previously calculated optimal solutions. In this case, adding such cutting planes would enable you to return a discrete approximation of the infinite set of optimal solutions. – Geoff Oxberry Jan 24 '12 at 10:40
• I would view that as a strange programming decision; why wouldn't you want to tell the user that the objective function was doing something weird in the neighborhood of the reported solution? For a nonlinear solver, I could see there being an issue with figuring out what's going on; but shouldn't that be easier to tell with a linear system? – aeismail Jan 24 '12 at 13:19
• I'd have to think about how one would detect degeneracy by actually constructing problems, but typically, users want an optimal solution, so the most important information for an LP is to return if the solution is optimal, feasible (but not optimal), infeasible, or unbounded. (These statuses are, in fact, what a solver like CPLEX would return.) Degeneracy is primarily a theoretical issue; the only reason it would be discussed in a numerical context is either in algorithm design or in practice, to note that degeneracy typically slows down a solver. – Geoff Oxberry Jan 24 '12 at 20:10