Traditionally, linear programming is used to find the one optimal solution to a set of constraints, variables and a goal (all described as linear relationships). Sometimes, when the objective is parallel to a constraint, there are infinite or many equally good optimal solutions. I am not asking about this latter case.

I am more interested in finding many solutions that are in the feasible region generated by my set of constraints. But I would like the solutions I find to be 'scattered' around the feasible region in the sense that they are maximally far from one another. Is there a known way to, without running a solver multiple times, generate multiple solutions and use the objective function to enforce that the solutions should be separated?

For example, any linear program with decisions a and b and constraints w <= a <= x and y <= b <= z can be 'duplicated' to find two solutions. Our new linear program has variables a1, a2, b1, and b2 and the constraints w <= a1 <= x and w <= a2 <= x and similar for b1, b2. However, when it comes to forming an objective function we run into trouble as we can not use norms other than the L1-norm without discarding linearity and we can not truly even use the L1 norm because it is not possible (so far as I know) to encode absolute values.

Perhaps I should look into convex optimization or semidefinite programming or something?

Is there a known way to generate a set of solutions to a linear program, and using an objective that enforces "distance" between the solutions?

  • 1
    $\begingroup$ Compute the smallest cube surrounding your feasible region (if it's unbounded, choose some bounded part), lay a hypergrid with the desired resolution over it and discard all points that don't fulfill the restrictions. Would that work for you? $\endgroup$
    – Raphael
    Apr 11, 2013 at 22:47
  • $\begingroup$ That might work for me, although it is not clear to me how I would go about computing the hypercube, and I think the feasible region I am investigating is highly non-trivial - I expect that many points would have to be discarded. My particular application has tens of thousands of variables/decisions and hundreds of constraints. $\endgroup$
    – Ross
    Apr 11, 2013 at 23:00
  • $\begingroup$ A basic feasible solution is in a vertex of the polytope. Can you not look at the normals of the incident faces to compute a direction "across" the polytope and follow it to the boundary of the feasible region? That should give you reasonably different solutions, but probably not the most different ones. $\endgroup$
    – adrianN
    Aug 10, 2013 at 16:22
  • $\begingroup$ Don't use a cube. Use an ellipsoid (specifically, a small ellipsoid containing the polytope, which can be found by the ellipsoid method). That way, you're guaranteed to find a reasonable number of points in the region. $\endgroup$
    – Peter Shor
    Sep 6, 2013 at 23:32

2 Answers 2


A heuristic, using linear programming

One approach might be to pick a random objective function, and maximize it. Then repeat, with a different set of objective functions each time.

In other words, suppose the unknowns are $x_1,x_2,\dots,x_n$, and you have some constraints $\mathcal{C}$. In each iteration you pick $c_1,c_2,\dots,c_n \in \mathbb{R}$ randomly, then search for a solution that maximizes $c_1 x_1 + \dots + c_n x_n$ subject to the constraints $\mathcal{C}$.

I would expect this heuristic might often find a somewhat scattered set of solutions -- not necessarily maximally scattered (maximally far from each other), but probably not too close to each other, either.

Maximizing the average pairwise L2 distance, using quadratic programming

Alternatively, use quadratic programming. For simplicity, let's look at the problem of finding two solutions. Suppose you want two solutions $x,y$ that are as far apart from each other as possible, under the $L_2$ norm (Euclidean distance). Then this can be formulated as a quadratic programming problem.

Basically, you want to maximize the squared distance $d(x,y)^2 = (x_1-y_1)^2+\cdots+(x_n-y_n)^2$ between $x$ and $y$, subject to the requirement that both $x$ and $y$ must satisfy the constraints. This is the problem of maximizing a quadratic objective function, with linear constraints -- i.e., quadratic programming.

If you want $k$ points that are maximally scattered, this is also possible. Say the points are $x^1,\dots,x^k \in \mathbb{R}^n$. Then you could maximize the objective function

$$\sum_{i<j} d(x^i,x^j)^2,$$

i.e., the function

$$\sum_{i<j} \sum_\ell (x^i_\ell - x^j_\ell)^2.$$

This is a quadratic function, and you have linear constraints $\mathcal{C}$ on each of the points $x^i$, so this is a quadratic programming instance. It finds you points that are maximally scattered in the sense that the average pairwise distance is maximized.

You can also formulate a greedy variant of this algorithm, where you already have $k$ solutions, and you want to find a $k+1$th solution that satisfies all the linear inequalities and also maximizes the average L2 distance from it to the other $k$ solutions. That too can be formulated as a quadratic programming problem.

Quadratic programming is harder than linear programming, but there are off-the-self solvers that will solve quadratic programming problems for you.

Maximizing the minimal pairwise L2 distance, using QCQP

Finally, let's say you want your $k$ points to be scattered in the sense that you want to maximize the minimum pairwise distance. In other words, let's say you want to find the largest possible threshold $t$ such that it is possible to find $k$ points $x^1,\dots,x^k \in \mathbb{R}^n$ that each satisfy the linear constraints, and such that each pair of points is at distance $t$ away from each other: $d(x^i,x^j) \ge t$ for all $i<j$. Then this can be formulated as a quadratic optimization program with quadratic constraints, i.e., QCQP. QCQP is even harder, but there are off-the-shelf solvers for QCQP you could try out, too.


I found an approach to generating absolute values.

Suppose we have the variables $a_1$, $a_2$, $b_1$ and $b_2$ and a bunch of constraints. Our objective functions looks something like: maximize $|a_1 - a_2| + |b_1 - b_2|$; the idea being we want to maximize the L1 norm of these two solutions (as per the original question).

We can introduce "slack variables" abs_a and abs_b and the constraints:

$$\mathrm{abs_a} + a_1 - a_2 \leq 0$$

$$\mathrm{abs_a} - a_1 + a_2 \leq 0$$

and similarly for $b_1$ and $b_2$. These constraints force $\mathrm{abs_a}$ to be at most the difference between $a_1$ and $a_2$, and possibly less. In other words $\mathrm{abs_a}$ can not be larger than the maximum difference between $a_1$ and $a_2$.

Then what is left is replace the objective function: maximize $\mathrm{abs_a} + \mathrm{abs_b}$.

  • $\begingroup$ Actually, these encodings only work for the minimization of absolute values. Therein it doesn't solve my problem. More information here: lpsolve.sourceforge.net/5.5/absolute.htm $\endgroup$
    – Ross
    Apr 12, 2013 at 18:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.