Parallel parametric linear programming solving, and application to polyhedral computations

TL;DR

This paper introduces a task-based parallel scheme for parametric linear programming, significantly improving computational speed for large-scale polyhedral problems and control applications.

Contribution

It presents a novel parallelization method for parametric linear programming that achieves near-linear speedup, enhancing efficiency in polyhedral computations.

Findings

Achieves quasi-linear speedup on large problems

Effective parallelization of parametric linear programming

Improves computational efficiency in control applications

Abstract

Parametric linear programming is central in polyhedral computations and in certain control applications.We propose a task-based scheme for parallelizing it, with quasi-linear speedup over large problems.

Equations

$$o=(3{\mu }_1+3{\mu }_2)+\left(-\frac{3}{8}{\mu }_1-\frac{1}{8}{\mu }_2\right)x_3+\left(-\frac{1}{8}{\mu }_1-\frac{3}{8}{\mu }_2\right)x_4$$

Full text

1 Introduction

A convex polyhedron, or polyhedron for short here, in dimension $n$ is the solution set over $\mathbb{Q}^{n}$ (or, equivalently, $\mathbb{R}^{n}$) of a system of inequalities (with integer or rational coefficients). Polyhedra in higher dimension are typically used to enclose the reachable states of systems whose state can be expressed, at least partially, as a vector of reals or rationals; e.g. hybrid systems or software [3].

The conventional approaches for polyhedral computations are the dual description (using both vertices and faces) and Fourier-Motzkin elimination. They both suffer from high complexity on relevant cases. We instead express image, projection, convex hull etc. as solutions to parametric linear programmings, where parameters occur linearly within the objective. A solution to such a program is a quasi-partition of the space of parameters into polyhedra, with one optimum associated to each polyhedron. The issue is how to compute this solution efficiently. In this article, we describe how we parallelized our algorithm.

2 Sequential algorithms

Here we are leaving out how polyhedral computations such as projection and convex hull can be reduced to parametric linear programming — this is covered in the literature [4, 7] — and focus on solving the parametric linear programs.

2.1 Non-parametric linear programming (LP)

A linear program with $n$ unknowns is defined by a system of equations $AX=B$, where $A$ is an $m\times n$ matrix; a solution is a vector $X$ such that $X\geq 0$ on all coordinates and $AX=B$. The program is said to be feasible if it has at least one solution, infeasible otherwise. In a non-parametric linear program one considers an objective $C$: one wants the solution that maximizes $C^{T}X$. The program is deemed unbounded if it is feasible yet it has no such optimal solution.

Example 1*.*

Consider the polygon $P$ defined by $x_{1}\geq 0$, $x_{2}\geq 0$, $3x_{1}-x_{2}\leq 6$, $-x_{1}+3x_{2}\leq 6$. Define $x_{3}=6-3x_{1}+x_{2}$ and $x_{4}=6+x_{1}-3x_{2}$. Let $X=(x_{1},x_{2},x_{3},x_{4})$, and then $P$ is the projection onto the first two coordinates of the solution set of $AX=B\land X\geq 0$ where $A=\left[\begin{smallmatrix}1&-3&0&-1\\ -3&1&-1&0\end{smallmatrix}\right]$ and $B=\left[\begin{smallmatrix}6\\ 6\end{smallmatrix}\right]$.

An LP solver takes as input $(A,B,C)$ and outputs “infeasible”, “unbounded” or an optimal solution. Most solvers work with floating-point numbers and their final answer may be incorrect: they may answer “infeasible” whereas the problem is feasible, or give “optimal solutions” that are not solutions, or not optimal.

In addition to a solution $X^{*}$, solvers also provide the associated basis: $X^{*}$ is defined by setting $n-m$ of its coordinates to [math] (known as nonbasic variables) and solving for the other coordinates (known as basic variables) using $AX^{*}=B$, and the solver provides the partition into basic and nonbasic variables it used. If a floating-point solver is used, it is possible to reconstruct an exact rational point $X^{*}$ using that information and a library for solving linear systems in rational arithmetic. One then checks whether it is truly a solution by checking $X^{*}\geq 0$.

The optimal basis also contains a proof of optimality of the solution. We compute the objective function $C^{T}X$ as $\sum_{i\in N}\alpha_{i}X_{i}+c$ where $N$ is the set of indices of the nonbasic variables and $c$ is a constant, and conclude that the solution obtained by setting these nonbasic variables to [math] is maximal because all the $\alpha_{i}$ are nonpositive. If $X^{*}$ is not a solution of the problem ($X^{*}\geq 0$ fails) or is not optimal, then we fall back to an exact implementation of the simplex algorithm.

Example 1* (continued).*

