# Non-sparse solution for a linear programming problem

I have formulated a linear program with equality, inequality and non-negativity constraints. My objective function is the minimization of a linear combination of decision variables (with different coefficients for different variables). I have 7803 decision variables, 51 inequality constraints and 111 equality constraints. I used the simplex method and found the solution. In the solution the values for most decision variables is 0 (like 98%), which is not surprising, based on the model formulation. Because of the nature of the problem, I need to have a way larger ratio of the decision variables with non-zero values in the solution. Is there any formulation technique that can do this for me? In other words, I need to modify my formulation (I guess probably in the OF) to find a less-sparse solution, at the cost of getting a little far from the optimum point. I am okay if it ends up in a non-linear optimization problem.

It is not surprising that your current solution has so many variables with value zero especially if your solver uses the simplex method instead of an interior point method. When you say "Because of the nature of the problem, ..." it means that your current model does not represent your actual problem (it is not the right model) and you should formulate the problem differently. Anyway, here is what you should consider. Given that your current problem is a linear program we can consider two cases:

1. Your linear program has a unique optimal solution.

2. Your linear program has alternative optimal solutions.

In the first case, it is impossible to increase the Number of Positive Variables (NPV) without making the optimal value worse. So you need to change your problem if you really want to increase the NPV. You need to solve a two-objective optimization problem to maximize NPV while optimizing your original objective function. There are different ways to formulate and solve multi-objective problems. For example, you may consider a convex combination (or weighted average) of the two objective function as your new objective function and solve the problem by typical methods used for single-objective problems.

In the second case, you may find another optimal solution for the problem with a better value of the NPV. To find such a solution, you may define a secondary objective for the problem. The new problem is considered as a lexicographic multi-objective problem which can be solved taking two approaches. The first approach is to define an objective function that maximizes the NPV without influencing the primary objective and then use typical single-objective methods to solve the problem (finding such an objective function is sometimes difficult but may be possible). The second approach is to solve a sequence of single-objective problems (once you solve the original problem, you transmit some information to a second problem whose objective is maximizing the NPV).

As you see, in both cases you need two objective functions. Given the number of variables and constraints of your problem it should be easy to solve even if you add some integer variables to your model. Note that your problem becomes a mixed integer program that may be solved by MILP solvers. I suggest to define binary variable $y_i\in\{0,1\}$ corresponding to each $x_i$ and add the constraints $l_i\, y_i \le x_i \le u_i\,y_i$ for each $i$, where $u_i$ is an upper bound on the variable $x_i$ (or a big-M parameter) and $l_i$ is any positive number as the lower bound of $x_i$ (e.g., $l_i = \epsilon$). Note that such a lower bound is effective only if $y_i=1$ and this is a decision that is made by the optimization problem itself. You can now consider $\max (\sum_{i} y_i)$ or $\min (-\sum_{i} y_i)$ as another objective of the problem. You can formulate the problem as a MILP (not a nonlinear program).

Note that you should still find the right single-objective model of your problem and what I mentioned here may give you some insight to do so.

• Thank you for your comprehensive answer. I think doing a two-objective optimization problem to maximize NPV while optimizing my original objective function is the exact thing I need. I found this Matlab function (GAMULTIOBJ) that solves these problems using GA but it seems it doesn't work with integer constraints (that I will have for yi's). Do you know how I can solve multi-objective mixed integer programs? Or are you familiar with a simple method for calculating weights for the two objective functions, so I can use the weights to make it a single-objective mixed integer program and solve it?
– Fred
Sep 21 '17 at 21:33
• I suggest using FPBH for solving multiobjective mixed integer linear programs. It is the most recent and trusted package that I know about. You can find the package and the instructions from github.
– user477602
Sep 21 '17 at 23:53
• I used Johan's trick and it is working so far. I will be able to verify my results later. If I find out the current solution is not good, I will come back and try your suggested method. Thank you for the time you spent on it.
– Fred
Sep 22 '17 at 20:00

A simple trick is to solve the LP resulting in the solution $x^{\star}$. As you don't like this solution, you now solve a new LP with the old constraints and additional constraint that $c^Tx = c^Tx^{\star}$ (i.e., just as good) but minimize the distance to what you think is a good solution (such as $\left\| x-1 \right\|_2^2$).

