Operations Research Stack Exchange is a question and answer site for operations research and analytics professionals, educators, and students. It only takes a minute to sign up.
Sign up to join this community
It is known that, in a linear program with $k$ constraints, there exists a basic feasible solution in which at most $k$ variables are non-zero. How can I find such a solution?
Is there a polynomial-time algorithm for this? I know that the simplex method finds these solutions, but it is not polynomial time.
In practice, how can I tell the common Python libraries for linear programming (particularly: cvxpy and scipy.linprog) to find me a solution with a small number of zeros?
As also mentioned by Philipp Christophel: polynomial-time does not mean that the algorithm is better in practice. I strongly suggest that you try the simplex method first, as it does exactly what you want, implementations are readily available, and it might very well be the fastest method for your purpose.
That being said: yes, you can obtain a basic feasible solution in polynomial time with the ellipsoid method or with interior point methods.
Khachiyan (1980) describes how the ellipsoid method can be used to find an optimal basic feasible solution in polynomial time. This result is interesting theoretically, but you definitely do not want to implement this, as practical performance is terrible.
Interior point methods by themselves return a feasible solution that is close to the optimum, but it need not be a basic feasible solution. To go from an interior point solution to a basic feasible solution, crossover algorithms can be used. These are implemented in most commercial solvers, e.g., see the crossover page in the CPLEX documentation and in the Gurobi documentation. So if simplex is not fast enough, you may try interior point/barrier with one of these solvers and enable the crossover.
Given optimal primal-dual solutions, crossover is possible in polynomial time. Megiddo (1991) gives a strongly polynomial algorithm to find an optimal basic feasible solution. However, I do not know how numerical imprecision in the numerical solution affects this result. And even then, I highly doubt that solvers implement a polynomial time algorithm, as algorithms without complexity guarantees may perform significantly better. This was observed by Bixby and Saltzman (1994), for example.
The more precise statement is that in a basic solution at most k out of the n variables are away from their bounds. If they have a non-zero lower or upper bounds, non-basic variables can also be non-zero. But that might not be the case in what you are working on.
Regarding (1): If you google "polynomial algorithm for linear programming" you can probably have a lot of fun reading. The short version: Interior point (barrier methods) and the ellipsoid method are polynomial, but only interior point methods are good in practice. The simplex method is not known to be polynomial but especially the dual simplex works very well in practice. I wrote more about this here already: State-of-the-art algorithms for solving linear programs
Regarding (2): I can't say anything specific about these packages, but interior point methods are also different from simplex methods in that they don't necessarily return basic solutions. If the problem is degenerate (and most problems are) they return a solution that sits on a face of the polyhedron instead of a vertex (basic solution). Thus interior point methods will return solutions with "a small number of zeros".
In practice, I have had great success in Python with LpSolve 5.5, I have never had to resort to interior point methods, yet. Dual simplex has been sufficient for me- even for integer problems.
But, for your second question, the global best solution depends on your cost function. If the best solution is reached with a large number of zeros - that is still the solution. In rare / integer cases you may find multiple solutions for the same objective value.
Instead, you can add an integer choice variable like 'Is $x_i$ allowed to be non-zero?' Depending on the number of decision variables you have, you can add an additional variable to each of them:
For each decision variable $x_i > 0$, add the constraint:
$$x_i \le Q_i \cdot \text{Inf},$$ where $Q_i$ is $\{0, 1\}$ and $Q_i=1$ allows $x_i$ to be greater than $0$.
Now sum $\sum\limits_iQ_i$ (for all $i$) to your objective function, with a weight that doesn't significantly affect your objective. And this should naturally lead to an objective function where fewer zeros are better.
If your problem has many decision variables this approach may just make the problem too large to solve in a reasonable time, especially since you're adding integer decision variables.