Assume the objective is $C=\begin{bmatrix}1&1&0&0\end{bmatrix}$, that is, $C^{T}X=x_{1}+x_{2}$. From $AX=B$ we deduce $x_{1}=3-\tfrac{3}{8}x_{3}-\tfrac{1}{8}x_{4}$ and $x_{2}=3-\tfrac{1}{8}x_{3}-\tfrac{3}{8}x_{4}$. Thus $x_{1}+x_{2}=6-\frac{1}{2}x_{3}-\frac{1}{2}x_{4}$.

Assume $x_{3}$ and $x_{4}$ are nonbasic variables and thus set to [math], then $X^{*}=(x_{1},x_{2},x_{3},x_{4})=(3,3,0,0)$. It is impossible to improve upon this solution: as $X\geq 0$, changing the values of $x_{3}$ and $x_{4}$ can only decrease the objective $o=6-\frac{1}{2}x_{3}-\frac{1}{2}x_{4}$. This expression of $o$ from the nonbasic variables can be obtained by linear algebra once the partition into basic and nonbasic variables is known.

While the optimal value $C^{T}X^{*}$, if it exists, is unique for a given $(A,B,C)$, there may exist several $X^{*}$ for it, a situation known as dual degeneracy. The same $X^{*}$ may be described by different bases, a situation known as primal degeneracy, happening when more than $n-m$ coordinates of $X^{*}$ are zero, and thus some basic variables could be used as nonbasic and the converse.

2.2 Parametric linear programming (PLP)

For a parametric linear program, we replace the constant vector $C$ by $C_{0}+\sum_{i=1}^{k}\mu_{i}C_{i}$ where the $\mu_{i}$ are parameters.111There exists another flavor of PLP with parameters in the right-hand sides of the constraints. When the $\mu_{i}$ change, the optimum $X^{*}$ changes. Assume temporarily that there is no degeneracy. Then, for given values of the $\mu_{i}$, the problem is either unbounded, or there is one single optimal solution $X^{*}$. It can be shown that the region of the $(\mu_{1},\dots,\mu_{k})$ associated to a given optimum $X^{*}$ is a convex polyhedron (if $C_{0}=0$, a convex polyhedral cone), and that these regions form a quasi partition of the space of parameters (two reegions may overlap at their boundary, but not in their interior) [4, 5, 7]. The output of the parametric linear programming solver is this quasi-partition, and the associated optima—in our applications, the problem is always bounded in the optimization directions, so we do not deal with the unbounded case.

Let us see in more details about how to compute these regions. We wish to attach to each basis (at least, each basis that is optimal for at least one vector of parameters) the region of parameters for which it is optimal.

Example 1* (continued).*

Instead of $C=\begin{bmatrix}1&1&0&0\end{bmatrix}$ we consider $C=\begin{bmatrix}\mu_{1}&\mu_{2}&0&0\end{bmatrix}$. Let us now express $o=C^{T}X$ as a function of the nonbasic variables $x_{3}$ and $x_{4}$:

[TABLE]

The coefficients of $x_{3}$ and $x_{4}$ are nonpositive if and only if $3\mu_{1}+\mu_{2}\geq 0$ and $\mu_{1}+3\mu_{2}\geq 0$, which define the cone of optimality associated to that basis and to the optimum $X^{*}=(3,3,0,0)$.

The description of the optimality polyhedron by the constraints obtained from the sign conditions in the objective function may be redundant: containing constraints that can be removed without changing the polyhedron. Our procedure [6] for removing redundant constraints from the description of a region $R_{1}$ also provides a set of vectors outside of $R_{1}$, a feature that will be useful.

Assume now we have solved the optimization problem for a vector of parameters $D_{1}$, and obtained a region $R_{1}$ in the parameters (of course, $D_{1}\in R_{1}$). We store the set of vectors outside of $R_{1}$ provided by the redundancy elimination procedure into a “working set” $W$ to be processed, choose $D_{2}$ in it. We compute the region $R_{2}$ associated to $D_{2}$. Assume that $R_{2}$ and $R_{1}$ are adjacent, meaning that they have a common boundary. We get vectors outside of $R_{2}$ and add them to $W$. We pick $D_{3}$ in $W$, check that it is not covered by $R_{1}$ or $R_{2}$, and, if it is not, compute $R_{3}$, etc. The algorithm terminates when $W$ becomes empty, meaning the $R_{1},\dots$ produced form the sought quasi-partition.

