## Wednesday, January 16, 2013

### Finding All Solutions (or Not)

People occasionally ask how to find all optimal solutions (or, less frequently, all feasible solutions) to a mathematical programming problem. To keep the length of this post reasonable (?), I'll restrict myself to linear and (mixed) integer linear programs, and make only some general observations. Specific methods will be left to future posts.

Let's consider a linear optimization problem of the form $\begin{array}{lrcl} \textrm{minimize} & & f(x,y)\\ \textrm{s.t.} & (x,y) & \in & \mathcal{P}\\ & x & \in & \mathbb{Z}^{m}\\ & y & \in & \mathbb{R}^{n} \end{array}$ where (without loss of generality) we are minimizing a linear function $f$ over a convex polytope $\mathcal{P}$. For a linear program (LP), omit the integer variables $x$; for a pure integer linear program (ILP), omit the continuous variables $y$. Again without loss of generality, we assume that the feasible region $\tilde{\mathcal{P}}=\mathcal{P}\cap\left(\mathbb{Z}^{m}\times\mathbb{R}^{n}\right)$ is not empty (else there is nothing to find) and that $\mathcal{P}$ (and therefore $\tilde{\mathcal{P}}$) is bounded. In any realistic situation, $\mathcal{P}$ can be bounded by adding "box constraints" (bounds) for $x$ and $y$. Unrealistic situations are left to the reader as an exercise. Since the continuous hull of $\mathcal{P}$ (relaxing the integrality constraints) is closed (courtesy of convexity) and bounded (by assumption), and $f$ is continuous (courtesy of linearity), and since the feasible region $\tilde{\mathcal{P}}$ is nonempty (by assumption), at least one optimal solution must exist.

The two questions are equivalent. Suppose that we are looking for all feasible solutions (i.e., the contents of $\tilde{\mathcal{P}}$). That's the same as finding all optimal solutions with a constant-valued objective function $f(x,y)=c$. On the other hand, finding all optimal solutions to the problem above is equivalent to finding all points in the modified feasible region $\left\{ x\in\mathbb{Z}^{m},y\in\mathbb{R}^{n}:(x,y)\in\mathcal{P}^{*}\right\}$ where $\mathcal{P}^{*}=\mathcal{P}\cap\left\{ (x,y):f(x,y)\le f^{*}\right\}$ and $f^{*}$ is the optimal objective value in the original problem. Since we obtained $\mathcal{P}^{*}$, the set of optimal solutions, by intersecting $\mathcal{P}$ with a half-space, $\mathcal{P}^{*}$ is another convex polytope.

Henceforth I'll restrict myself to the question of finding all optimal solutions.

Uncountably infinite is a lot. Start with a linear program (no $x$). If $\hat{y}$ and $\tilde{y}$ are two optimal solutions with objective value $f^{*}$, then any convex combination $\lambda\hat{y}+(1-\lambda)\tilde{y}$ with $0\le\lambda\le1$ belongs to $\mathcal{P}$ (by convexity) and has the same objective value $f^{*}$ (by linearity of $f$). So right away we have an uncountable number of optima. With a pure integer program (no $y$), the situation is rosier: since we have assumed that $\mathcal{P}$ is bounded, the feasible region $\tilde{\mathcal{P}}$ will contain a finite number of integer lattice points, and so the number of optimal solutions will be finite. (Note, however, that "finite" does not preclude "annoyingly large".)

With a mixed integer linear program (MILP), things are trifle more complex. The projection of the set of optimal solutions onto $\mathbb{Z}^{m}$, i.e., the set of $x$ for which a $y$ exists such that $(x,y)$ is optimal, will be finite. If for each such $x$ there is exactly one $y$ yielding an optimal solution, the number of optima is finite. As soon as we find one $x$ such that $(x,y)$ and $(x,y')$ are both optimal with $y\neq y'$, we're back to having an uncountable number of optima.

Putting all that together, for ILPs we can ask to find all optima. For LPs, the best we can do is ask to find all extreme points (vertices) of $\mathcal{P}$ that are optimal. For MILPs, we can ask for all $x^{*}$ that correspond to optimal solutions, and for each such $x^{*}$ we can ask for a single $y^{*}$, or at most all extreme points of $\left\{ y\in\mathbb{R}^{n}:\left(x^{*},y^{*}\right)\in\mathcal{P}^{*}\right\}$ (yet another polytope).

This could take a while. Let's consider the LP case again (no $x$). Per the Wikipedia entry, the Avis-Fukuda algorithm for enumerating the vertices of a polytope $\mathcal{P}$ expressed as the solution to a set of linear inequalities in $y$ has time complexity of $O(knv)$, where $k$ is the number of nonredundant inequalities defining $\mathcal{P}$, $n$ is the dimension of $y$ and $v$ is the number of vertices to be enumerated. (Sorry, my notation is slightly different from the Wikipedia's. Deal with it.) This does not look too bad ... until you stop to think about $v$. If the polynomial is a hypercube, $k=2n$ and $v=2^{n}$, so we're looking at $O\left(n^{2}2^{n}\right)$ time. That's fine for small $n$, but it blows up fairly quickly as $n$ increases.

What about the integer case? Well, ILPs and MILPs are NP-hard in general,
and that's just to find one optimal solution. Finding all optimal
solutions would be, what, NP-NP-hard?

Is this really what you want? This brings me around to motivation. Why do you want all optimal (let alone all feasible) solutions? One possible answer is to offer a decision maker some choices, charitably assuming that more than one optimum exists. (Here is a related post.) Another is that there are criteria not included in the model that distinguish the "quality" of various solutions (so not all optimal solutions are equally good, even though they are all equally optimal).