• isn't the new constraint exactly what old equality constraints were? This is what I understand from cTx=cTx⋆. Also for the new OF, can you make an example please? Like, if the original OF is A1x1 + A2x2 + ... AnXn, what will be the new OF?
– Fred
Sep 21 '17 at 20:46
• You keep all old constraints, but add an equality on the objective, to get a solution with same objective but hopefully more non-zeros (old objective is $c^Tx$, new objective is quadratic here, but you could use linear norm such as sum of absolute values) Sep 21 '17 at 20:51
• Okay, now I understand that I'll need to add the old OF as a equality constraint to the new problem. The new OF will be then (x1-1)^2+(x2-1)^2+...(x7803-1)^2? Is this what you mean? Can you explain how this OF will make the number of non-zero variables in the solution fewer? Thanks.
– Fred
Sep 21 '17 at 21:22
• Because the only thing we are trying to do is to force all variables towards 1, which is not zero. It is a heuristic. There is absolutely no guarantee that it will work Sep 22 '17 at 5:55
• What you probably should do is to tweak solver options. A significant effort is placed in the solver to achieve precisely what you don't want. A final step in advanced solvers is the cross-over, when one tries to move to a vertex of the feasible set, in case the final solution is on a plane or line or something, or in the interior. My advice would be to turn on an interior-point algorithm in the solver, turn off cross-over, and reduce the optimality tolerances to force the solver to terminate prematurely when it is still far into the interior, and thus should have more non-zeros. Sep 22 '17 at 5:59

You could just add constraints of the form $x_i \ge \epsilon$ for some of the variables, where $\epsilon$ is the smallest positive value you would allow as "non-zero". You might choose the variables with smallest shadow prices in your optimal solution of the original problem.

• That enters subjectivity in the problem and I don't know how to justify that. Also, then won't all the decision variables that were originally zero and are forced to be greater than epsilon, equal epsilon in the solution? That would not be helpful. I need an alternate formulation that is more flexible and splits the quantities among more variables, rather than assigning them all to a limited number of variables.
– Fred
Sep 20 '17 at 22:06
• From your description, it's really hard to guess what the actual problem is. Why do you need more nonzero values? Why are equal values not helpful? Sep 21 '17 at 0:34
• Because the problem's solution is not values used for decision making, but are supposed to be approximates of actual real world values. I can simply observe that in the real world case, the number of non-zero variables are way higher than in the solution of my LP. In other words, what is actually taking place is not the optimal solution.
– Fred
Sep 21 '17 at 21:38
• So perhaps this is something that should not be modeled as an optimization problem at all. Without knowing anything about whatever it is you're trying to model, there's not much else we can say. Sep 24 '17 at 7:58
• There are many examples of usage of optimization for modeling behavior of users of a system, assuming that the users try to minimize their costs, such as the traffic assignment problem, based on Wardrop's equilibrium.
– Fred
Sep 24 '17 at 15:17

There is an non convex optimisation based approach used in signal processing : you can use NMF with sparse constraints. Use a multiplicative gradient descent with l1 regularisation step, and project your inequality constraint at each step.

EDIT: Roughly the idea of Lee, DD & Seung, HS (1999) is to use a "multiplicative" gradient update to maintain nonegativity during the factorisation : split the gradient into two positive components $g=g_+−g_−g$ and update using $g_+/g_−$ direction.

Then you can include constraint inside the gradient descent, that would be the fuzzy part of my proposal. See & Seung algorithm as very little mathematical guarantees but somehow performs very well. There are a lot of contrained variants of NMF, you can check J. Le Roux, J. R. Hershey, and F. Weninger. "Sparse NMF – half-baked or well done ?" for sparse constraints.

• Can you provide a reference? I googled the terms you mentioned, but it is hard to fully understand. Also, is there a Matlab function that do this (or any other language)?
– Fred
Sep 20 '17 at 17:02
• sklearn and matlab have a NMF function, but i guess you should implement the iterations yourself. Roughly the idea of Lee, DD & Seung, HS (1999) is to use a "multiplicative" gradient update : split the gradient into two positive components $g = g_+ - g_-$ and update by $g_+/g_-$. Then you can include constraint inside the gradient descent, that would be the fuzzy part of my proposal. See&Seung algorithm as very little mathematical guarantees but somehow performs very well. There are a lot of contrained variants of NMF, you can check J. Le Roux, J. R. Hershey, and F. Weninger. Sep 21 '17 at 12:02
• Thank you for your help. It sounds a little too complicated for me now. I will try other suggestions first, and if they did not work, I will come back and try to understand this.
– Fred
Sep 21 '17 at 21:35