This simplistic algorithm can fail to work because it assumes that it is discovering the adjacency relation of the graph. The problem is that, if we move from a region $R_{i}$ to a vector $D_{j}\notin R_{i}$, it is not certain that the region $R_{j}$ generated from $D_{j}$ is adjacent to $R_{i}$ — we could miss some intermediate region. We modify our traversal algorithm as follows. The working set contains pairs $(R,D^{\prime})$ where $R$ is a region and $D^{\prime}\notin R$ a vector (there is a special value none for $R$). The region $R^{\prime}$ corresponding to $D^{\prime}$ is computed. If $R$ and $R^{\prime}$ are not adjacent, then a vector $D^{\prime\prime}$ in between $R$ and $R^{\prime}$ is computed, and $(R,D^{\prime\prime})$ added to the working set. This ensures that we obtain a quasi-partition in the end. Additionally, we obtain a spanning tree of the region graph, with edges from $R$ to $R^{\prime}$.

The last difficulty is degeneracy. We have so far assumed that each optimization direction corresponds to exactly one basis. In general this is not the case, and the interiors of the optimality regions may overlap. This hinders performance. The final result is no longer a quasi-partition, but instead just a covering of the parameter space—enough for projection, convex hull etc. being correct.

3 Parallel parametric linear programming

Our algorithms are designed in a task-based execution model. The sequential algorithm executes tasks taken from a working set, which can themselves spawn new tasks. In addition, it maintains the set $\mathit{regions}$ of regions already seen, used:

i) for checking if a vector $D$ belongs to a region already covered ($\mathit{is\_covered}$);

ii) for checking adjacency of regions;

iii) for adding new regions found.

Therefore, in a parallel task model, this algorithm is straightforwardly parallel. The regions are inserted into a concurrent array. We investigated two task scheduling strategies. A static approach starts all the available tasks, waits for them to complete and collects all the new tasks $(R,D)$ into the working set, until no new task is created and the working set is empty. A dynamic approach allows adding new tasks to the working set dynamically and runs the tasks until that set is empty.

The number of tasks running to completion (not aborted early due to a test) is the same as the number of generated regions. The $\mathit{is\_covered}(D,\mathit{regions})$ loop can be easily parallelized as well. We opted against it as it would introduce a difficult-to-tune second level of parallelism.

We implemented these algorithms using Intel’s Thread Building Blocks (TBB [8]) and OpenMP tasks [1], both providing a task-based parallelism model with different features.

The dynamic task queue can be implemented using TBB’s tbb::parallel_do, which dynamically schedules tasks from the working set on a number of threads. The static scheduling approach can simply be implemented by a task synchronization barrier (such as OpenMP’s barrier).

That first implementation of the dynamic task scheduling approach was slow. The working set often contained tasks such that the regions generated from them were the same, leading to redundant computations. The workaround was to add a hash table storing the set of bases (each being identified by the ordered set of its basic variables) that have been or are currently being processed. A task will abort after solving the floating-point linear program if it finds that its basis is already in the table.

4 Performance evaluation

We implemented our parallel algorithms in C++, with three alternate schemes selectable at compile-time: no parallelism, OpenMP parallelism or TBB.

All benchmarks were run on the Paranoia cluster of Grid’5000 [2] and on a server called Pressembois. Paranoia has 8 nodes, each with 2 Intel® Xeon® E5-2660v2 CPUs (10 cores, 20 threads/CPU) and 128 GiB of RAM. Code was compiled using GCC 6.3.1 and OpenMP 4.5 (201511). The nodes run Linux Debian Stretch with a 4.9.0 kernel. Pressembois has 2 Intel Xeon Gold 6138 CPU (20 cores/CPU, 40 threads/CPU) and 192 GiB of RAM. It runs a 4.9 Linux kernel, and we used GCC 6.3. Every experiment was run 10 times. The plots presented in this section provide the average and standard deviation. Paranoia was used for the OpenMP experiments, whereas Pressembois was used for TBB.

We evaluated our parallel parametric linear programming implementation by using it to project polyhedra, a very fundamental operation. We used a set of typical polyhedra, with different characteristics: numbers of dimensions, of dimensions to be projected and of constraints, sparsity. Here we present a subset of these benchmarks, each comprising 50 to 100 polyhedra.

On problems that have only few regions, not enough parallelism can be extracted to exploit all the cores of the machine. For instance, Figure 1 presents two experiments on 2 to 36 regions using the OpenMP version. It gives an acceptable speed-up on a few cores (up to 10), then the computation does not generate enough tasks to keep the additional cores busy. As expected, when the solution has many regions, computation scales better. Figure 2 presents the performance obtained on polyhedra made of 24 constraints, involving 8 to 764 regions, using the OpenMP version. The speed-up is sublinear, especially beyond 20 cores.

On larger polyhedra, with 120 constraints and 50 variables, the speedup is close to linear with both OpenMP and TBB (Fig. 3). The parallelism extracted from the computation is illustrated by Fig. 4, on a polyhedron involving 29 constraints and 16 variables. Figure 4b shows the number of parallel tasks.