2

I have a large linear program, which I need to solve fast. I'm willing to sacrifice optimality but not feasibility. I have an initial feasible point x0. Any suggestions for an algorithm?

For example, I could use the simplex algorithm (using my x0 as warm-start, although it's not necessarily bfs/vertex) to jump from vertex to vertex and stop it after 10 iterations or when the time is up.

Another example is to use a barrier method that improves the objective while staying in the feasible domain.

Again, the point is to start with a feasible solution and improve it in the time limit.

Zohar Levi
  • 253

2 Answers2

2

The best approach approach to a limited approach is to collect as many extreme points and extreme directions via the Simplex method as possible through a finite amount of iterations, then given those finite amount of extreme points and directions, we can construct other feasible points with the following definition:

We can represent $X = (x_0, \ldots, x_n)^$ as a convex combination of the extreme points of model plus a conical combination of the extreme directions of a model.

The reason for this is further explained here. However, with a collection of extreme points and directions, we can find other feasible points within the model as such:

$$\begin{bmatrix} x_{00} & \cdots & x_{0n}\\ \vdots & \ddots & \vdots\\ x_{m0} & \cdots & x_{mn} \end{bmatrix} \begin{bmatrix} \lambda_0 \\ \vdots \\ \lambda_n\end{bmatrix} + \begin{bmatrix} d_{00} & \cdots & d_{0n}\\ \vdots & \ddots & \vdots\\ d_{m0} & \cdots & d_{mn} \end{bmatrix} \begin{bmatrix} l_0 \\ \vdots \\ l_m \end{bmatrix} = \begin{bmatrix} X_0 \\ \vdots \\ X_m\end{bmatrix} $$

With $\lambda_m,l_m\in\Bbb R^+$ $\forall m$, and $x$ represents the column/row of each chosen extreme point, $d$ represents the column/row of each chosen extreme direction, and $X$ represents the values that correspond to a new feasible point in the model. An example of this in action has been done here on another LP problem.

Of course, this isn’t the only way to find interior points of a model, and thus if above isn’t your fancy, there’s the option of other interior-point methods approaches (and more barrier methods). However, collecting as much information from Simplex is your best call before moving forward with most other approaches.

1

The primal-dual algorithm may be a useful approach in this scenario. A one-sentence summary: starting from an initial feasible solution, we find an optimal direction to go in to increase the objective function until we run into a boundary, then find a new direction to go in.

Compared to the simplex method, this has more complicated iterations (each one will need several steps of the simplex method) but each iteration can yield a greater improvement. As with many other methods, we can stop early to get a feasible but not optimal solution. I think it's a good choice in this setting because:

  • It makes more progress at the beginning: early iterations are faster and often produce more improvement than later ones, whereas the simplex method is more uniform.
  • It takes better advantage of the initial solution, which is feasible but not a corner point.

Here's a slightly more detailed description.

To be consistent with notation other people use, we're going to think of the problem you want to solve as the "dual program", and assume it is in the form $$(D): \quad \max\{\mathbf u^{\textsf T}\mathbf b : \mathbf u^{\textsf T}A \le \mathbf c^{\textsf T}\}$$ with $\mathbf u_0$ as the initial feasible solution (and $\mathbf u_1, \mathbf u_2, \dots$ as later feasible solutions).

Our goal at each iteration is to find a direction $\mathbf v$ such that for not-too-large $t > 0$, $\mathbf u_i + t \mathbf v$ is still feasible, and increases the objective function quickly. Then we will find the largest $t^*$ that preserves feasibility, and let $\mathbf u_{i+1} = \mathbf u_i + t^* \mathbf v$.

To find $\mathbf v$, we write down a smaller linear program called the "restricted dual". For this, let $J$ be the set of dual constraints that are tight at $\mathbf u_i$ (that is, $(\mathbf u_i)^{\textsf T}A_J = \mathbf c_J^{\textsf T}$). We find $\mathbf v$ by solving $$ (RD) : \quad \max \{ \mathbf v^{\textsf T}\mathbf b : \mathbf v^{\textsf T}A_J \le \mathbf 0^{\textsf T}, \mathbf v \le \mathbf 1\} $$ where the constraint $\mathbf v \le \mathbf 1$ is added because we're looking for a direction, and don't care about magnitude.

The reason that these are all duals is that the dual of $(RD)$, the "restricted primal" $(RP)$, is the one we actually solve to find $\mathbf v$ (using the simplex method). This is $$ (RP) : \quad \min\{\mathbf 1^{\mathsf T}\mathbf y : A_J \mathbf x_J + I \mathbf y = \mathbf b, \mathbf x_J \ge \mathbf 0, \mathbf y \ge \mathbf 0\}. $$ (Note that from the optimal dictionary/tableau for $(RP)$, we can find the optimal $\mathbf v$ using the reduced costs of $\mathbf y$, since the $\mathbf y$-variables are dual to the $\mathbf v \le \mathbf 1$ constraints.)

The reason to use $(RP)$ instead of $(RD)$ is that between different iterations, the previous optimal solution $(RP)$ can be taken as a starting point for the next $(RP)$: it is still feasible, even though the set of variables changes. Intuitively, $(RP)$ can be thought of as a version of a bigger problem with constraints $A\mathbf x + I\mathbf y = \mathbf b$, but most of its variables (the ones outside $J$) are "frozen" and cannot be made positive.

For more resources on this topic: I learned about the primal-dual algorithm from the book Combinatorial Optimization by Papadimitriou and Steiglitz (it's in Chapter 5). I have some old lecture notes online about the algorithm as well, which can be found here; the direct links to three lectures' worth of notes on this topic are 1, 2, 3.

Misha Lavrov
  • 142,276