My suggestion in both those cases is to find one optimal solution $\left(x^{*},y^{*}\right)$ with objective value $f^{*}$, then find "epsilon-optimal" solutions satisfying some diversity criterion. We pick our favorite $\epsilon>0$, add the constraint $f(x,y)\le f^{*}+\epsilon$, and change the objective to something else. Possibilities include (but are by no means limited to) maximizing the $L_{1}$ distance of the next solution from $\left(x^{*},y^{*}\right)$ or from the centroid of previously found solutions, maximizing the minimum distance from all previously found solutions, maximizing or minimizing one particular variable (e.g., $\max x_{i}$ or $\min y_{i}$), or optimizing a randomly generated linear function of $x$ and $y$.

If you have some alternate criterion, not part of the original model, that you were going to use to pick among competing optima, and if that criterion is linear (or can be approximated by a linear function), it makes an obvious choice for the second stage objective.

1. "Finding all optimal solutions would be, what, NP-NP-hard?"

#P-hard, I think. http://en.wikipedia.org/wiki/Sharp-P http://en.wikipedia.org/wiki/Sharp-P-complete

1. Looks right, but I glaze over quickly when it comes to complexity terminology.

2. Hello Prof Rubin,
I'm not able to understand the first two sentences of the section "Uncountably infinite is a lot.". I certainly did not understand how the convex combination with the lambda part has the same objective value. Can you please help me out in understanding that part?

1. Suppose your LP is $$\min c'y$$ s.t. $$Ay \le b$$ with $y\ge 0$. Let $\hat{y}$ and $\tilde{y}$ be two optimal solutions with objective value $f*$, and let $\lambda \in [0,1]$. Then $$c' [\lambda\hat{y}+(1-\lambda)\tilde{y}]=\lambda(c'\hat{y})+(1-\lambda)(c'\tilde{y})=\lambda f*+(1-\lambda)f*=f*$$ and $$A[\lambda\hat{y}+(1-\lambda)\tilde{y}]=\lambda(A\hat{y})+(1-\lambda)(A\tilde{y})\le\lambda b+(1-\lambda)b=b.$$

3. Ahh...thank you very much Prof Rubin. This cleared my doubts and you were very fast in replying me back. Thanks for that as well.
One off-topic question I have is that, I found quite a lot of posts very helpful for my thesis work and I want to cite your work but I can't find a proper Bibtex format for my references. Do you think you assist me regarding this issue?

1. Maybe the method suggested by the Wikipedia: https://en.wikipedia.org/w/index.php?title=Special:CiteThisPage&page=LaTeX&id=413720397 (subsection labeled "BibTeX entry")?

4. using a commercial solver, is it possible to get all optimal vertices?

1. If you're asking whether commercial solvers contain routines/commands to find all solutions, I've not aware of any that do.

5. thanks for the post. One reason why you may want all the optimal solutions is that you know a set of optimal solutions exist and you are interested in the properties of the set. For example, when working with partially identified parameters in statistical estimators.

1. Fair point. What I might describe as theoretical or algorithmic investigations could legitimately require identifying the full set of optima. My comments should have been targeted at people solving problem instances for the purposes of making decisions (which may not benefit from a panoply of optimal choices, and will almost sure benefit from reaching a conclusion in a timely manner).

6. Hello, I'm using CPLEX Optimization Studio 12.6 to find the optimum solution of a CP. I would like to find a way to get all the optimal solutions for that problem, because the default settings seem to return just one optimal solution. Can you help me or give me an example ? thanks

1. I'm not sure, since I don't use CP Optimizer. For constraint satisfaction problems, this is covered in the user manual for CP Optimizer. Look in "Search in CP Optimizer > Searching for solutions > Solving a satisfiability problem". You can use solve() to find the optimal objective value, add a constraint limiting the objective to that value, then do a new search as a constraint satisfaction problem and call next() repeatedly. That may not be the most efficient way, though. You probably should search the CP Optimizer support forum and, if you don't find it (I suspect it's a FAQ), ask there.

2. It works,
really thank you for this idea!
Nice to meet you Helpful Prof Rubin.

7. Having a similar problem with lp_solve. I know my models contain a few (well not too many) optimal solutions but the java wrapper I am using does not expose any methods for finding more than one.

If I find the first one, add the extra constraint using the optimal value, can I then start toggling the variables one by one (they are binary variables) while I resolve? That'll give me 2^n possible iterations? which just seems like a hack :-(

1. I'm not familiar with lp_solve (I know what it is, but I don't recall ever using it), so I won't be much help with specifics. If your model has binary variables, you can try adding a no-good constraint to eliminate the first optimal solution, then solving again. If the new optimal solution has a worse objective value, you're done; if the objective is the same, record it, add another no-good constraint, and iterate. It's also a hack, but less brute force than toggling binaries.

2. Thanks Paul, Yes that sounds a lot better. Will try that.

3. Yep. That works great for the type and size of problem instances I am dealing with. Many thanks Paul :-)

8. Thanks for the post!

One other potential approach:
If you restrict yourself to finding all optimal solutions to an LP, you could run an interior point method (without crossover), which would find the analytic center of the optimal face. If the IPM finds an extreme point then you have the unique optimal solution, otherwise the set of binding constraints at the analytic center defines the optimal face and you can extend this set of constraints by combining it with any of the non-binding constraints to form a basis, which would be optimal wherever it is feasible.

Due to intermittent spamming, comments are being moderated. 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 Operations Research Stack Exchange.