Capital letters (e.g. A) denote matrices, small bold letters (e.g. x x x) denote vectors, and small letters (e.g. b) denote scalars. The ith row of A is a a a i• , its jth column is a a a • j . For simplification, we use a a a i to abridge a a a i• . The element at the ith row and the jth column of the matrix A is denoted by a i j . The superscript T (e.g. A T ) represents the transpose of a matrix or a vector. P denotes a polyhedron and C a constraint. The ith constraint of P is C i . Q denotes the field of rational numbers, and F is the set of finite floating-point numbers, which is a subset of Q. In the pseudo-code the floating point (or rational) variables are denoted by name F (or name Q ); name Q×F means that the variables are stored in both rational and floating point numbers.

Contributions

 in several aspects and provided a more efficient implementation.

. We provide checkers to guarantee the soundness of the results, and thus the floating-point solver will always return a sound solution. When it is necessary, a rational solver is invoked which guarantees the precision of our result.

Degeneracy We resolve some difficulties due to geometric degeneracy, which previously resulted in redundant computations. We aim at obtaining a partition of the space without overlapping areas because of two reasons: i) we want to avoid the redundant computations; ii) when there is no overlaps, we are able to check the precision of our result. But the occurrence of degeneracy leads to overlaps, and thus we need to deal with the degeneracy. We proved that there is no dual degeneracy in our PLP solver, and we proposed a method to avoid the primal degeneracy. Then we can check the precision of our result. If the result is not precise, i.e., there are missing solutions of PLP, a rational solver will be launched for computing these missing solutions.

. The details will be explained in Section 5.3.

In 2012 Verimag started implementing a library using constraint-only description, which is called VPL (Verified Polyhedra Library)2 [START_REF] Fouilhé | Revisiting the abstract domain of polyhedra: constraints-only representation and formal proof[END_REF][START_REF] Maréchal | New Algorithmics for Polyhedral Calculus via Parametric Linear Programming[END_REF]. There are several reasons for using only constraints. We have already cited the high complexity of the generator description of some polyhedra commonly found in abstract interpretation. Besides, the high cost of Chernikova's algorithm is able to be avoided if we use the single description. Another reason is to be able to certify the results of the computation, in particular that the obtained polyhedra includes the one that should have been computed, i.e., the soundness of the result.

In the first version of VPL, which was implemented by Alexis Fouilhé in 2015 [START_REF] Fouilhé | Revisiting the abstract domain of polyhedra: constraints-only representation and formal proof[END_REF], all main operations boil down to projection, which is performed using Fourier-Motzkin elimination [START_REF] George B Dantzig | Fourier-motzkin elimination and its dual[END_REF]. This method generates many redundant constraints which must be eliminated at a high cost. Besides, for projecting out multiple variables, it computes all intermediate steps. For each variable λ i to be eliminated, the inequalities are divided into three groups depending on the sign (positive, negative and zero) of the coefficient of λ i . Then it computes the combination of each inequality from the positive group and that from the negative group. These intermediate steps may have high description complexity. Jones et al. explained the relationship between polyhedral projection and parametric linear programming (PLP) in [START_REF] Jones | On polyhedral projection and parametric programming[END_REF]. Based on this work, [START_REF] Jacob | Polyhedral analysis using parametric objectives[END_REF] proposed a new approach to compute polyhedral projection using PLP, which can replace Fourier-Motzkin elimination. Inspired by [START_REF] Jones | On polyhedral projection and parametric programming[END_REF] and [START_REF] Jacob | Polyhedral analysis using parametric objectives[END_REF], Alexandre Maréchal implemented a novel algorithm for computing the projection and convex hull, which both boil down to PLP, and released version 2.0 of the VPL. VPL 2.0 was implemented in arbitrary precision arithmetic in OCaml [START_REF] Maréchal | Scalable minimizing-operators on polyhedra via parametric linear programming[END_REF]. The algorithm of the PLP solver is intrinsically parallelable, however it is developed with OCaml, which does not well support parallelism programming. Our work focuses on improving the efficiency of the algorithm of PLP solver, and parallelizing it for obtaining a better performance.

In recent years some state-of-the-art methods are presented that aim at improving the efficiency of polyhedra domain. The expensive cost of polyhedra domain is partially caused by the "dimension curse" of double description. Thus there appears to be different ideas to solve the problem. What we did was getting rid of the double description and developing new algorithms for the operators that are inefficient using constraint-only description. There are other methods which focus on reducing the dimension of the polyhedra [START_REF] Halbwachs | Some ways to reduce the space dimension in polyhedra computations[END_REF][START_REF] Singh | Fast polyhedra abstract domain[END_REF]. In the polyhedra library ELINA 3 [START_REF] Singh | Fast polyhedra abstract domain[END_REF] the variables are divided into small groups, where the variables in different groups are not related to each other. The authors compared the performance of ELINA with the NewPolka library in Apron, and PPL library. By decreasing the dimension of polyhedra, i.e, the number of variables, the ELINA library has a speed up of dozens of times to thousands of times.

In other aspects, most of the existing polyhedra libraries use rational arithmetic, and thus some new attempts use floating point arithmetic to improve the efficiency. Apart from our approach, the work in [START_REF] Chen | A sound floating-point polyhedra abstract domain[END_REF] also presented an polyhedra domain in floating point numbers. The authors used constraint-only description, and implemented a sound approach based on Fourier-Motzkin elimination. As the floating point arithmetic is not precise, they obtain a set of linear interval inequalities after eliminating one variable. Then they apply a sound linearization approach to these inequalities. Comparing with their method, our approach has the advantage in several aspects: i) our approach gets rid of the intermediate steps of Fourier-Motzkin elimination, which produce amount of redundant constraints; ii) we can project out multiple variables each time; iii) our results are more precise as we reconstruct the rational solutions.

Chapter 2 Preliminaries

In this chapter, we will present some important preliminaries including convex polyhedra, floating point arithmetic, non-parametric and parametric linear programming, etc.

The polyhedra abstract domain is firstly presented by P. Cousot and N. Halbwachs in [START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF]. A convex polyhedron is composed of a set of linear constraints (equalities or inequalities). In static analysis, a linear constraint represents the relationship between the variables of a program. Compared with other numerical abstract domains, polyhedra domain is the most precise but most expensive one. There are dozens of polyhedral operators in the polyhedra domain. In this paper we will mainly talk about minimization, projection and convex hull.

A floating point number is the representation of a rational number in the computers, and the representation may be approximate because of the rounding errors. When we use floating point arithmetic, it is necessary to measure the round errors to guarantee the correctness of the results. The study of the consequence of the rounding errors is a topic of numerical analysis. In this paper, we used some basic knowledge of numerical analysis, which is explained in [START_REF] Nicholas | Accuracy and stability of numerical algorithms[END_REF].

Linear programming is a method to compute the optimal value of a linear objective function subjecting to a set of linear constraints. In the past decades several methods to solve linear programming problems have been presented [START_REF] George B Dantzig | Maximization of a linear function of variables subject to linear inequalities[END_REF][START_REF] George | Linear programming 2: theory and extensions[END_REF]. In this paper we will only use Dantzig's Simplex algorithm. The difference between non-parametric and parametric linear programming is that the objective function of the latter contains parameters. By setting the parameters to various values we obtain a set of nonparametric linear programs whose objective functions are different. The parametric linear program can be solved by computing these non-parametric linear programs and generalizing the results. We will also present the method to solve a PLP problem using the Simplex algorithm directly. Both methods will be used in the following chapters.

Convex polyhedra and polyhedra domain

As we aim to implement an algorithm to manipulate the operations on convex polyhedra, in this section we firstly present the definition and the descriptions of a convex polyhedron. We also introduce some operators in the polyhedra domain.

Convex polyhedra

A polyhedron is made of faces, edges and vertices. A polyhedron is bounded if it does not contain unbounded faces. All the polyhedra we will mention in this paper are convex polyhedra, and thus we will omit the word "convex" in the following context. According to [START_REF] Grünbaum | Convex polytopes[END_REF], a convex set is defined as in Definition 2.1.

Double description of polyhedra

A polyhedron can be represented by constraints or generators. A constraint defines a linear relation between variables. It can be an equality, a strict or non-strict inequality. A constraint C i on variables λ λ λ in Q n is represented by Equation 2.1.

a a a T i λ λ λ ≤ b i (2.1)
where a a a i = [a i0 , • • • , a in ] T are the coefficients of constraint C i and b i is a constant. In the equation ≤ can also be <, >, = or ≥, and we use ≤ as an example here.

Equation 2.2 shows a polyhedron defined by the conjunction of m constraints.

P = {λ λ λ ∈ Q n | m i=1 a a a T i λ λ λ ≤ b i } (2.2)
The Equation 2.2 can be represented in the matrix form. Give A ∈ Q m×n and b b b ∈ Q m , the same polyhedron can be represented by Equation 2.3.

P = {λ λ λ ∈ Q n |Aλ λ λ ≤ b b b} (2.3)
A polyhedron can also be represented by generators, i.e. vertices, rays, and lines, as it is shown in Equation 2. [START_REF]OpenMP Application Programming Interface[END_REF]. A bounded polyhedron is represented by a set of vertices, and an unbounded polyhedron is represented by the combination of vertices, rays and lines. The rays and lines define the unbounded directions. A line can be considered as two rays which are collinear and in opposite directions.

P = {λ λ λ ∈ Q n |λ λ λ = |V| ∑ i=1 α i v v v i + |Y| ∑ i=1 β i r r r i + |L| ∑ i=1 γ i l l l i , α i , β i ≥ 0, |V| ∑ i=1 α i = 1} (2.4)
where V is a set of vertices, Y is a set of rays, L is a set of lines, and α i , β i , γ i are constants.

Example 2.2. We can represent the unbounded polyhedron in Figure 2.2 by constraints as:

P = {2λ 2 -1 ≥ 0, λ 1 -λ 2 ≥ 0}
The corresponding generator description is:

V = {( 1 2 , 1 2
)} Y = {(1, 0), (1, 1)}

Operators in polyhedra domain

The polyhedra domain contains dozens of operators. We will present the most frequently used ones.

Elimination of redundancy A polyhedron may contain constraints that are redundant, meaning that removing these constraints will not change the polyhedron. The process of removing redundant constraints is called minimization. A polyhedron whose redundant constraints are removed is called a minimized polyhedron. We will present this operator in Section 3.2.

Meet

The meet of polyhedra P 1 and P 2 , which is denoted by P 1 ⊓ P 2 , is the conjunction of the constraints of these two polyhedra, and then the result needs to be minimized. Join The join operator computes the convex hull of polyhedra. The convex hull of P 1 and P 2 , i.e. P 1 ⊔ P 2 , is the smallest convex set that contains P 1 and P 2 . This operator will be introduced in Section 3.3.

Inclusion

The polyhedron P 1 is included in P 2 if ∀x x x ∈ P 1 , x x x ∈ P 2 . The inclusion test can be implemented by testing redundancy: P 1 ⊑ P 2 if all the constraints of P 2 are redundant with respect to P 1 .

Equality The polyhedra P 1 = P 2 if they define the same convex set, i.e. the two polyhedra contain equivalent constraints. The equality test can be implemented by inclusion test: P 1 = P 2 if P 1 ⊑ P 2 and P 2 ⊑ P 1 .

Projection The projection operator is used to eliminate variables and project the polyhedron to lower dimension. It will be presented in Section 3.3.

Widening The widening operator is used for the loops. It is used to enforce the convergence of iteration sequences or accelerate converging steps of the iterations.

Floating point arithmetic

As floating point arithmetic accounts for a large proportion in our computation, we will illustrate the fundamental of floating point numbers. We start from the representation of a floating point number, and then introduce the rounding errors. We also present some basic knowledge of numerical analysis [START_REF] Nicholas | Accuracy and stability of numerical algorithms[END_REF] which is used to measure the rounding errors.

Floating point numbers

In computers real numbers are represented approximately by floating-point numbers in the form dd • • • d × β e , where dd • • • d is called significand, e is exponent, and β is the base. To unify the form of floating point representation, IEEE (Institute of Electrical and Electronics Engineers) [START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF] defined a technical standard for floating point numbers in 1990s, which is called IEEE Standard for Floating-Point Arithmetic (IEEE 754). There are three specifications in IEEE 754:

-the base β , which is always 2 ;

-the precision, i.e., the number of digits in dd • • • d -an exponent range from e min to e max , with e min = 1e max

According to IEEE 754 standard, the format of a floating point number is shown in Table 2.1. The sign bit is either 0 for positive or 1 for negative. The significand of a normalized number is always in the range of 1 to 2, i.e. in the form 1.d • • • d, and thus there is no need to represent the leading digit 1 explicitly. So the significand 1.d • • • d is stored as d . . . d. A bias needs to be added to the exponent for representing negative numbers. The infinity +∞ and -∞ are denoted by setting the exponent bits all ones and the significand all zeros. The numbers in form 0.d • • • d, which are called denormalized numbers, are represented by setting the exponent all zeros. total sign exponent significand single precision 32 bits 1 bit 8 bits 24 bits (one bit is implicit) double precision 64 bits 1 bit 11 bits 53 bits (one bit is implicit) Table 2.1 -floating point numbers format in IEEE 754 standard Example 2.3. We use the subscript to represent the base, which is either 2 or 10. Let us see how to store the decimal number 0.40625 10 as single precision binary floating point number. 0.40625 10 = 0.25 10 + 0.125 10 + 0.03125 10 = 2 -2 10 + 2 -3 10 + 2 -5 10 = 0.01101 2 = 1.101 2 × 2 -2 . The sign bit stores 1, as it is positive. There are 8 bits for exponent, and the range is -127 10 to 128 10 , so we need to add a bias 127 10 to the exponent. Thus we have -2 10 + 127 10 = 125 10 = 01111101 2 . For the significand, we only store 101. Then we obtain the binary number, which is shown in Table 2.2.

Rounding errors

Floating point numbers are not precise. The most common reason is that some real numbers cannot be represented by limited binary numbers. As floating point numbers are stored as finite binary numbers, the infinite binary numbers or the numbers with too many digits will be rounded. Therefore, there is a difference between the floating point number and the exact value, which is called rounding errors. The assert will fail, because the decimal number 0.1 cannot be finitely represented by a binary number: 0.1 10 = 0.00011 2 which means the digits 0011 would repeat forever in an infinite representation.

Then the binary number will be truncated to limit bits. The stored value could be less or greater than 0.1 depending on the number of significant bits.

Machine epsilon

Because of the rounding errors, we cannot compare two floating-point numbers directly. Instead, we need to use a threshold. To find out an appropriate threshold, it is necessary to measure the error between the floating-point number and the real number it represents. The measurement of the quantity of rounding errors is a problem in the field of numerical analysis. The rounding errors can be measured by the machine epsilon, whose definition is given by Definition 2.2 according to [START_REF] Nicholas | Accuracy and stability of numerical algorithms[END_REF].

Definition 2.2 (Machine epsilon). The machine epsilon ε is the difference between 1.0 and the next representable floating point number greater than 1.0.

Given the precision p, according to IEEE 754 standard, we store p -1 digits in the significand part, and the smallest representable number is obtained by filling the significand bits with 00 . . . 01, which represents 2

Basic arithmetic model

In the basic arithmetic model, the rounding errors between the floating point and the exact computation of an operator op is given by Equation 2.5 [START_REF] Nicholas | Accuracy and stability of numerical algorithms[END_REF].

f l(x op y) = x op y 1 + δ op ∈ {+, -, ×, ÷} (2.5)
where |δ | < u, u = 1 2 ε and ε is the machine epsilon. From Equation 2.5 we have:

| f l(x op y) -(x op y)| < ε| f l(x op y)| (2.6)
where | • | represents the absolute value.

Error bounds of inner product

The only rounding errors we need to measure in our algorithm is that of inner product of floating-point numbers. According to [START_REF] Nicholas | Accuracy and stability of numerical algorithms[END_REF], the error between the inner product of two vectors a a a T r r r and its floating-point value f l(a a a T r r r) is:

| f l(a a a T r r r) -a a a T r r r| ≤ | f l(a a a)| T | f l(r r r)| αu 1 -αu (2.7)
where α = ⌈log 2 m⌉ + 1, m is the dimension of a a a and r r r, and u = 1 2 ε. Furthermore, Equation 2.8 can be derived from Equation 2.7.

| f l(a a a T r r r) -a a a T r r r| ≤ 1.01mu| f l(a a a)| T | f l(r r r)| (2.8)
where m is the size of the vector. As u < ε, we can use ε to replace u for simplifying implementation. To simplify the explanation, we use t(a a a, r r r) to represent the rounding error of the inner product of a a a and r r r.

t(a a a, r r r

) = 1.01mε| f l(a a a)| T | f l(r r r)| (2.9)
In most part of our procedure we do not consider the rounding errors during computation since we reconstruct the final result in rational numbers and verify that it is a solution of the initial rational problem. In several steps we need to consider and measure the rounding to guarantee the correctness of the floating point arithmetic.

Linear programming

The linear programming (LP) is the core of our algorithm. In order to develop an efficient parametric linear programming solver, we need to solve a set of linear programming problems. Therefore, we will present the form of linear programming in this section.

Definition of an LP problem

Definition 2.3 (Linear function). A linear function f (x x x) is a function that satisfies additivity and homogeneity, i.e., f (x x x + y y y) = f (x x x) + f (y y y) and f (x x x) = α f (x x x), where x x x are variables. A linear equation is an equation in the form ∑ n i=0 a i x i + b = 0, where a i and b are constants.

Given a function f (x x x) = ∑ n i=0 a i x i + b, when b ̸ = 0, f (x x x) does not satisfy Definition 2.3, and it is called linear affine function. But sometimes, as what we will do in this paper, a linear affine function is simply called linear function.

Linear programming (LP) is a method to find the optimal value of a linear function subject to a set of linear constraints. Given variables

λ λ λ = [λ 1 , . . . , λ n ] T , an LP problem in standard form is expressed in Problem 2.10. maximize Z(λ λ λ ) = n ∑ j=1 o j λ j subject to n ∑ j=1 a i j λ j = b i ∀i ∈ {1, • • • m} and λ λ λ ≥ 0 (2.10) where Z(λ λ λ ) is the linear objective function, ∑ n j=1 a i j λ j = b i (∀i ∈ {1, • • • m})
is the set of linear constraints, and a i j and o j are coefficients.

Standard form of an LP problem

There are several ways to express an LP problem. In this paper, except in a few cases, we will use the standard form to describe the LP problems, as it is shown in Problem 2.10. All the LP problems can be transformed into standard form. For the objective function, minimizing Z(λ λ λ ) is equivalent to maximizing -Z(λ λ λ ). The constraint a a a i λ λ λ ≤ b i can be transformed into equality by adding a slack variable s i ≥ 0 for each constraint C i : a a a i λ λ λ + s i = b i . We can also express the constraint a a a i λ λ λ ≤ b i as r i = a a a i λ λ λ , r i ≤ b i , where r i is called row variable. In other words, there are two alternative ways to transform the constraint C i : a a a i λ λ λ ≤ b i into standard form: C i : a a a i λ λ λ + s i = b i , s i ≥ 0 and r i = a a a i λ λ λ , r i ≤ b i , which are equivalent. If a decision variable λ i does not have the lower bound 0, we can split it into two variables:

λ i = λ ′ i -λ ′′ i , λ ′ i , λ ′′ i ≥ 0.
Example 2.5. Assuming we have an LP problem in Problem 2.11, and we will transform it into standard form.

minimize

Z(λ λ λ ) = 2λ 1 + λ 2 subject to λ 1 -2λ 2 ≤ -4 2λ 1 + λ 2 ≥ 1 (2.11)
First we transform the objective function into:

maximize Z(λ λ λ ) = -2λ 1 -λ 2
Then we unify the constraints with ≤ operator and obtain {λ 1 -2λ 2 ≤ -4, -2λ 1 -λ 2 ≤ -1}. Then to make λ λ λ ≥ 0, we split the variables as

λ 1 = λ ′ 1 -λ ′′ 1 , λ 2 = λ ′ 2 -λ ′′ 2 , where λ ′ 1 , λ ′′ 1 , λ ′ 2 , λ ′′ 2 ≥ 0. The orig- inal constraints become {λ ′ 1 -λ ′′ 1 -2λ ′ 2 + 2λ ′′ 2 ≤ -4, -2λ ′ 1 + 2λ ′′ 1 -λ ′ 2 + λ ′′ 2 ≤ -1, λ ′ 1 , λ ′′ 1 , λ ′ 2 , λ ′′ 2 ≥ 0}.
By adding the slack variables, we have the LP in standard form, which is shown in Problem 2.12.

maximize Z(λ λ λ ) = -2λ 1 -λ 2 subject to λ ′ 1 -λ ′′ 1 -2λ ′ 2 + 2λ ′′ 2 + s 1 = -4 -2λ ′ 1 + 2λ ′′ 1 -λ ′ 2 + λ ′′ 2 + s 2 = -1 and λ ′ 1 , λ ′′ 1 , λ ′ 2 , λ ′′ 2 , s 1 , s 2 ≥ 0 (2.12)
Alternatively Problem 2.13 shows the standard form using row variables.

maximize The LP problem is optimal at λ λ λ * , and the optimal value is Z * such that:

Z(λ λ λ ) = -2λ 1 -λ 2 subject to r 1 = λ ′ 1 -λ ′′ 1 -2λ ′ 2 + 2λ ′′ 2 + s 1 r 1 ≤ -4 r 2 = -2λ ′ 1 + 2λ ′′ 1 -λ ′ 2 + λ ′′ 2 r 2 ≤ -1 and λ ′ 1 , λ ′′ 1 , λ ′ 2 , λ ′′ 2 ≥ 0 (2.
Z * = Z(λ λ λ * ) (2.15)
A solution is said to be feasible if it fulfills all the constraints; and a feasible solution is optimized when the variables appear in the objective function cannot increase/decrease anymore.

Example 2.6. Assume we have an LP problem shown in Problem 2. [START_REF] Fouilhé | Revisiting the abstract domain of polyhedra: constraints-only representation and formal proof[END_REF].

maximize Z(λ λ λ ) = λ 1 -3λ 2 + 2λ 3 -λ 4 subject to 3λ 1 -λ 2 + λ 3 = 5 λ 1 -3λ 2 + 3λ 4 = -3 and λ λ λ ≥ 0 (2.16)
The objective function will be optimized at λ λ λ * = (0, 1, 6, 0), and Z * = 9. The objective function is equivalent to Z ′ (λ λ λ ) = - 16 3 λ 1 -2λ 4 + 9 by substituting

λ 2 = 1 3 λ 1 + λ 4 + 1, λ 3 = -8 3 λ 1 + λ 4 + 6.
As λ λ λ ≥ 0, λ 1 , λ 4 cannot decrease anymore when they reach 0, and thus the value of Z ′ (λ λ λ ) cannot increase anymore, i.e, it reaches optimum.

Simplex algorithm

We use the Simplex algorithm [START_REF] George B Dantzig | Maximization of a linear function of variables subject to linear inequalities[END_REF][START_REF] George | Linear programming 2: theory and extensions[END_REF][START_REF] Chvatal | Linear programming[END_REF] to solve the LP problems, so we explain it in details in this section. We first introduce some useful denominations, and then explain each step of the Simplex algorithm with an example. There are primal Simplex algorithm and dual Simplex algorithm, but we will only present the former. The later is a process of applying the primal Simplex algorithm to the dual of the original problem. In the following context, when we mention Simplex, we mean the primal Simplex algorithm.

The standard Simplex algorithm contains two phases: the first phase looks for a feasible solution, and the second searches the optimum. The constraints of an LP problem define a convex polyhedron. The Simplex algorithm searches the solution λ λ λ * among the vertices of this polyhedron. The first phase either starts from a vertex if the original problem is already feasible, or starts from exterior of the polyhedron otherwise. Once it reaches a vertex, the second phase starts.

Dictionaries An LP problem in the Simplex algorithm can be described by a dictionary, which is a system of linear equations. It substitutes variables by their equivalent counterpart expressed in terms of other variables. Consider Problem 2.10 which contains m constraints with n variables. If n > m, then m variables can be expressed by the left nm variables, as it is shown in Problem 2.17.

Z(λ λ λ ) = n ∑ j=m+1 o ′ j λ j + c λ i = b i - n ∑ j=m+1 a ′ i j λ j (i = i, • • • , m) (2.17)
where c is a constant. If we set the right-hand side variables to zero, the value of left-hand side variable λ i will be b i . As all the variables are non-negative, if all the b i has non-negative value, we will obtain a feasible solution. The corresponding dictionary is called feasible dictionary.

Basic and non-basic variables

The variables appear on the left-hand side of a dictionary are called basic variables, and that on the right-hand side are nonbasic variables. The set of basic variables is called a basis of the dictionary. The basic and non-basic variables form a partition of the variables, and the new objective function is obtained by substituting the basic variables with non-basic variables.

Example 2.7. Consider the Example 2.6. When the objective function reaches optimum, the basic variables are λ 3 , λ 4 , and nonbasic variables are λ 1 , λ 2 . We have

λ 3 = -3λ 1 + λ 2 + 5, λ 4 = -1 3 λ 1 + λ 2 -1. The new obtained objective function is Z ′ (λ λ λ ) = -14 3 λ 1 -2λ 2 + 11.
Entering and leaving variables The process of constructing a new linear system by switching a basic variable with a nonbasic variable is called pivoting. The chosen basic variable is called leaving variable (as it will leave the basis), and the chosen nonbasic variable is called entering variable.

First phase of the Simplex algorithm

If the initial LP problem is feasible, the first phase will be skipped; otherwise an auxiliary problem will be solved for obtaining a feasible solution. Assuming Problem 2.10 is not feasible, its auxiliary problem is shown in Problem 2.18.

maximizew

subject to n ∑ j=1 a i j λ j -w = b i (∀i ∈ {1, • • • , m}) and λ j ≥ 0, w ≥ 0 (∀ j ∈ {1, • • • , n}) (2.18)
The auxiliary problem can always be transformed into feasibility by choosing w as the entering variable. The original LP problem is feasible if and only if the auxiliary problem has optimal solution and the optimum is -w = 0. In other words, if the auxiliary problem reaches optimum but w is in the basis, the original problem is still infeasible. The feasible solution will be obtained by eliminating w, which is the nonbasic variable in the final system. Let us illustrate the process with Example 2.8.

Example 2.8. maximize Z(λ λ λ ) = 3λ 1 + λ 2 -2λ 3 -3λ 4 subject to λ 1 + 3λ 2 -2λ 3 = 4 2λ 1 + λ 3 -3λ 4 = -2 and λ j ≥ 0 (∀ j ∈ {1, • • • , 4})
First we substitute λ 1 , λ 2 by λ 3 , λ 4 and obtain:

Z(λ λ λ ) = - 8 3 λ 3 + λ 4 - 4 3
λ 1 = -1 - 1 2 λ 3 + 3 2 λ 4 λ 2 = 5 3 + 5 6 λ 3 - 1 2 λ 4
By setting λ 3 , λ 4 to zero, we can obtain λ 1 = -1, which means that the dictionary is infeasible. To obtain a feasible solution, we need to perform the first phase of the Simplex algorithm by solving the following auxiliary problem. maximize Z(λ λ λ , w) = -w

subject to λ 1 + 1 2 λ 3 - 3 2 λ 4 -w = -1 λ 2 - 5 6 λ 3 + 1 2 λ 4 -w = 5 3 and λ j ≥ 0, w ≥ 0 (∀ j ∈ {1, • • • , 4})
By substituting λ 1 and λ 2 we can obtain the corresponding dictionary:

Z(λ λ λ , w) = -w λ 1 = -1 - 1 2 λ 3 + 3 2 λ 4 + w λ 2 = 5 3 + 5 6 λ 3 - 1 2 λ 4 + w (2.19)
To obtain a feasible dictionary, in the first pivoting we need to choose w as the entering variable and λ 1 as the leaving variable. Then we will obtain the dictionary as follows:

Z(λ λ λ , w) = -1 -λ 1 - 1 2 λ 3 + 3 2 λ 4 w = 1 + λ 1 + 1 2 λ 3 - 3 2 λ 4 λ 2 = 8 3 + λ 1 + 4 3 λ 3 -2λ 4
(2.20)

Then we choose λ 4 to enter the basis, and w to leave. Finally we obtain a feasible dictionary:

Z(λ λ λ , w) = -w λ 4 = 2 3 + 2 3 λ 1 + 1 3 λ 3 - 2 3 w λ 2 = 4 3 - 1 3 λ 1 + 2 3 λ 3 + 4 3 w (2.21)
As w ≥ 0, the optimal function reaches the optimum 0. By eliminating w, we obtain the feasible dictionary of the original problem, in which λ 2 and λ 4 are basic variables.

Z(λ λ λ ) = 2 3 λ 1 - 7 3 λ 3 - 2 3 λ 4 = 2 3 + 2 3 λ 1 + 1 3 λ 3 λ 2 = 4 3 - 1 3 λ 1 + 2 3 λ 3 (2.22)

Second phase of the Simplex algorithm

The second phase of the Simplex will find the optimal solution or report that the problem does not have optimum, meaning that the objective function is unbounded.

In the standard form, the objective function needs to be maximized, and all the variables must satisfy λ i ≥ 0. Thus when all the coefficients of the nonbasic variables are non-positive, the objective function reaches optimum; otherwise the objective function can be increased. Then one nonbasic variable whose coefficient is positive should enter the basis. Given the entering variable λ j , the leaving variable is

λ k = b ′ i -∑ j a ′ i j λ j , which minimizes the ratio b ′ i a ′ i j
with a ′ i j > 0, where j ̸ = k and λ j are nonbasic variables.

Example 2.9. Let us consider Problem 2.22 of Example 2.8. We obtained a feasible dictionary, but the variable λ 1 in the objective function has a positive coefficient 2 3 , meaning that the objective function does not reach the optimum. It can be increased by augmenting the value of λ 1 . Thus we choose the variable whose coefficient is positive, i.e. λ 1 , as the entering variable in the next pivoting. Recall that the constraints are in the form λ k = b ′ i -∑ j a ′ i j λ j . The only choice of leaving variable is λ 2 , as in the first constraint λ 4 = 2 3 -(-2 3 λ 1 -1 3 λ 3 ) the coefficient of λ 1 is negative. By doing the pivoting, we obtain the final dictionary shown in Problem 2.23, in which the objective function cannot be furthermore increased since the maximal value is only obtained by setting the nonbasic variables λ 2 and λ 3 to 0, and any other choice would decrease the value of Z(λ λ λ ).

Z(λ λ λ ) = -2λ 2 -λ 3 + 2 λ 4 = 10 3 -2λ 2 + 5 3 λ 3 λ 1 = 4 -3λ 2 + 2λ 3 (2.23)
All the nonbasic variables are set to their lower bound 0, and we obtain the optimal vertex (4, 0, 0, 10 3 ) with the optimal value 2.

Simplex tableau

In the Simplex algorithm, the dictionary is normally represented in Simplex tableau. In the follow sections, we will always use the Simplex tableau to illustrate a dictionary. The Simplex tableau of Problem 2.10 is shown in Table 2.3, where o i are the coefficients of the objective function.

λ 1 • • • λ n = constant . . . • • • • • • • • • • • • row i a i1 • • • a in b i . . . • • • • • • • • • • • • objective o 1 • • • o n -Z * Table 2.3 -Simplex tableau.
The last dictionary of Example 2.8 which gives the optimum can be written in the Simplex tableau as shown in Table 2.4.

λ 1 λ 2 λ 3 λ 4 = constant row 1 0 2 -5 3 1 10 3 row 2 1 3 -2 0 4 objective 0 -2 -1 0 -2
Table 2.4 -Simplex tableau example.

In the Simplex tableau, the objective function is written in the form Z(λ λ λ ) -Z * = 0. In our example, the objective function row in Table 2.4 is -2λ 2λ 3 -Z * = -2, so we have Z * = 2 by setting λ 2 , λ 3 to 0.

Parametric linear programming

In this section we will catch a glimpse of the parametric linear programming (PLP), using which we compute the polyhedral projection and convex hull. We will see the definition of the PLP, and also the processes of solving a PLP problem.

Definition of PLP

A PLP problem is a linear optimization problem on the variables

λ λ λ = [λ 1 , • • • , λ m ] T whose objective function contains parameters x x x = [x 1 , • • • , x n ] T . Problem 2.24 shows the standard form of a PLP problem 1 . maximize Z(λ λ λ , x x x) = m ∑ i=1 ( n ∑ j=1 o i j x j -c i )λ i subject to m ∑ j=1 a i j λ j = b i ∀i ∈ {1, • • • k} and λ λ λ ≥ 0 (2.24)
The problem can be reformulated in a more concise way by introducing an extra dimension to the vector of parameters for representing the constant part

c i . Given O ∈ Q (n+1)×m , A ∈ Q k×m , b b b ∈ Q k and x x x = [x 1 , • • • , x n , 1] T , Problem 2.

Methods to solve PLP problems

The PLP problems can be solved directly by Simplex algorithm, or by solving a set of LP problems which are obtained by instantiating the parameters x x x. We will use both methods in our algorithm of PLP solver. The optimal solution of a PLP problem is a set of couples (R i , Z * i (x x x)), where R i is the region of the space of parameters x x x, in which the basis does not change and Z * i (x x x) is the optimal function corresponding to R i , meaning that any instantiation point x x x in R i leads to the optimal function Z * i (x x x). The region is obtained from the sign condition from the objective function, i.e., all the coefficients should be non-positive (or non-negative) if the objective is to be maximized (or minimized). An example of the Simplex tableau of PLP is shown in Table 2.5. The coefficients of the nonbasic variables in the objective function are linear expressions, and that of the basic variables are 0. Assuming we are computing the maximization of the objective function, the coefficients should be non-positive. By this condition we obtain the region i∈{1,

••• ,q} ∑ n k=1 o ik x k + c i ≤ 0.
Solving PLP via Simplex algorithm First we apply the first phase of the Simplex algorithm (Section 2.4) for obtaining a feasible solution with the basic and nonbasic variables. Then we express the original objective function of the PLP problem in terms of the nonbasic variables for obtaining an optimal function with the corresponding region. The next step is to pivot the Simplex tableau for obtaining other optimal functions. Considering the example in Table 2.5, we select one from the nonbasic variables, say λ q , to enter the basis. The leaving variable is λ r , which makes the ratio b j a jq be minimum for all j, where a jq > 0. If ∀ j, a jq ≤ 0, then λ q cannot enter the basis. We will take an example to explain this method in Example 2.10.

non-basic variables basic variables

= constants

λ 1 • • • λ q • • • • • • λ r • • • λ m . .
x k + c 1 • • • ∑ n k=1 o qk x k + c q • • • • • • 0 • • • 0 -Z * (x x x)
Table 2.5 -Simplex tableau.

Solving PLP via LP problems We randomly select an instantiation point in the parameter space, and substitute the parameters x x x by this instantiation point for obtaining an LP problem. If the LP is optimal at the vertex λ λ λ * , the PLP will also reach optimum at the same vertex λ λ λ * . By solving this LP problem we can obtain an optimal solution and the corresponding basic and nonbasic variables, if it is not unbounded. By expressing the objective function in terms of the nonbasic variables, we obtain an optimal solution (R i , Z * i (x x x)) of the PLP. Then we select another instantiation point outside R i , and repeat the same steps for obtaining other optimal solutions. This method will be shown on an example in Example 2.11.

=        -1 2 1 3 2 1 3 -2 0 -1 3 2       
, and

x x x = x 1 x 2 1 T . maximize Z(λ λ λ ) = (-x 1 + 2x 2 )λ 1 + (2x 1 + x 2 -1)λ 2 + (x 1 + 3x 2 + 3)λ 3 + (3x 1 -2x 2 + 2)λ 4 subject to 3λ 1 -λ 2 + λ 3 = 5 λ 1 -3λ 2 + 3λ 4 = -
-λ 2 + λ 3 -w = 5 λ 1 -3λ 2 + 3λ 4 -w = -3
and λ λ λ ≥ 0, w ≥ 0 (2.27)

In the next pivoting, we choose λ 3 as the entering variable. We can see that the coefficients of λ 3 in row 1 and row 2 are both positive, so we need to compare their ratios of the constant to the coefficient. As 7 4 / 1 8 > 9 4 / 3 8 , we choose λ 1 as the leaving variable. Then we obtain the new basis λ 2 , λ 3 , and the Simplex tableau is shown in Table 2.7. We obtain the optimal solution Z *

2 (x x x) = 8x 1 + 19x 2 + 17 with the region R 2 = {-3x 1 -17 3 x 2 -25 3 ≤ 0, 6x 1 + 2x 2 + 4 ≤ 0}.
As the coefficients of λ 4 in all rows are negative, λ 4 cannot enter the basis. Now all the possible bases have been found, and the algorithm terminates.

λ 1 λ 2 λ 3 λ 4 = constant row 1 0 1
λ 1 λ 2 λ 3 λ 4 = constant row 1 -1 3 1 0 -1 1 row 2 8 3 0 1 -1 6 objective -3x 1 -17 3 x 2 -25 3 0 0 6x 1 + 2x 2 + 4 -8x 1 -19x 2 -17

Preliminaries

The set of LP problems which are obtained by instantiating the parameters x x x with any point in R 1 will reach optimum with the same basis λ 3 , λ 4 , and thus the PLP will obtain the same optimal function. To obtain other optimal functions, we need to instantiate x x x with the points outside R 1 . With the instantiation point (-2, -2), we obtain the optimal solution (R 1 , Z * 1 (x x x)). Until now all the optimal solutions have been found, as any other point outside R 1 and R 2 will lead to unbounded LP problems.

Summary of the chapter

We presented some fundamental knowledge which is useful for the explanation of our efficient PLP solver in the following chapters. We saw the convex polyhedra and its representations. In our approach we use the constraint-only representation. We introduced some operators in the polyhedra domain. Our approach aimed at computing the polyhedral projection and convex hull. As the minimization operator is required in our PLP solver, we also provided an efficient minimization algorithm based on the work in [START_REF] Maréchal | Efficient elimination of redundancies in polyhedra by raytracing[END_REF].

We introduced the floating point arithmetic and rounding errors. In the paper we use Equation 2.6 and Equation 2.9 to compute the thresholds which are used for comparing floating point numbers.

We illustrated the processes of solving an LP problem using the Simplex algorithm. The computation processes can be shown in a dictionary or in the Simplex tableau form. In the following context we mainly use the latter. The methods to solve a PLP problem are also been presented. The solution of a PLP is a set of tuples (R i , Z * i (x x x)), where Z * i (x x x) is an optimal function and R i is the corresponding region. Please be note that for two different regions R i and R j , their optimal functions Z * i (x x x) and Z * j (x x x) could be the same. This point will be explained into details in Chapter 4.

Chapter 3

Sequential Parametric Linear Programming Solver

In this chapter we will present the overview of the sequential algorithm of the PLP solver. This work was presented in our prior paper [START_REF] Yu | An efficient parametric linear programming solver and application to polyhedral projection[END_REF].

We provide a flowchart that illustrates the processes of computing the projection (or convex hull) using our PLP solver. We will see that the minimization operator is required at several steps in the PLP solver, so we will present the minimization algorithm before explaining the algorithm of the PLP solver. We adapted the algorithm of raytracing minimization [START_REF] Maréchal | Efficient elimination of redundancies in polyhedra by raytracing[END_REF] to our floating point arithmetic.

Our PLP solver is based on the work of Maréchal et al. in [START_REF] Maréchal | Scalable minimizing-operators on polyhedra via parametric linear programming[END_REF]. Different from their work, we mainly used floating point arithmetic, and the rational result will be constructed. We aimed to obtain the rational result using floating point computations without impact on soundness and precision. Soundness means that we need to obtain an over-approximate polyhedron of the correct one; precision means that the resulting polyhedron should be the same as that obtained by rational computation using the same algorithm.

Due to the imprecision of floating point computations, the result may be incorrect and the algorithm may be not able to terminate. To avoid these problems, we use floating point thresholds to test the result and use rational numbers to check and recompute the result when it is necessary.

The whole algorithm of the PLP solver is divided into several steps, and each step will be explained into details. During the explanation of the algorithm, our focus will be on the cooperation of rational and floating point arithmetic, and some specifics of floating point computations.

Flowchart of the algorithm

The Figure 3.1 shows the flow chart of our algorithm. The rectangles are processes and diamonds are decisions. The processes/decisions colored in orange are computed in floating-point arithmetic, and those in blue use rational numbers. The dotted red frames show cases that rarely happen, which means that most computations in our approach use floating-point numbers.

Minimization

In the PLP algorithm, the minimization is required to minimize the input polyhedron and the obtained regions, so we will present the algorithm in this section. We will present two methods of minimization, and both of them will be used in our PLP solver.

Redundant Constraints

In our computation, there are several steps at which redundant constraints must be removed by minimization of the polyhedron. Definition 3.1 (Redundant constraints). A constraint is said to be redundant if it can be removed without changing the shape of the polyhedron. 3.2, the red constraint C 4 : 2x 1 + 4x 2 ≥ 7 is redundant.

Example 3.1. Considering the constraints {C 1 : -x 1 + 2x 2 ≥ 2, C 2 : 2x 1 -x 2 ≥ 1, C 3 : -x 1 -x 2 ≥ -8, C 4 : 2x 1 + 4x 2 ≥ 7} in Figure

Minimization using Farkas' lemma

The redundancy of a constraint can be checked by Farkas' lemma (Theorem • There exists a y y y ′ ∈ R m such that A T y y y ′ ≥ 0 and b b b T y y y ′ < 0.

=            a 11 • • • a m1 0 . . . a 1n • • • a mn 0 b 1 • • • b m 1            , d d d = [a i1 , • • • , a in , b i ] T (3.
Example 3.2. Let use consider the same example as Example 3.1. To determine if the constraint C 4 of Example 3.1 is redundant with respect to {C 1 , C 2 , C 3 }, we need to solve the satisfiability problem ∃y 1 , y 2 , y 3 , y 0 ≥ 0 such that a a a 4 x x xb 4 = y 1 (a a a 1 x x xb 1 ) + y 2 (a a a 2 x x xb 2 ) + y 3 (a a a 3 x x xb 3 ) + y 0 , i.e., Problem 3.4.

-

y 1 + 2y 2 -y 3 = 2 2y 1 -y 2 -y 3 = 4 -2y 1 -y 2 + 8y 3 + y 0 = -7 y 0 , y 1 , y 2 , y 3 ≥ 0 (3.4)
We obtain y y y = [y 1 , y 2 , y 3 , y 0 ] T = [ 10 3 ,8 3 , 0, 7 3 ] T , such that 2x 1 + 4x 2 -7 = 10 3 (-

x 1 + 2x 2 -2) +
On the contrary, C 1 is irredundant, as the second statement of Farkas' lemma is true. By solving Problem 3.5, we can find y y y ′ = [4, 0, 1], and A T y y y

′ = [7, 4, 1, 1] T > 0, b b b T y y y ′ = -6 < 0. 2y 1 -y 2 + 2y 3 = -1 -y 1 -y 2 + 4y 3 = 2 -y 1 + 8y 2 -7y 3 + y 0 = -2 y 0 , y 1 , y 2 , y 3 ≥ 0 (3.5)

Raytracing method for minimization

When there are much more irredundant constraints than redundant ones, using Farkas' lemma is not efficient. To find out the irredundant constraints more efficiently, Maréchal et al. presented a novel method using raytracing to minimize polyhedra [START_REF] Maréchal | Efficient elimination of redundancies in polyhedra by raytracing[END_REF]. We improved this approach by using floating point arithmetic instead of rational computations. The use of floating-point numbers here will not cause a soundness problem for static analysis, because in the worst case we eliminate constraints that should not be removed. That results in a larger polyhedron which includes the correct one, i.e., the abstract polyhedron is over approximated. At several steps we need to check the results of floating point computation for guaranteeing the precision and soundness. We first review the minimization algorithm, and then present adaptations we made for floating point computation.

Overview of the raytracing minimization

There are two phases in raytracing minimization. The first phase aims at finding as many irredundant constraints as possible. Rays are launched from a point inside the polyhedron toward each boundary. The constraint whose boundary is firstly hit by a ray is irredundant. The remaining constraints whose redundancy status is not solved by the first phase will be determined by the second phase: if we can find an irredundancy witness point, then the corresponding constraint is irredundant.

Definition 3.2 (Irredundancy Witness

). An irredundancy witness of a constraint C i is a point that violates C i but satisfies the other constraints. A ray r r r i was launched in the first phase from the interior point toward the boundary of C i : a a a i x x x ≤ b i in the direction a a a i . To determine the redundancy status of C i , we need to solve the problem ∃x x x, a a a i x x x > b i , j̸ =i a a a j x x x ≤ b j , where the constraints C j : a a a j x x x ≤ b j are the set of constraints whose boundaries hit by r r r i i i before the boundary of C i is reached. If the satisfiability problem is infeasible, then C i is redundant as the irredundancy witness point cannot be found; otherwise a point p p p will be obtained, but p p p is not necessarily an witness of the irredundancy of C i (see Figure 3.4). Another ray r r r ′ i will be launched from the interior point to the point p p p. Then C i will be first hit by r r r ′ i if it is irredundant, and p p p is its irredundancy witness; or another point p p p ′ will be obtained by solving a a a i x x x > b i , j̸ =i a a a j x x x ≤ b j , k̸ =i a a a k x x x ≤ b k , where the boundaries of C k : a a a k x x x ≤ b k are first hit by r r r ′ i . By adding constraints incrementally and solving a set of satisfiability problems, either the boundary of C i will be first hit by a ray if it is irredundant, or an infeasible solution will be obtained otherwise.

Adaptations for floating point computation

Now that we review the raytracing minimization algorithm of [START_REF] Maréchal | Efficient elimination of redundancies in polyhedra by raytracing[END_REF], let use detail the technicalities required for adapting the algorithm to floating point arithmetic.

Solving satisfiability and optimization problems

We use GLPK 1 (GNU Linear Programming Kit) to solve the optimization problems in our approach. For solving the satisfiability problems, we just set the objective function to 0. GLPK provides several routines for solving the LP problems, which are written in ANSI C, but we only used the Simplex algorithm.

Launching rays from a central point The raytracing algorithm uses an interior point from which we launch rays in directions perpendicular to the boundaries of the constraints. If we simply ask for any point inside a bounded polyhedron, we can construct an LP problem to compute a point that satisfies the constraints and does not locate on the boundaries, i.e., satisfies i a a a i x x x < b i with strict inequalities. Since GLPK does not support inequalities, one could be tempted to shift the non-strict constraints a little to simulate strict ones. Considering the strict constraint C i : a a a i x x x < b i , we can obtain a non-strict constraint C ′ i by shifting it:

C ′ i : a a a i x x x ≤ b i -∥a a a∥δ (3.6)
where δ is a positive constant.

However choosing an appropriate quantity for δ is difficult: if it is small we will obtain a point which is too close to the boundaries; if it is large, we may end up with an infeasible problem. Hence we prefer to compute a point located in a "central position" of the polyhedron, which is far away from the boundaries.

Considering a point p p p which is inside a polyhedron P, d i = b i -a a a i p p p ∥a a a i ∥ is the distance from p p p to the boundary of constraint C i of P, where ∥a a a i ∥ is the Euclidean norm of the vector a a a i . We aim at pushing the point p p p to a position which is as far as possible from the boundaries. Problem 3.7 shows the LP for finding the central point by adjusting the point p p p and the variable l. The variables d i only appear for the sake of readability. maximize l

subject to d i = b i -a a a i p p p ∥a a a i ∥ , l ≤ d i (∀C i ∈ P) (3.7)
When the LP is optimized, we get an optimal vertex [p p p l], from which we can obtain the expected central point p p p (Figure 3.5(a)). For unbounded polyhedra, there is no central point, as the polyhedron does not have a "central" position. In this case it is safe to apply the shifting strategy presented in Problem 3.6: no matter how large δ is, we can always obtain an interior point. The Simplex algorithm will obtain a satisfiability point at one vertex of the shifted polyhedron. As it is shown in Figure 3.5(b), one point of the two will be found as the interior point. One exception is the case where the unbounded polyhedron contains parallel constraints (Figure 3.5(c)). In this case if we shift the constraints with a large threshold, no point could be found. But the LP in Problem 3.7 is suitable for this special situation.

We do not need to test the shape of the polyhedron for choosing the method to find the interior point. We always try to find a central point using Problem 3.7. If the optimization problem is unbounded, we then try to obtain an interior point by shifting constraints. Proof. If the constraints of the polyhedron are unsatisfiable, we cannot find a point p p p inside the polyhedron. This case can be easily detected. We here consider the problem of detecting the "flat" polyhedron whose interior is empty. It means that the inequality constraints entail an equality. Consider a satisfiable polyhedron i a a a i x x x ≤ b. There exists at least one point p p p that satisfies this polyhedron, and we must have b ia a a i p p p ≥ 0, i.e., d i ≥ 0. As l is maximized, there must exist d i such that l = d i , otherwise l can still be increased. Then we know that l ≥ 0. (If p p p is outside the polyhedron, we must have l < 0, and l is not maximized.)

If the polyhedron has empty interior, the point p p p must locate on the boundary of the polyhedron, then one of the d i is equal to 0, so l must be 0. Now we need to prove that if l = 0, P must have empty interior. It is equivalent to: if P has non-empty interior, then l > 0. Assume p p p is on the boundary of C i , i.e. d i = l = 0. As the polyhedron has non-empty interior, if we move the point p p p a little towards the interior of the polyhedron, we will have l = d i > 0. That means when l = 0 the objective function l is not maximized. Therefore, if the polyhedron has non-empty interior, l must be lager than 0.

We just proved that when l = 0 the polyhedron has empty interior, but the Simplex solver in GLPK uses floating-point arithmetic, and thus we cannot perform the test l = 0 if l is a floating point number. Therefore we use a relative large threshold 10 -2 for comparison. If l < 10 -2 , we will check the satisfiability of the obtained interior point by a a a i p p p ≤ b + t(a a a i , p p p), where t(a a a i , p p p) = 1.01nε|a a a i | T |p p p| is the threshold (Equation 2.9). If the satisfiability test fails, the algorithm will report that the floating point solver cannot find a central point. Then we switch to a rational solver to deal with the problem. We do not need to test the correctness of the optimization, as we just try to find a point at the center, and any point near the center is fine.

Computing an irredundancy witness point In the process of testing if a constraint C i : a a a i x x x ≤ b i is redundant with respect to other constraints C j , where j ̸ = i, we need to solve a satisfiability problem: ∃x x x, a a a i x x x > b i , j̸ =i a a a j x x x ≤ b j for computing the irredundancy witness point of the constraint C i (Section 3.2.3.1). For efficiency, we solve this problem in floating point using GLPK. However, GLPK does not support strict inequalities, thus we need tricks to deal with them.

One way is to shift the strict inequality constraint C i a little for obtaining a non-strict inequality C ′ i : a a a i i i x x x ≥ b i + ε, where ε is a positive constant. This method is however requires to provide a suitable ε, which is the same problem that we faced for computing the interior point. One exception is that when the polyhedron is a cone, it is always possible to find a point that satisfies a a a i x x x ≥ b i + ε, j̸ =i a a a i x x x ≤ b by shifting the constraint C i , no matter how large ε is.

We designed another method for non-conic polyhedra. Instead of solving a satisfiability problem, we solve the optimization problem shown in Problem 3.8. The found optimal vertex x x x * is the solution we are looking for if it satisfies a a a i x x x * > b i + t(a a a i , x x x * ) (see Equation 2.8 about t). maximize a a a i x x x subject to a a a j x x x ≤ b j (∀ j ̸ = i)

P = {C 1 : x 2 ≤ 2, C 2 : x 1 + x 2 ≤ 7, C 3 : -x 1 + x 2 ≤ 0},
whose interior point is p p p. We want to test the redundancy status of the constraint C 1 . In Figure 3.6(a), the ray launched from p p p to C 1 hit C 2 first, and found the point p p p ′ which violates C 1 and satisfies C 2 . Then in Figure 3.6(b) the ray launched from p p p to p p p ′ hit C 3 first, so we need to find a point that violates C 1 and satisfies C 2 , C 3 , i.e., solving the satisfiability problem

{x 2 > 2, x 1 + x 2 ≤ 7, -x 1 + x 2 ≤ 0}.
In

Preprocess of raytracing minimization

Before giving the whole floating point raytracing minimization algorithm, there is another point to be explained. The polyhedron to be minimized may contain implicit equalities, i.e., an equality resulting from the combination of several inequalities, for instance

x 1 ≤ x 2 , x 2 ≤ x 3 , x 3 ≤ x 1 producing x 1 = x 2 = x 3 .
In this case the polyhedron to be minimized has an empty interior. As we have explained, the raytracing minimization seeks for a point inside the polyhedron to be minimized. When there are implicit equalities the raytracing minimization method cannot find the interior point, and thus we need to eliminate these implicit equalities before applying the raytracing minimization.

Eliminating implicit equalities Consider the polyhedron P = i a a a i x x x ≤ b i . The way to remove the implicit constraints is that we select the constraints that may construct the implicit constraints in floating point computation, and then from these selected constraints we use a rational solver to search for the combination of inequalities constructing the equality. The purpose of using the floating point solver is to decrease the number of constraints that will be given to the rational solver and to improve the efficiency.

As what we have mentioned in Section 3.2.3.2, the minimization algorithm contains the step of testing the interior emptiness of the polyhedron P. If the floating point solver reports that P has an empty interior, there may exist implicit inequalities. Then we try to find a point p p p which fulfills all the constraints C i : a a a i x x x ≤ b i of P. If p p p can be found, P probably contains implicit equalities. In this case we select all the constraints C i fulfill a a a i p p p ≤ b i + t(a a a i , p p p), i ∈ {1, • • • , k}, i.e, the constraints closed to p p p (see Equation 2.8 about t). Among these selected constraints we look for a combination producing the equation a a a i x x xb i = -∑ α j (a a a j x x xb j ), α j ≥ 0, i ̸ = j. If this combination cannot be found, it means that there is no implicit equalities; otherwise we will obtain a set of equalities. By substituting the equalities into the inequalities, we will obtain the reduced inequalities, which is explained in Example 3.5. All the implicit equalities will be eliminated by repeating this process.

Example 3.5. Given a polyhedron P = {C 1 :

2x 1 + x 2 ≤ 5, C 2 : x 1 -3x 2 ≤ 2, C 3 : -x 1 + 3x 2 ≤ -2},
the constraints C 2 and C 3 construct an equality x 1 -3x 2 = 2. As the equality does not appear in the polyhedron explicitly, it is an implicit equality. Once obtain the equality x 1 -3x 2 = 2, we substitute x 1 by x 2 , i.e. x 1 = 3x 2 + 2. By substituting x 1 of the remaining inequality 2x 1 + x 2 ≤ 5, we obtain 7x 2 ≤ 1.

The floating point raytracing minimization algorithm

The pseudo-code of the algorithm of the raytracing minimization is shown in Algorithm 1. As we said, Q denotes the field of rational numbers, and F is the set of finite floating-point numbers. The variables stored as floating point numbers are denoted by name F .

Termination There is another point to be explained about Algorithm 1. During the process of computing the witness point, Algorithm 1 may loop forever because of imprecision of floating point In exact computation the algorithm should report that C 1 is irredundant. But if we compute in floating point numbers, when the gray area is extremely small, the distance l 1 from p p p to C 1 and the distance l 2 from p p p to C 2 are nearly equal. Because of rounding errors, the floating point solver may misjudge that l 2 < l 1 , i.e., the ray r r r hit C 2 first. Then it would find the same point p p p ′ , which violates C 1 and satisfies C 2 . Then the same ray r r r would be launched from p p p to p p p ′ . Because of the misjudge of the comparison of l 1 and l 2 , the process would be repeated, and Algorithm 1 would loop for ever.

p F = GetCentralPoint(poly F ) if LP is unbounded then // if
To avoid this problem, we need to maintain a list of constraints which are hit before reaching C 1 . If the hit constraint exists in the list, the process of checking redundancy should stop and report the constraint C 1 as redundant. Soundness The algorithm of raytracing minimization using floating point arithmetic is sound, meaning that the polyhedron which is minimized by the floating point minimization algorithm always contains the exact minimized polyhedron. In other words, the minimization algorithm should always report a constraint as redundant when it cannot determine the redundancy status of this constraint.

In the first phase of raytracing minimization, the determination of irredundancy consists in comparing the distances from the interior point to each constraint, which are floating point numbers. A constraint C i is irredundant only if d i + t < d j , ∀ j ̸ = i, where d i (or d j ) is the distance from the interior point to the constraint C i (or C j ), t = max(|d i |, |d j |)ε is the threshold and ε is the machine epsilon. If there are several constraints C j satisfy |d id j | < t, then C i will be marked as undetermined and will be remained to the second phase. In the second phase, if GLPK cannot find an irredundancy witness point w w w that fulfills C i : a a a i x x x ≤ b i , we can believe it and mark C i as redundant; if GLPK found the irredundancy witness point, we need to test that the point really violates C i with a threshold t: a a a i w w w > b i + t(a a a i , p p p) (see Equation 2.8 about t). If the test fails, the constraint C i will be marked as redundant.

We have seen that in the second phase, a constraint is reported as redundant if: i) GLPK cannot find the witness point; ii) the thresholded satisfiability test a a a i w w w > b i + t(a a a i , p p p) fails; iii) there is numerical errors during the computation of witness point (the same hit constraint appears repeatedly).

As a result, some irredundant constraints may be misjudged as redundant, which results in overapproximation, but not unsoundness in static analysis.

Projection and Convex hull

We present in this section how the PLP solver can be used to compute projection and convex hull, following the work of Maréchal et al. [START_REF] Maréchal | Scalable minimizing-operators on polyhedra via parametric linear programming[END_REF]. Then we will be able to explain our efficiency improvement using floating point computations.

Projection

Given a polyhedron P in d dimensions, the polyhedral projection can eliminate 1 to d -1 dimensions and obtain a polyhedron P ′ in lower dimensions such that

P ⊏ P ′ × Q 1 × • • • × Q e for elimination of dimensions 1 to e.
Example 3.6. Figure 3.9 shows an example of projecting a 3D polyhedron to 2D. The polyhedron to be projected (in gray) is

P = {-15x 1 -x 2 + 6x 3 ≥ 0, -3x 1 -4x 2 + 5x 3 ≥ 0, -2x 2 + x 3 ≥ -3, -x 2 + 2x 3 ≥ 0, 6x 1 +x 2 +10x 3 ≥ 0, x 1 +4x 2 -24x 3 ≥ -59, 2x 1 +8x 2 +11x 3 ≥ 0}.
We eliminate the variable x 3 and obtain the projected polyhedron (in purple) Fourier-Motzkin elimination [START_REF] Khachiyan | Fourier-motzkin elimination method[END_REF] is a classic method to project polyhedra in constraint-only description. This method has two drawbacks: i) it can only eliminate one variable each time, ii) many redundant constraints could be produced during the process. Those drawbacks are the motivation to compute projection using PLP [START_REF] Jacob | Polyhedral analysis using parametric objectives[END_REF][START_REF] Jones | On polyhedral projection and parametric programming[END_REF][START_REF] Maréchal | Scalable minimizing-operators on polyhedra via parametric linear programming[END_REF].

P ′ = {-x 1 ≥ -1, 77x 1 + 32x 2 ≥ -259, x 1 - 44x 2 ≥ -131, x 1 + 4x 2 ≥ -11}.
PLP for projection The polyhedron to be projected is P: Ax x x + b b b ≥ 0. Assume we will eliminate x p , • • • , x q , where 1 ≤ p ≤ q ≤ n. To compute the projection, following [START_REF] Maréchal | Scalable minimizing-operators on polyhedra via parametric linear programming[END_REF] we construct a PLP shown in Problem 3.9, in which x x x are the parameters, and λ λ λ are decision variables, with

x x x = [x 1 , • • • , x n ] T , λ λ λ = [λ 0 , • • • , λ m ] T . minimize m ∑ i=1 λ i (a a a i x x x + b i ) + λ 0 subject to m ∑ i=1 λ i (a a a i p p p + b i ) + λ 0 = 1 (*) n ∑ i=1 a i j λ i = 0 (∀ j ∈ {p, • • • , q})
and

λ i ≥ 0 (∀i ∈ {0, • • • , m}) (3.9)
where

p p p = [p 1 , • • • , p n ]
T is a point inside the projected polyhedron and it is called normalization point.

The (*) equation in Problem 3.9 is called normalization constraint, which makes sure that: i) at least one of the λ i is not zero; ii) all the constraints of regions intersect at the same point p p p; iii) the obtained projected polyhedron is free of redundancy (see [START_REF] Maréchal | New Algorithmics for Polyhedral Calculus via Parametric Linear Programming[END_REF] for explanation). Let us consider the projected polyhedron shown in Figure 3.9. Figure 3.10(a) shows the regions obtained from the normalized PLP problem, i.e., Problem 3.9. If we replace the normalization constraint with ∑ m i=0 λ i = 1, which only guarantees that not all the variables are zero, we obtain the regions without normalization in Figure 3.10(b). In these two figures, the regions with the same color have the same optimal function. We can see that in

Convex hull

The convex hull of two polyhedra is the smallest convex set that contains the two polyhedra. Given two polyhedra P 1 : A ′ x x x ′ + b b b ′ ≥ 0, P 2 : A ′′ x x x ′′ + b b b ′′ ≥ 0, the convex hull of P 1 and P 2 is:

P 1 ⊔ P 2 = {x x x|x x x = αx x x ′ + (1 -α)x x x ′′ , A ′ x x x ′ + b b b ′ ≥ 0, A ′′ x x x ′′ + b b b ′′ ≥ 0, 0 ≤ α ≤ 1}.
The Figure 3.11 shows an example of the convex hull of two polyhedra in 2D. The orange areas are the polyhedra, and the red frame is the boundary of the convex hull. PLP for convex hull We recall here the method designed by Maréchal et al. [START_REF] Maréchal | Scalable minimizing-operators on polyhedra via parametric linear programming[END_REF] to compute the convex hull via PLP. Let P = P 1 ⊔ P 2 with P 1 : A ′ x x x ′ + b b b ′ ≥ 0 and P 2 :

A ′′ x x x ′′ + b b b ′′ ≥ 0, where A ′ ∈ Q m ′ ×n , b b b ′ ∈ Q m ′ , A ′′ ∈ Q m ′′ ×n and b b b ′′ ∈ Q m ′′ .
As P 1 ⊑ P and P 2 ⊑ P, each constraint of P is redundant with respect to both P 1 and P 2 . According to what we have explained in Section 3.2.2, for each constraint

C i : a a a i x x x + b i ≥ 0 of P, we have ∃λ λ λ ′ ≥ 0, a a a i x x x + b i = ∑ m ′ j=1 λ ′ j (a a a ′ j x x x ′ + b ′ j ) + λ ′ 0 = ∑ m ′′ j=1 λ ′′ j (a a a ′′ j x x x ′′ + b ′′ j ) + λ ′′ 0 .
Thus we have the constraints (**) in Problem 3.10. The constraint (*) is the normalization constraint, as in Problem 3.9.

minimize m ′ ∑ i=1 λ ′ i (a a a ′ i x x x ′ + b ′ i ) + λ ′ 0 subject to m ′ ∑ i=1 λ ′ i (a a a ′ i p p p + b ′ i ) + λ ′ 0 = 1 (*) A ′T λ λ λ ′ -A ′′T λ λ λ ′′ = 0 b ′T λ λ λ ′ + λ ′ 0 -b ′′T λ λ λ ′′ -λ ′′ 0 = 0 (**) and λ ′ i ≥ 0, λ ′′ j ≥ 0 (∀i ∈ {0, • • • , m ′ }, ∀ j ∈ {0, • • • , m ′′ }) (**) (3.10)
where p p p = [p 1 , • • • , p n ] T is the normalization point.

The sequential algorithm of PLP solver for computing projection and convex hull

In this section we present our sequential algorithm of PLP solver. We will see how to solve the PLP problem in floating point numbers and reconstruct the rational results. Some checkers are provided to ensure the soundness and precision of the results. When it is necessary, rational solvers can be invoked for computing exact solutions.

Overview of the sequential algorithm

In the pseudo-code of the sequential algorithm of the PLP solver (Algorithm 2), we annotate data with name type , where name is the name of data and type is either Q or/and F. Q × F means that the data is stored in both rational and floating-point numbers. The PLP problem is stored in both floating point and rational numbers. The optimal functions are in rational numbers, and the corresponding regions are in both floating point and rational numbers.

We use a work list to store a set of tasks. A task is stored as a triple (w w w, R f rom , F f rom ), in which w w w is generated by the region R f rom and it is the irredundancy witness point of the constraint whose boundary is F f rom . The task points w w w are in the parameters' domain, and we use them to compute unexplored regions. An initial task will be added into the work list, and more tasks will be appended during the execution of the algorithm.

Given a PLP problem, the solution goes through several steps. We first list the steps, and then explain each step into details in the following sections.

Step 1 Construct the PLP for computing projection or convex hull.

Step 2 Pick up a random valuation of the parameters for obtaining a initial task (w w w, none, none) and push it into the work list W .

Step 3 Fetch a task (w w w, R f rom , F f rom ) from W .

Step 4 Check if the point w w w is covered by any known region. If it is not, go to Step 5; otherwise go to Step 10.

Step 5 Instantiate the parametric objective with w w w to obtain a standard (non-parametric) LP problem.

Step 6 Solve the LP problem obtained in Step 5 using GLPK to obtain a basis B.

Step 7 Reconstruct the PLP problem which is equivalent to the original one associated to B.

Step 8 Minimize the region, and store the current solution associated to B.

Step 9 Add more tasks of which the task points are outside the obtained region into the work list.

Step 10 Check if the obtained region and the region R f rom are adjacent.

Step 11 Go to Step 2 and repeat the processes until W = / 0.

Construction of PLP problems (Step 1)

The first step of our approach is to construct the PLP problem for computing projection (Problem 3.9) or convex hull (Problem 3.10). For constructing the PLP problem, we need to choose a normalization point. By instantiating the parameters with this point, we obtain the normalization constraint. The set of constraints of the PLP problem are all equalities, and thus we need to eliminate the redundant equalities, i.e. the equalities which are equivalent to others, which can be performed by Gaussian Elimination. Then we substitute the equalities into the inequalities λ λ λ ≥ 0 for reconciling our approach with GLPK. After that we minimize the set of inequality constraints. The details are explained below. At the end of Step 1, we store the constraints of the PLP problem in both rational and floating point numbers. The rational constraints are used for reconstructing rational results, and the floating point constraints are given to GLPK for solving the LP problems. The constraints in floating point numbers are inequalities, and then the Simplex algorithm in GLPK will add slacks variables to each inequalities, i.e., the inequalities will be transformed into equalities in GLPK. Thus we store the rational constraints as equalities such that they are in the same form as in GLPK.

Choosing a normalization point Consider the projection operator implemented via PLP and P ′ = P/{x p , • • • , x q } is the resulting polyhedron, i.e., the polyhedron obtained by projecting out the variables {x p , • • • , x q } of P. The point p p p ′ used during PLP computation for normalization of P ′ must be in the interior of P ′ . p p p ′ is computed from a point p p p in the interior of P by p p p ′ = p p p/{x p , • • • , x q }, i.e., the components {x p , • • • , x q } of p p p are discarded. For convex hull, P ′ = P 1 ⊔ P 2 , the point p p p is chosen either in P 1 or in P 2 . As P 1 ⊑ P ′ and P 2 ⊑ P ′ , p p p ′ = p p p must be in the interior of the obtained polyhedron P ′ .

The normalization point p p p ′ in the projected polyhedron is computed in floating-point numbers. As the constraints of the PLP problem will be stored in both floating-point and rational numbers, we need to convert p p p ′ into rational for constructing the rational constraints, and make sure that the rational point is still inside the resulting polyhedron. The point p p p ′ ∈ P ′ if p p p ∈ P, so we compute p p p, convert it into rational numbers, and then check if p p p is still inside P. This is straightforward: firstly we keep k digits in decimal part and convert the value into rational by computing its continued fraction. If the rational point does not satisfy P, we keep k + 1 digits and do that repeatedly until a satisfiable point is obtained.

Applying Gaussian Elimination A PLP problem may contain constraints that are equivalent to each other, i.e., redundant equalities, which should be removed. We apply Gaussian Elimination to the augmented matrix (Definition 3.3) of the constraints of the PLP problem for obtaining the augmented matrix in reduced row echelon form2 , and thus all the redundant equalities will be removed. Definition 3.4 (Reduced row echelon form of a matrix). A matrix is in reduced row echelon form if it fulfills the conditions below:

• All rows with at least one nonzero element are above rows of all zeros.

• Each nonzero row contains an element 1.

• Each column containing the element 1 has zeros everywhere else. Reconciling with GLPK To solve LP problems, GLPK always adds row variables to the constraints, no matter they are equalities or inequalities. Assuming we have an equality a a a i λ λ λ = b i , which is expressed in GLPK as an equality r i = a a a i λ λ λ with a bound r i = b i , where r i is a row variable. The systematic introduction of row variables can introduce spurious dimension in the results of our PLP solver, and this will lead to wrong regions, which will be explained in Example 3.8. The idea to reconcile GLPK and our procedure is to prevent GLPK from adding useless row variables to equality constraints. Hence we transform the original problem into a reduced form without equality constraints. Let us illustrate it on an example.

Example 3.8. Consider a polyhedron {-40x 2 ≤ -7, -x 1 ≤ 20, -x 1 + 5x 2 + 3x 3 ≤ 20}, from which we want to eliminate the variable x 1 using our PLP solver. We construct the PLP problem for projection in Problem 3.11 with the normalization point (0, 1.2, 0).

minimize λ 1 (40x 2 -7) + λ 2 (x 1 + 20) + λ 3 (x 1 -5x 2 -3x 3 + 20) + λ 0 subject to 41λ 1 + 20λ 2 + 14λ 3 + λ 0 = 1 λ 2 + λ 3 = 0 and λ 1 , λ 2 , λ 3 , λ 0 ≥ 0 (3.11)
The PLP solver instantiates the parameters (x 1 , x 2 ) of the objective function with a point, say (1, 1, 1), to obtain a standard (non-parametric) LP problem shown in Problem 3.12. GLPK transforms this LP problem into Problem 3.13. )

minimize 33λ 1 + 21λ 2 + 13λ 3 + λ 0 subject to 41λ 1 + 20λ 2 + 14λ 3 + λ 0 = 1 λ 2 + λ 3 = 0 and λ 1 , λ 2 , λ 3 , λ 0 ≥ 0 (3.12) minimize subject to and 33λ 1 + 21λ 2 + 13λ 3 + λ 0 r 1 = 41λ 1 + 20λ 2 + 14λ 3 + λ 0 , r 1 = 1 r 2 = λ 2 + λ 3 , r 2 = 0 λ 1 , λ 2 , λ 3 , λ 0 ≥ 0 (3.
λ 1 = 1 41 r 1 -20 41 λ 2 -14 41 λ 3 -1 41 λ 0 , r 1 = 1 r 2 = λ 2 + λ 3 , r 2 = 0 λ 1 , λ 2 , λ 3 , λ 0 ≥ 0 (3.14)
We know that λ 2 = -λ 3 = 0, as λ 2 , λ 3 ≥ 0, and r 1 = 1, which means the objective function in Problem 3.14 is equivalent to ( 40 41 x 2 -7 41 ) + λ 0 (- 40 41 x 2 + 48 41 ), where λ 2 and r 1 are replaced by their fixed bound. Since λ 0 has a lower bound, the objective function reaches minimum when the parametric coefficient of λ 0 is non-negative, i.e. in the region {- 40 41 x 2 + 48 41 ≥ 0}. This is the expected result. But the information λ 2 = -λ 3 is not straightforward for GLPK because of the introduce of λ 2 . When our region reconstruction procedure requests GLPK to get the bound of the non-basic variables r 1 , λ 2 and λ 3 , it obtains fixed bound r 1 = 1 but lower bound λ 2 ≥ 0, λ 3 ≥ 0. Based on this information λ 2 and λ 3 remain in the objective function, which leads to the region {x 1 -800 41 x 2 + 960 41 ≥ 0, x 1 -765 41 x 2 -3x 3 + 918 41 ≥ 0, - 40 41 x 2 + 48 41 ≥ 0} in Q 3 , whereas we expected a region in Q 2 . The error is obvious, as x 1 should have been eliminated.

To avoid this issue, we transform the equalities into inequalities before giving them to GLPK. We substitute λ 1 , λ 2 in the inequalities and the objective function by λ 1 = 6 41 λ 3 -1 41 λ 0 + 1 41 , λ 2 = -λ 3 , and obtain the reduced form of the constraints of Problem 3.11 is shown in Problem 3.15. The substitution can be computed by Gaussian Elimination. The augmented matrix of the constraints

{41λ 1 + 20λ 2 + 14λ 3 + λ 0 = 1, λ 2 + λ 3 = 0} is    41 20
λ 1 = 6 41 λ 3 -1 41 λ 0 + 1 41 , λ 2 = -λ 3 .
subject to 6 41 λ 3 -

1 41 λ 0 + 1 41 ≥ 0 -λ 3 ≥ 0 and λ 3 , λ 0 ≥ 0 (3.15)
The next step is to remove the redundant constraints. As what explained in Section 3.2.3.3, we need to eliminate the implicit equalities before using the raytracing minimization. We found the implicit equality λ 3 = 0, substituted it into the inequality constraints and the objective function, and obtained our final LP problem shown in Problem 3.16. Then we need to use the raytracing minimization to remove the redundant constraints. In our example, the two constraints are both irredundant. GLPK will add a row variable to the constraint and obtain r 1 = λ 0 , r 1 ≥ 0. To keep the equivalent representation of the constraint in GLPK and in our approach, we add a slack variable λ 1 to the constraint and obtain the equality 1 41 λ 0 + λ 1 = 1 41 with λ 0 , λ 1 ≥ 0. We store this equality as the constraint of the PLP problem. The final PLP problem is shown in Problem 3.17

Initial tasks (Step 2)

In Step 2 of the PLP solver, one initial point, which is different from the normalization point, will be randomly selected to trigger the algorithm. The initial task point will be inserted into the work list firstly, and then more task points will be inserted during the computation approach.

Checking belonging of points (Step 3-4)

Before exploring a region using a task point (Step 5), we need to check if the point is covered by any known region for avoiding recomputation. Assume we fetch a task (q q q, R, F) from the work list, and test if the task point q q q is inside the region i o o o i x x x ≤ c i . As there is no strict inequality in floating-point arithmetic, we cannot simply test o o o i q q q < c i , and hence we test o o o i q q q ≤ c it(o o o i , q q q), where t is defined by Equation 2.8. If the test succeeds, it guarantees that q q q locates in the interior of the region. However, if it reports the point is not covered, the point may still be inside the region. This may lead to recomputation of the corresponding LP problem at a low probability. This approach guarantees that no task will be missed. The termination is guaranteed, as we maintain a set of bases (stored in a hash table) which are obtained from GLPK. If a basis is found to be existed in the hash table, the computation will terminate.

Obtaining an optimal basis by solving LP using GLPK (Step 5-6)

If the point q q q taken from the work list is not covered by any known region, we will construct the objective function of the LP problem by instantiating the parameters with q q q. The LP problem is:

minimize n ∑ i=1 λ i (a a a i q q q + b i ) + λ 0 subject to n ∑ i=1 λ i (a a a i p p p + b i ) + λ 0 = 1 (*) n ∑ i=1 a i j λ i = 0 (∀ j ∈ {p, • • • , q})
and

λ i ≥ 0 (∀i ∈ {0, • • • , n}) (3.18)
Then we solve the LP problem with the Simplex method in GLPK. GLPK reconstructs the constraints into the form (λ andc i is a constant. A new objective function is obtained by substituting the basic variables with nonbasic variables. We extract the indices of basic and non-basic variables from GLPK for reconstructing the rational results.

B ) i = ∑ n j=1 a i j (λ N ) j + c i , where (λ B ) i is a basic variable, (λ N ) i is a non-basic variable,

Reconstructing rational matrix and extracting result (Step 7-9)

We solved Problem 3.18 in floating-point numbers using GLPK. Had the solving been done in exact arithmetic, one could retain the optimal point λ λ λ * , but here we cannot use it directly. Instead, we exploit the partition of the variables into basic and non-basic variables, and from this partition we can recompute exactly, in rational numbers, the optimal function, as well as a certificate that it is feasible.

Consider the Simplex tableau below in Table 3

x k + c 1 • • • ∑ n k=1 o mk x k + c m -Z * (x x x) = ∑ n k=1 z k x k + z 0 Table 3.1 -Simplex tableau.
the matrix of the PLP objective function in Table 3.1. The ith column corresponds to the coefficients of ith variable, where the variables

λ λ λ = [λ 1 , • • • , λ m ]
; the last column of the each matrix represents the constants.

M =        a 11 • • • a 1m b 1 . . . . . . . . . . . . a n1 • • • a nm b n        O =            o 11 • • • o m1 z 1 . . . . . . . . . . . . o 1n • • • o mn z n c 1 • • • c m z 0            (3.19)
To generalize the result of the LP, we need to reconstruct the matrices M and O associated to the basis B provided by GLPK. We extract the columns corresponding to the basic variables from matrices in Problem 3.19, and use M B and O B to denote the sub-matrices of M and O only containing these extracted columns. By linear algebra in rational arithmetic 3 we compute a matrix Θ, representing the substitution performed by the Simplex algorithm, which result in the coefficients of the objective function of the PLP problem being 0. We compute the matrix Θ by Θ = O B M -1

B , where M -1 B denotes the inverse of M B (actually, we do not inverse that matrix but instead call a solver for linear equation systems). Then Detecting optimum in PLP We detect the optimal function and the corresponding region from the matrix O ′ . The variables of our PLP problem for computing projection have lower bound 0, which means that the objective function of the PLP problem reaches the optimal when all the non-basic variables reach their lower bound and their coefficients should be non-negative. Each non-zero column in O ′ represents a function in x x x, which is the coefficient of a non-basic variable, say λ j . The conjunction of constraints j / ∈B (O ′

• j ) T x x x ≥ 0 defines a region of x x x where j belongs to the indices of non-basic variables. Note that the region itself is a polyhedron. Therefore, the conjunction of constraints may contain redundancy and we must apply the minimization procedure over it. The minimization algorithm produces irredundancy witness points of all the irredundant constraints, which are points located over each constraint. Thus they can be used for instantiating the PLP objective and exploring new part of the parameter space. We insert these witness points into the work list and use them to do the next explorations.

Checking adjacency regions (Step 10)

Once we found a new region, the witness points of this region are added into the work list as new task points. Then more regions will be explored by these tasks. But there is a subtlety: the region found by the witness point may not be adjacent to the original region, i.e. they do not share the same boundary (see Figure 3.12). Adjacent regions Consider the example in Figure 3.12. Assume we found a new region R 1 , and one of its constraints is C 1 . The exploration using the witness point w w w 1 of C 1 finds a region R 2 which is not adjacent to R 1 , i.e., there is a gap between the boundary of C 1 belonging to R 1 and that of C 2 belonging to R 2 . It is possible that the witness point w w w 2 of C 2 will not find the region between R 1 and R 2 (w w w 2 could be in R 1 or beyond). Hence there is a risk of missing a region. To prevent missing regions, we need to check adjacency of regions once we found a new region. Considering the same example in Figure 3.12: once found the region R 2 with the point w w w 1 , we check adjacency of R 1 and R 2 . The adjacency test is designed to provide a point w w w ′ between the two regions if adjacency does not hold. Thus w w w ′ is inserted into the work list for exploring the missing region between R 1 and R 2 . The method to check adjacency To check adjacency we compare each constraint of the two regions syntactically in rational numbers. We can use a syntactic test thanks to the normalization constraint (*) of Problem 3.9. If two regions have a common boundary, they are likely to be adjacent, but this is not a sufficient condition. We need to consider the situation in Figure 3.13, in which the constraints C 1 and C 2 have a common boundary, which is denoted by F, but their regions R 1 and R 2 are not adjacent.

R1

To detect this situation we need an additional test. If two regions are adjacent, they should have the same optimal value along the direction of the common boundary. We choose a point randomly on the boundary F, and evaluate the optimal functions of R Finding a point between non-adjacent regions When the previous test reports that regions R 1 and R 2 are not adjacent, we compute an additional point between them, as it is shown in Figure 3.14. We compute the length l 0 between the two points w w w 1 and w w w 2 . The length from w w w 1 to the nearest frontier of C 1 is l 1 and that from w w w 2 to the frontier of C 2 is l 2 . The segment from w w w 1 to w w w 2 intersects the boundaries of C 1 and C 2 at w ′ 1 and w ′ 2 respectively. The new instantiation point w w w 3 is the middle of the segment w ′ 1 w ′ 2 . The adjacency test is performed in rational computation. The new point w w w 3 is computed in floating point arithmetic, as w w w 1 and w w w 2 are in floating point numbers.

Generally speaking, we do not need to consider the imprecision of floating point arithmetic here, as the point w w w 3 can locate at any place inside the missed region. What matters is that w w w 3 should not be too closed to the boundary, otherwise the exploration of regions may loop forever. Considering the same example in Figure 3.14, if the added point w w w 3 in R 3 is extremely closed to C 1 , it may find the region R 1 (although the algorithm will terminate once found that the basis has already been found). Then R 1 is checked to be not adjacent to R 2 , and the same point would be added. Repeatedly R 1 would be found, and so on. To avoid this situation, we check l 0 > l 1 + l 2 + t, where t = ε max(|l 0 |, |l 1 |, |l 2 |). If this test fails, we discard the point w w w ′ , and rely on a rational checker to make up the missed region (Chapter 4.4). Note that missing a region cannot break the soundness of the result. Input: poly Q : the polyhedron to be projected [x p , ..., x q ]: the variables to be eliminated n: number of initial points Output: (optimums Q ,regions Q×F ) the set of optimal functions with their regions Function Plp(poly Q , [x p , ..., x q ], n) /* convert the polyhedron to be projected which is stored in rational numbers into floaing point numbers */ poly F = GetFloat(poly Q ) // minimize the polyhedron to be projected irredundantIdx = Minimize(poly F ) minimizedPoly Q×F = GetSubPoly(poly Q , irredundantIdx) // construct the PLP problem for computing projection pl p Q×F = ConstructPlp(minimizedPoly Q×F , [x p , ..., x q ]) // take any point different from the normalization point as an initial task point

worklist F = GetInitialTasks(minimizedPoly Q , n) (optimums Q ,regions Q×F ) = none while worklist F ̸ = none do // take a task in the work list (w F , R f rom Q , F f rom ) = GetTask(worklist F ) /* R curr Q is none if w F is not inside any known region; otherwise R curr Q is the region w F located in */ R curr Q = CheckCovered(regions F , w F ) // compute
if AreAdjacent(R curr Q , R f rom Q , F f rom ) then F curr = GetCrossBoundary(R curr Q , R f rom Q , F f rom ) StoreAdjacencyInfo(R f rom Q , F f rom , R curr Q , F curr ) else AddExtraPoint(worklist, R curr Q , R f rom Q )

Checkers and Rational Solvers

We implemented most part of our procedure in floating-point numbers for the search of a solution followed by a reconstruction of the exact solution using efficient algorithm of linear algebra. But there are several cases where we should switch to exact checkers (and solvers if necessary) during the search in order to ensure the soundness and precision of the result. In the following sections we will address these rational processes.

Detecting flat regions

Our regions are obtained by reconstructing the rational matrix of the PLP objective function. They are converted into floating-point representation when the solver builds an LP problem for GLPK by instantiating the parameters in the PLP objective with a point (in floating point representation) picked up in the work list. As the regions are normalized and intersect at the same point, they are in the shape of cones. During the conversion, the constraints will lose accuracy, and thus a cone could be misjudged as flat, meaning it has empty interior. For instance, we have a cone {C 1 : -100000001 10000000 x 1 + x 2 ≤ 0, C 2 : 100000000 10000000 x 1x 2 ≤ 0}, which is not flat. After conversion, C 1 and C 2 will be represented in floating-point numbers as {C 1 : -10.0x 1 + x 2 ≤ 0, C 2 : 10.0x 1x 2 ≤ 0}, and the floating-point cone is flat. Thus the flat region will be discarded which results in imprecision of the result.

Our algorithm of raytracing minimization will report that the region is flat if it cannot find a point inside the floating point region. In this case we need to invoke a rational checker to test if the region is really flat. In The rational checker, all the constraints are shifted to the direction in the interior of the region. As the region is a cone, the constraints of a flat region will become infeasible after shifting, and a non-flat region will still be a cone. Thus we invoke our pure rational Simplex solver to test the satisfiability of the shifted constraints. If the shifted constraints becomes infeasible, i.e., the region is flat, we launch a rational minimization algorithm, which is implemented using Farkas' lemma, to obtain the minimized region.

Verifying feasibility of the result provided by GLPK

GLPK uses a threshold (10 -7 by default) when checking the feasibility of the solution. It may thus report a feasible result when the problem is in fact infeasible. Assume that we have an LP problem whose constraints are C 1 : λ 1 ≥ 0,C 2 : λ 2 ≥ 0,C 3 : λ 1 + λ 2 ≤ 10 -8 , GLPK will obtain (λ 1 , λ 2 ) = (0, 0) as a solution, whereas this point does not fulfill the constraints.

A solution is feasible if all the nonbasic variables reach their lower bound 0, and the basic variables fulfill their bound λ λ λ ≥ 0. That means when the solution is infeasible, at least one of the basic variable is negative, i.e., ∃λ i ∈ B, λ i < 0, where B is the basis provided by GLPK. In this case the objective function of the PLP is not a non-negative combination of the constraints of the polyhedron to be projected, which results in an unsound result (see Example 3.9). Therefore we must test the feasibility of the result provided by GLPK in rational numbers.

Example 3.9. Consider the example in Figure 3.15: we have a polyhedron P = {x 1 + x 2 -5 ≥ 0, -x 1 + x 2 + 3 ≥ 0, x 1x 2 + 2 ≥ 0, -x 1x 2 + 8 ≥ 0}, and we want to eliminate x 1 . The PLP problem for computing this projection is shown below.

minimize λ 1 (x 1 + x 2 -5) + λ 2 (-x 1 + x 2 + 3) + λ 3 (x 1 -x 2 + 2) + λ 4 (-x 1 -x 2 + 8) + λ 0 subject to λ 1 + 3λ 2 + 2λ 3 + 2λ 4 + λ 0 = 1 λ 1 -λ 2 + λ 3 -λ 4 = 0 and λ i ≥ 0 (∀i ∈ {0, • • • , 4}) (3.20)
The expected projected polyhedron is P ′ = {x 2 ≥ 1, x 2 ≤ 5}. If we choose λ 1 , λ 3 as basic variables, we will have λ 1 = 5λ 2 + 4λ 4 + λ 0 -1, and λ 3 = -4λ 2 -3λ 4λ 0 + 1. By setting λ 2 , λ 4 , λ 0 to 0, we obtain λ 1 = -1, and thus the basis λ 1 , λ 3 leads the LP problem to infeasibility. If we use the infeasible basis λ 1 , λ 3 to reconstruct the objective function, we will obtain a new objective function:

λ 2 (10x 2 -30) + λ 4 (6x 2 -18) + λ 0 (2x 2 -7) + (-2x 2 + 7)
from which we can obtain a face x 2 ≤ 7 2 . This face narrows the expected polyhedron P ′ to {x 2 ≥ 1, x 2 ≥ 7 2 }, and thus the obtained polyhedron is not an over-approximation of the correct polyhedron. We end up with an unsound result. We just need to verify that all coordinates in b b b ′ are non-negative. If it is not in this case, GLPK must encounter a precision problem. Then we call our Simplex algorithm in rational arithmetic which follows the standard textbook implementation. But it is not pure-rational. We construct the objective function of the LP problem with a point in floating point numbers, and thus we still use floating-point numbers for checking if all the coefficients of non-basic variables in the objective function are non-negative, i.e. checking if o j ≥ ε, where o j is the floating-point coefficient of the objective function and ε is the machine epsilon. But we do the pivoting in the Simplex tableau in rational numbers. This mixing of rational and floating point computations may provide a fake optimal solution, i.e., a solution which is feasible but not optimized, but in the context of projection and convex hull of polyhedra for static analysis it will not lead to unsoundness. A fake optimal solution may lead to less precise solution, i.e., some faces are missing, but our adjacency checker (see Section 4.4) will reconstruct all the missing faces.

Algorithm 3: Rational feasibility checker and solver.

Input: pl p Q : the PLP problem to be solved w F : the current task point basicIndices: the indices of basic variables obtained from GLPK Function FeasibilityChecker(pl p Q , basis) currOb jective F = GetObj(w F , pl p Q , feasibleBasis) optimal = IsOptimal(currOb jective F ) // compute an optimal solution while optimal == false do enteringIdx = GetEnteringVar(currOb jective F ) newBasis = RationalSimplexPivot(pl p Q , enteringIdx) currOb jective F = GetObj(w F , pl p Q , feasibleBasis) optimal = IsOptimal(currOb jective F ) return newBasis else return basis

reconstructMatrix Q = Reconstruct(pl p Q , basicIndices) feasible = IsFeasible(reconstructMatrix Q ) /* if GLPK

Obtaining sound and exact solution

We presented the rational checkers and solvers. In this section we summarize the mean for obtaining a sound and exact solution. We consider in turn all the processes (1 to 8) of our PLP solver depicted in the flowchart of Figure 3.1 and for each one we study the potential impact of floating point computations on the exact solution.

The result should always be sound if all the variables fulfill λ λ λ ≥ 0, i.e., if the PLP problem is feasible, as the objective is a non-negative combination of the constraints of the polyhedron to be projected. As we explained in Section 3.5.2, we used a rational checker (and a rational solver) to ensure the feasibility of the PLP, the result must be sound.

We have explained in Section 3.2.3 that in raytracing minimization all the constraints which cannot be determined by the floating point solver are reported as redundant, and meanwhile they are added into a undetermined list. As the regions are stored in both floating point and rational numbers, all the constraints with an undetermined status (redundant or irredundant) can be resolved using rational solver which looks for the existence of a Farkas' combination expressing the undetermined constraint. So the results of processes 1 (minimize the input polyhedron using raytracing minimization) and 8 (minimize the region using raytracing minimization) can be fixed by rational solvers and exact solutions will always be produced.

The process 2 (add initial instantiation point of parameters into the work list) does not affect the result, as any point in the parameter space can be chosen for the initial exploration. We explained in Section 3.4.4 that the process 3 (test if the point belongs to a known region) will not discard any instantiation point, as a point which is not covered by any region can never be reported to be covered. Thus the process 3 cannot lead to impression.

Considering the processes 4 (construct an LP problem with the instantiation point) and 5 (solve the LP problem with GLPK), in Section 3.5.2 we explained that GLPK may obtain a solution which is infeasible, and thus we need to verify the feasibility of the result provided by GLPK. GLPK may also report a solution which is not optimal for the given direction but optimal for a slightly different LP problem. We did not provide a checker for this situation, as in any case we will be able to use the returned basis to construct a region. This step may lead to missing regions, but they can be reconstructed by a rational procedure (Section 4.4).

The process 6 (add an instantiation point between the two regions) may discard some points in extreme cases (explained in Section 3.4.7), and this may cause the miss of regions. But the missing regions will be made up by the adjacency checker at the end, as what we presented in Section 4.4.

The process 7 (check if a region is flat in floating point numbers) has been presented in Section 3.5.1: a rational checker (and a rational solver) will be launched when the floating point solver reports that the region is flat.

As a summary, all the processes running in floating point arithmetic and having an impact on soundness and precision are verified by rational checkers and the result is repaired by rational processes when it is necessary. The adjacency checker at the end guarantees completeness of the exploration. Hence our approach is able to obtain the same solution as a pure rational approach.

Dealing with Equalities

Until now all the constraints we mentioned are inequalities. In this section we will talk about how to deal with equalities.

Minimization with equalities

Initially the equalities and inequalities are divided and stored separately. We choose a variable in each equality and substitute it in the inequalities and equalities by its definition in terms of other variables. Then replace these selected variables in the inequalities by the others.

Example 3.10. Given the polyhedron

P = {3x 1 + x 2 + 2x 4 + 4x 5 ≥ 5, -x 1 + x 3 + 4x 4 + 2x 5 ≥ -2, 6x 1 -2x 3 = 5, 8x 1 -3x 2 = 2, -2x 1 + 3x 2 -2x 3 = 3}, the constraints are divided into P eqn = {6x 1 -2x 3 = 5, 8x 1 -3x 2 = 2, -2x 1 + 3x 2 -2x 3 = 3} P ineq = {3x 1 + x 2 + 2x 4 + 4x 5 ≥ 5, -x 1 + x 3 + 4x 4 + 2x 5 ≥ -2}
As x 1 = 1 3 x 3 + 5 6 , we can obtain

x 1 = 1 3 x 3 + 5 6 -3x 2 + 8 3 x 3 = - 14 3 3x 2 - 8 3 x 3 = 14 3
As x 2 = 8 9 x 3 + 14 9 , we have

x 1 = 1 3 x 3 + 5 6
x 2 = 8 9

x 3 + 14 9

(3.21)

The third equality becomes 14 3 = 14 3 , so it is eliminated. The Equation 3.21 shows the equalities without redundancy. Now we substitute x 1 and x 2 with x 3 in the inequalities. The obtained constraints are shown in Equation 3.22.

34x 3 + 36x 4 + 72x 5 ≥ 17 4x 3 + 24x 4 + 12x 5 ≥ -7 (3.22)
Then we minimize the inequalities using either method presented in Section 3.2. The minimized polyhedron is:

P ′ = {6x 1 -2x 3 = 5, 9x 2 -8x 3 = 14, 34x 3 + 36x 4 + 72x 5 ≥ 17, 4x 3 + 24x 4 + 12x 5 ≥ -7}.

Projection of a polyhedron with equalities

Before computing projection of the polyhedron, we select variables defined by the equalities among the ones to be eliminated. Then we apply the projection operators to the equalities and inequalities separately. The variable x e to be eliminated no more appears in the inequalities and its value is defined by an equality in terms of other variables. Thus this equality does not restrict the projected polyhedron and it can simply be removed.

Example 3.11. Consider the polyhedron P = {x 1x 2 = 0, x 1 + x 3 ≥ 5}. If we want to eliminate x 1 , then the polyhedron will be transformed into P = {x 1 = x 2 , x 2 + x 3 ≥ 5}. The variable x 1 no more appears in the inequalities, and only in equality x 1 = x 2 . The projected polyhedron is then P\{x 1 } = {x 2 + x 3 ≥ 5}.

Convex hull of polyhedra with equalities

To compute the convex hull of two polyhedra using PLP, we need to pick up an interior point for minimization in one of the polyhedra. In the case where both polyhedra contain equalities (thus they both have empty interior), it is impossible to find such an interior point.

One strategy to solve this problem is that instead of selecting a point inside one of the polyhedra, we can find a point in the interior of their convex hull, as long as the convex hull has non-empty interior. By the definition of the convex hull

P 1 ⊔ P 2 = {x x x|x x x = αx x x ′ + (1 -α)x x x ′′ , x x x ′ ∈ P 1 , x x x ′′ ∈ P 2 , 0 ≤ α ≤ 1},
we can choose a point x x x ′ inside P 1 and another point x x x ′′ inside P 2 for computing a point x x x inside their convex hull.

Another strategy to deal with this situation is to compute the convex hull via projection. Benoy et al. presented in [START_REF] Benoy | Computing convex hulls with a linear solver[END_REF] a method to get the convex hull by computing a polyhedral projection. Recall the convex hull of two polyhedra P 1 and P 2 :

P 1 ⊔ P 2 = {x x x|x x x = αx x x ′ + (1 -α)x x x ′′ , A ′ I x x x ′ + b b b ′ I ≥ 0, A ′ E x x x ′ + b b b ′ E = 0, A ′′ I x x x ′′ + b b b ′′ I ≥ 0, A ′′ E x x x ′′ + b b b ′′ E = 0, 0 ≤ α ≤ 1} (3.23)
To obtain the convex hull, we need to eliminate x x x, x x x ′ and α. As the equation x x x = αx x x ′ + (1α)x x x ′′ is not linear, we need to change variables. By introducing y y y ′ = αx x x ′ and y y y ′′ = (1α)x x x ′′ , we obtain Equation 3.24.

P 1 ⊔ P 2 = {x x x|x x x = y y y ′ + y y y ′′ , A ′ I y y y ′ + αb b b ′ I ≥ 0, A ′ E y y y ′ + αb b b ′ E = 0, A ′′ I y y y ′′ + (1 -α)b b b ′′ I ≥ 0, A ′′ E y y y ′′ + (1 -α)b b b ′′ E = 0, 0 ≤ α ≤ 1} (3.24)
By substituting y y y ′′ = x x xy y y ′ into the inequalities, we obtain the encoding of convex hull which is shown in Equation 3.25.

P 1 ⊔ P 2 = {x x x|A ′ I y y y ′ + αb b b ′ I ≥ 0, A ′ E y y y ′ + αb b b ′ E = 0, A ′′ I (x x x -y y y ′ ) + (1 -α)b b b ′′ I ≥ 0, A ′′ E (x x x -y y y ′ ) + (1 -α)b b b ′′ E = 0, 0 ≤ α ≤ 1} (3.25)
To obtain the polyhedron which represents the convex hull, i.e. the set of x x x, we need to eliminate y y y ′ and α.

makes the corresponding basic variable equal to 0. In this situation the vertex could be specified in 4 different ways, i.e., we may have 4 bases if any three constraints of the four could be chosen as non-basic variables. The optimal vertex remains the same, but the partition of basic and nonbasic variables varies. This is a case of primal degeneracy. Let us take Example 4.1 to figure out how primal degeneracy occurs during computations. We will see that given a entering variable, there could be multiple choices of leaving variables as their ratios (the constant over the coefficient) are the same. We choose one of them to leave the basis, and in the next pivoting the constants corresponding to the remaining candidates will become zero.

= {x 3 ≥ 0, 2x 1 -x 3 ≥ 2, -2x 1 + 4x 2 -3x 3 ≥ 2, -2x 1 -x 3 ≥ -6, -2x 2 -x 3 ≥ -8}
, and we compute the maximal value in the direction (0, 0, 1). The LP problem is shown in Problem 4.1.

maximize x 3 subject to 2x 1 -x 3 ≥ 2 -2x 1 + 4x 2 -3x 3 ≥ 2 -2x 1 -x 3 ≥ -6 -2x 2 -x 3 ≥ -8 and x 1 , x 2 , x 3 ≥ 0 (4.1)
We transform this problem into the standard form by adding slack variables, and obtain Problem 4.2.

maximize x 3 subject to -2x 1 + x 3 + x 4 = -2 2x 1 -4x 2 + 3x 3 + x 5 = -2 2x 1 + x 3 + x 6 = 6 2x 2 + x 3 + x 7 = 8 and x i ≥ 0 ∀i ∈ {1, 2, 3, 4, 5, 6, 7} (4.2)
We solve this LP problem by Simplex algorithm. Let us use the slack variables x 4 , x 5 , x 6 , x 7 as the initial basic variables. We can see that the problem is infeasible currently, and we need to find a feasible dictionary by applying the first phase of Simplex algorithm.

After several pivoting, we obtained the feasible dictionary with basic variable x 1 , x 2 , x 6 , x 7 , and the corresponding Simplex tableau is shown in Table 4.1. Now we start the second phase of the Simplex algorithm. We choose x 3 as the entering variable, as its coefficient in the objective function is positive. We have two choice of leaving variable, as both row 3 and row 4 bound x 3 by 2. The primal degeneracy occurs as we have multiple choices of leaving variable. We choose x 6 to leave the basis. By doing one pivoting, we obtain Table 4.2. The obtained objective function does not contain positive coefficient, so the problem is maximized, and the optimal value is 2. The basic variables are x 1 , x 2 , x 3 , x 7 . By setting the nonbasic variables to 0, we obtain the value of basic variables, and the optimal vertex is (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 ) = (2, 3, 2, 0, 0, 0, 0). Geometrically the optimal point is (2, 3, 2) (Figure 4.2(a)). In Table 4.2 x 7 is a basic variable, but its value is 0. As we said, the degeneracy occurs when there are basic variables have value 0. Indeed exchanging the basic variable x 7 with any of the nonbasic variables x 4 , x 5 , x 6 will bound x 8 to the same value 0, and thus the value of the objective function will remain the same.

If we choose x 7 as the leaving variable instead of x 6 , the optimal value and the optimal vertex will be the same, but the basis will become x 1 , x 2 , x 3 , x 6 . If we do one more pivoting, another basis x 1 , x 2 , x 3 , x 5 at the same optimal vertex will be found. Another possible construction of the basis is x 1 , x 2 , x 3 , x 4 . We force x 4 to enter the basis, and obtain Table 4.3. However, although it leads to the optimal value 2, the objective function contains a positive coefficient. Thus the Simplex algorithm believes that the LP did not reach the optimum. We obtained 3 possible bases which lead the LP problem to the optimum.

x 1 x 2 x 3 x 4 x 5 x 6 x 7 = constant row 1 1 0 -1 2 -1 2 0 0 0 1 row 2 0 1 -1 -1 4 -1 4 0 0 1 row 3 0 0 2 1 0 1 0 4 row 4 0 0 3 1 2 1 2 0 1 6 objective 0 0 1 0 0 0 0 0 Table 4.1 x 1 x 2 x 3 x 4 x 5 x 6 x 7 = constant row 1 1 0 0 -1 4 0 1 4 0 2 row 2 0 1 0 1 4 -1 4 1 2 0 3 row 3 0 0 1 1 2 0 1 2 0 2 row 4 0 0 0 -1 1 2 -3 2 1 0 objective 0 0 0 -1 2 0 -1 2 0 -2 Table 4.2 x 1 x 2 x 3 x 4 x 5 x 6 x 7 = constant row 1 1 0 0 0 -1 8 5 8 -1 4 2 row 2 0 1 0 0 -1 8 1 8 1 4 3 row 3 0 0 1 0 1 4 -1 4 1 2 2 row 4 0 0 0 1 -1 2 3 2 -1 0 objective 0 0 0 0 -1 4 1 4 -1 2 -2 Table 4.3

Dual degeneracy

When dual degeneracy occurs, the objective function reaches the same optimum at multiple vertices.

The optimal value remains the same, but the optimal vertices are different. Let us illustrate the dual degeneracy on Example 4.2.

Example 4.2. Considering the polyhedron in Figure 4.1, let us compute the maximization in the direction -x 3 . By applying the first phase and several pivoting in the second phase of the Simplex algorithm, we obtain Table 4.4, in which the basis is x 1 , x 2 , x 6 , x 7 , and the corresponding optimal vertex is (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 ) = (1, 1, 0, 0, 0, 0, 4, 6). We could also obtain three other optimal bases: x 1 , x 2 , x 4 , x 7 with optimal vertex (3, 2, 0, 2, 0, -2, 4), x 1 , x 2 , x 4 , x 5 with optimal vertex (3, 4, 0, 4, 8, 0, 0) and x 1 , x 2 , x 5 , x 6 with optimal vertex (1, 4, 0, 0, 12, 4, 0). In all these situations, the optimal value are the same which is 0, but the vertices are different. That is the difference from the primal degeneracy.

Geometrically we could obtain the four vertices at the bottom in Figure 4.1, which are (1, 1, 0), (3, 2, 0), (3, 4, 0) and (1, 4, 0) (Figure 4.2(b)).

x 1 x 2 x 3 x 4 x 5 x 6 x 7 = constant row 1 1 0 -1 2 -1 2 0 0 0 1 row 2 0 1 -1 -1 4 -1 4 0 0 1 row 3 0 0 2 1 0 1 0 4 row 4 0 0 3 1 2 1 2 0 1 6 objective 0 0 -1 0 0 0 0 0 Table 4.

Degeneracy in linear programming

The problem of degeneracy in LP has been known for decades. In this section we talk about the primal degeneracy firstly, and then the dual degeneracy. Several methods to deal with the primal degeneracy for avoiding cycling will be introduced.

Primal degeneracy

The occurrence of primal degeneracy leads to inefficiency of the algorithm, and moreover it may cause the Simplex algorithm to cycle. Many methods have been proposed to avoid cycling, for example Bland's rule [START_REF] Robert | New finite pivoting rules for the simplex method[END_REF] and the perturbation method [START_REF] George B Dantzig | Application of the simplex method to a transportation problem[END_REF].

Bland's rule guarantees the termination of the Simplex algorithm by selecting the entering and leaving variables which have the smallest indices. Bland's rule is easy to implement but the value of objective function increases slowly.

Theorem 4.1 (Bland's rule). Among all candidates to enter the basis, select the variable having the lowest index; among all candidates to leave the basis, select the variable having the lowest index.

Another well-known method is the perturbation method. To avoid multiple choices of the leaving variable, this method perturbs the constraints by adding a small value ε i to each constraint C i and make sure that

ε 1 ≫ ε 2 ≫ • • • ≫ ε n-1 ≫ ε n > 0.
But it is difficult to choose the value of ε ε ε, and there will be arithmetic problem. A common way to choose the value for ε ε ε is that

ε 1 = c, ε 2 = c 2 , • • • , ε n = c n ,
where c is a very small positive constant. But it still cannot ensure that ε 1 > kε 2 after amount of additions, where k is a large positive constant. A better way to implement the perturbation method is the lexicographic method [START_REF] George B Dantzig | The generalized simplex method for minimizing a linear form under linear inequality restraints[END_REF], in which ε ε ε are kept symbolic. Let us present the lexicographic perturbation with an example.

x 3 subject to -2x 1 + x 3 + x 4 = -2 + ε 1 2x 1 -4x 2 + 3x 3 + x 5 = -2 + ε 2 2x 1 + x 3 + x 6 = 6 + ε 3 2x 2 + x 3 + x 7 = 8 + ε 4 and x i ≥ 0 ∀i ∈ {1, 2, 3, 4, 5, 6, 7} (4.3)
After perturbation, the Table 4.1 (in Page 62) becomes Table 4.5. Given the entering variable x 3 , the only choice of leaving variable is x 7 , as 2

+ 1 2 ε 1 + 1 2 ε 3 > 2 + 1 6 ε 1 + 1 6 ε 2 + 1 3 ε 4 .
In other words, the dictionary in Table 4.2 will become infeasible as we can obtain

x 7 = -ε 1 + 1 2 ε 2 -3 2 ε 3 < 0.
For the same reason, the basis x 1 , x 2 , x 3 , x 4 will also lead to infeasibility.

Geometrically the vertex v v v is split into two vertices v v v 1 and v v v 2 , as it is shown in Figure 4.3.

v v v 1 = (2 -5 12 ε 1 + 1 12 ε 2 + 1 6 ε 4 , 3 -1 12 ε 1 -1 12 ε 2 + 1 3 ε 4 , 2 + 1 6 ε 1 + 1 6 ε 2 + 1 3 ε 4 )
is the optimal vertex without primal degeneracy.

x 1 x 2 x 3 x 4 x 5 x 6 x 7 = constant row 1 1 0 -1 2 -1 2 0 0 0 1 -1 2 ε 1 row 2 0 1 -1 -1 4 -1 4 0 0 1 -1 4 ε 1 -1 4 ε 2 row 3 0 0 2 1 0 1 0 4 + ε 1 + ε 3 row 4 0 0 3 1 2 1 2 0 1 6 + 1 2 ε 1 + 1 2 ε 2 + ε 4 objective 0 0 1 0 0 0 0 0 Table 4.5
An alternative way to implement the perturbation method is to compare the perturbation terms in lexicographic order, i.e., the lexicographic perturbation, which will be explained in details later in this chapter.

v v1 v2

Dual degeneracy

It is proposed in [START_REF] Horst | Linear programming and its applications[END_REF] that: cycling can only occur in the presence of primal degeneracy. That means the dual degeneracy will not lead to cycling. When there are multiple optimal vertices, the Simplex algorithm will reach one of them. If we only care about the optimal value the dual degeneracy will not affect the result.

Degeneracy in Parametric Linear Programming

In this section we show the impact of the degeneracy on the PLP. We will first prove that there is no dual degeneracy our PLP problems for computing projection and convex hull thanks to the normalization constraint. Then we will propose a method to deal with the primal degeneracy for avoiding overlapping regions.

Overlapping regions

Ideally, the parametric linear programming outputs a partition of the space of parameters, meaning that the produced regions do not have overlap except on their boundaries (we shall from now on say "do not overlap" for short) and cover the whole parameter space. This may not be the case due to two reasons: i) geometric degeneracy, which leads to overlapping regions; ii) imprecision due to floating point arithmetic, which result in insufficient coverage. The latter will be dealt with by our rational checker, which will be explained in Section 3.5.

When the regions do not overlap, it is possible to check if the space of parameters is fully covered. It needs to verify that each boundary of a region is shared by one of its adjacent regions (proof in Section 3.5). If a boundary has no region on the other side, the space is not fully covered. This simple test cannot be used if there exist overlapping regions. We thus need to get rid of overlapping regions.

In a non-degenerate parametric linear program, for a given optimization objective, there is only one optimal vertex (no dual degeneracy), and this optimal vertex is described by only one basis (no primal degeneracy), i.e., the optimal solution corresponds to a unique partition of variables into basic and non-basic. Thus, in a non-degenerate parametric linear program, for a given instantiation of parameters there is one single optimal basis (except at boundaries), meaning that only one basis corresponds to one region. However when there is degeneracy, there will be multiple bases corresponding to one optimal vertex/function, and each of them computes a region. These regions may be overlapping. We call the regions corresponds to the same optimal function degenerated regions. Theorem 4.2. There will not be overlapping regions if there is no degeneracy.

Proof. Recall that the solution of a PLP problem is a set of tuples (R i , Z * i (x x x)), each of which associates to a basis obtained from GLPK. Each optimal function may correspond to multiple regions, i.e., ∃i, j such that

R i ̸ = R j but Z * i (x x x) = Z * j (x x x).
In parametric linear programming, the regions are yielded by the partition of variables into basic and non-basic, i.e., each region is defined by a unique basis. The parameters within one region lead the PLP problem to the same partition of variables. If there are overlapping regions, say R i and R j , the PLP problem will be optimized by multiple bases when the parameters belong to R i ⊓ R j . In this case there must be degeneracy: these multiple bases may lead to multiple optimal vertices when we have dual degeneracy (in fact we will prove that there is no dual degeneracy in our algorithm of PLP solver), or the same optimal vertex when we have primal degeneracy. By transposition, we know that if there is non-basic variables basic variables constants

λ 1 • • • λ q • • • • • • λ r • • • λ s • • • λ n . . .
objective f 1 • • • f q • • • • • • 0 • • • 0 • • • 0 -Z * (x x x) Table 4.6 -Simplex tableau.
no degeneracy, the PLP problem will always obtain a unique basis with given parameters, and there will be no overlapping regions.

Dual degeneracy in PLP

Theorem 4.3. For projection and convex hull, the parametric linear program exhibits no dual degeneracy.

Proof. We shall see that the normalization constraint (the constraint ( * ) in Problem 3.9) in the parametric linear programs defining projection and convex hull prevents dual degeneracy.

Assume that at the optimum Z * (x x x) we have the Simplex tableau in Table 4.6. Let λ k denote the decision variables with λ k ≥ 0, and let f k = o o o k x x x + c k be the parametric coefficients of the objective function. Assuming the basic variable leaving the basis is λ r and the nonbasic variable entering the basis is λ q . Currently λ r is defined by the jth row as ∑ j a j p λ p + λ r = b j , where the λ p are nonbasic variables. That means λ r = b j when the nonbasic variables are set to their lower bound, which is 0 here. We look for another optimum in another direction by doing one pivoting. As the current system is feasible, we have b j ≥ 0. To maintain the feasibility, we must choose λ q such that a jq > 0. As we only choose the non-basic variable who has negative objective coefficient to enter the basis, then we know f q < 0. By pivoting we obtain the new objective function Z ′ (λ λ λ , x x x) = Z(λ λ λ , x x x)f q b j a jq . The new basis defines a new solution λ * and a new optimal function shown in Equation 4.4.

Z * ′ (x x x) = Z * (x x x) -f q b j a jq (4.4)
Let us assume that we are faced with the dual degeneracy, which means that we obtain the same objective function after the pivoting, i.e., Z * ′ (x x x) = tZ * (x x x), where t is a positive scalar. Let us denote the normalization point by p p p. Due to the normalization constraint at the point p p p enforcing Z * ′ (p p p) = Z * (p p p) = 1, we have t = 1. Hence we obtain Equation 4.5.

Z * ′ (x x x) = Z * (x x x) (4.5)
Considering Equation 4.4 and Equation 4.5 we obtain f q b j a jq = 0 (4.6)

Since f q ̸ = 0, b j must equal to 0, which means that we in fact faced a primal degeneracy.

Primal degeneracy in PLP

Many methods to prevent the non-parametric LP problem from cycling are known [START_REF] Robert | New finite pivoting rules for the simplex method[END_REF][START_REF] George | Linear programming 2: theory and extensions[END_REF][START_REF] George B Dantzig | Application of the simplex method to a transportation problem[END_REF]. These methods can be applied to the PLP problems. But in PLP problems we also need to avoid overlapping regions besides avoiding the cycling. We implemented an approach to avoid overlapping regions based on the work of Jones et al. [START_REF] Jones | Lexicographic perturbation for multiparametric linear programming with applications to control[END_REF], which used the lexicographic perturbation method [START_REF] George | Linear programming 2: theory and extensions[END_REF]. This method adds the perturbation terms to the right-hand side of the constraints. We need to modify their approach for two reasons: i) GLPK does not support the lexicographic perturbation method, and thus we cannot solve an LP problem with the perturbation matrix; ii) we need to adapt the parallel algorithm.

Our strategy is the following: once entering a new region, we check if there is primal degeneracy: it occurs when one or several basic variables equal zero. In this case we will explore all degenerated regions for the same optimal function, using, as explained below, a method avoiding overlaps. Our algorithm is shown in Algorithm 4.

How to avoid the primal degeneracy in PLP? First of all, let us see a case where the regions are overlapping.

Example 4.4. Figure 4.4 shows two regions in 3D, and p p p is the normalization point. We cut the regions and obtain a 2D view of them, which is shown in blue and pink. Figure 4.5 shows the 2D view of the regions, and these regions correspond to the same optimal function. Figure 4.5(a) and Figure 4.5(b) are two possible compositions of degenerated regions that cover the whole corresponding space. Figure 4.5(c) shows the problematic case where the regions are overlapping. The reason is that when the parameters locate in the purple area, two different bases will lead the constructed LP problem to optimum. We aim to avoid the overlap and obtain the result either in Our solution against overlaps is to avoid the primal degeneracy, as we have proved that if there is no degeneracy there will be no overlapping regions (Theorem 4.2). The method to avoid primal degeneracy is adding perturbation terms to the right-hand side of the constraints [START_REF] Jones | Lexicographic perturbation for multiparametric linear programming with applications to control[END_REF]. These perturbation terms are "infinitesimal", meaning that the right-hand side, instead of being a vector of rational scalars, becomes a matrix where the first column corresponds to the original vector b b b, the second column corresponds to the first infinitesimal ε 1 , the third column to the second infinitesimal ε 2 , etc. Instead of comparing scalar coordinates using the usual ordering on rational numbers, we compare row vectors of rationals using the lexicographic ordering. After the perturbation, there will be no primal degeneracy. As the perturbed vectors

b ′ + (k ′ 1 ε 1 , • • • , k ′ n ε n ) and b ′′ + (k ′′ 1 ε 1 , • • • , k ′′ n ε n ) in two different rows cannot be equal.
Implementation of the perturbation method The initial perturbation matrix is a k * k identity matrix:

M p = I,
where k is the number of constraints. Then the perturbation matrix will be updated during the reconstruction of the constraint matrix. After adding this perturbation matrix, the right-hand side becomes

B = [b b b|M p ]. The new constants are vectors in the form v v v i = [b i 0 • • • 1 • • • 0]. We compare the vectors by lexico-order: v v v > l v v v ′ if the first non-zero element of v v v is greater than that of v v v ′ , where > l represents lexicographically greater. For instance, v v v 1 = [1 0 0] is greater than v v v 2 = [0 1 1].
To obtain a new basis, in contrast to working with non-degeneracy regions, we do not solve the problem using floating point solver, as GLPK does not support features for perturbation method. Instead, we pivot directly on the perturbated rational matrix. non-basic variables basic variables constants

λ 1 • • • λ q • • • • • • λ r • • • λ n . . .
objective • • • • • • ∑ n k=1 o qk x k + c q • • • • • • 0 • • • 0 -Z * (x x x) Table 4.7 λ 1 • • • λ q • • • • • • λ r • • • λ n . . . . .
• • • • • • 0 • • • • • • -1 a jq (∑ n k=1 o qk x k + c q ) • • • 0 -Z * (x x x) Table 4.8
In each pivoting, one non-basic variable will be chosen for entering the basis. All the non-basic variables will be considered for finding out all the adjacent regions. Given an entering variable λ j , we select λ l in C i as the leaving variable from all the constraints in which b i = 0 such that the ratio v v v i a i j is the smallest using the comparison > l . If such a leaving variable with b i = 0 does not exist, a new optimal function will be obtained by crossing the corresponding boundary and it will not be treated by Algorithm 4, but will be computed with a task point by Algorithm 2. On the other hand if a leaving variable with b i = 0 exists, we will obtain a degenerated region (the optimal function will not change). In this case we store the basis associated to that region and we maintain a set of bases that have been explored. After each pivot, we will check if the obtained basis exists in the set of bases, and thus the already explored regions will not be recomputed. We search all the degenerated regions by crossing all the boundaries of the already explored degenerated regions. Once all the boundaries have been crossed, all the degenerated regions corresponding to the same optimal function are found.

We benefit from Algorithm 4 in another aspect: there is no need to check adjacency of the degenerated regions. As we compute an degenerated region by doing only one pivoting across a boundary, the obtained region must be adjacent to the original one, as it is shown in Theorem 4.4. Proof. Considering the Table 4.7, we want to cross the boundary ∑ n k=1 o qk x k + c q ≥ 0 of the region associated to the variable λ q , i.e., we choose λ q as entering variable. Assuming the leaving variable is λ r , by doing one pivoting we obtain the tableau in Table 4.8. The obtained region is -1 a jq (∑ n k=1 o qk x k + c q ) ≥ 0. As 1 a jq ≥ 0 (otherwise the dictionary becomes infeasible), the region is equivalent to -(∑ n k=1 o qk x k + c q ) ≥ 0, which has the same boundary of the previous region ∑ n k=1 o qk x k + c q ≥ 0 but refer to the complementary half-space, and the two regions are adjacent.

As a summary, we transform from Algorithm 2 to Algorithm 4 when there are degenerated regions, which is detected by checking the existence of b i = 0. For preventing the regions from overlapping, we need to avoid the primal degeneracy by adding perturbation terms to the right-hand side of the constraints of the PLP problem. When solving the PLP problem we do pivoting in the rational Simplex tableau instead of solving the instantiated LP problems using GLPK, as GLPK does not support the perturbation method.

Adjacency checker

We have proved that if no degeneracy occurs, there will be no overlapping regions. Thus all the regions have their adjacencies. We shall now prove that no face of the projected polyhedron will be missed if each region can find all the adjacent regions across each of its boundaries.

We will show that a situation as the one shown in Figure 4.6 cannot happen: the four regions correspond to different optimal functions. R 1 , R 2 and R 3 all found their adjacent regions, but the region of the space denoted by R 4 has not been explored. In Figure 4.6 there exist two adjacent regions for boundaries of R 4 . We here show that this is an impossible situation in a PLP problem. Proof. Consider the situation in Figure 4.7: we cross the boundary F of the region R i , and the adjacent regions related to F are R j and R k with corresponding optimal functions Z j and Z k , and Z j ̸ = Z k . We need to do one pivoting to cross F from R i . Assume λ q is the entering variable associated to F in the objective. If there are two adjacent regions, there will be two possible leaving variables, say λ r and λ s . In the Simplex algorithm we always choose the variable with the smallest ratio of the constant and the coefficient as the leaving variable. When there are two possible leaving variables, the value of these two ratios must be equal, that is b j a jq = b k a kq . In this case we face the primal degeneracy, and

R1 R2

f * (x x x) - b j a jq f q = f * (x x x) -b k
a kq f q . This is a contradictory to the assumption Z j ̸ = Z k . Hence the situation can not happen. Theorem 4.5 shows that to find out all the faces, we just need to ensure that all the regions have their adjacencies. Although we tried to pick up task points between the regions which are not adjacent, in practice there may be still missed region because of floating point arithmetic, We invoke an adjacency checker in the last phase of the algorithm. If a boundary is not shared by two regions we will use our rational solver to explore beyond that boundary. We detail that last phase. Note that the information of adjacency is saved in Algorithm 2: if the regions R i and R j are adjacent by crossing the boundaries F m of R i and F n of R j , we set the cell of (R i , F m ) and (R j , F n ) to true in the adjacency table. When the adjacency checker finds out a pair (R k , F p ) whose cell in the adjacency table contains false, we cross the boundary F p and use a rational procedure to compute the missing region and the corresponding optimal function. If the region r does not have its adjacency beyond the boundary b, we compute the missing region by crossing b. In other words, we choose the variable associated to b as the entering variable and do one pivoting for obtaining the missing region. Then we check if all the boundaries of the new obtained region have known adjacent regions. The main algorithm terminates when all the new obtained regions have complete adjacencies. Algorithm 5 summarizes this process. status: irredundant, redundant, or undetermined. Initially all constraints are undetermined. In the first phase, some constraints are determined as irredundant. In the second phase, each thread will take one undetermined constraint and check its redundancy status. The list of the state is preassigned and has fixed size, which is the same as the number of constraints. So there is no concurrency problem of access and reallocation, i.e., the container will not grow and the stored elements will not be moved. Another problem we need to take care is the race condition. In this algorithm we do not need any extra work to protect the shared data because of two reasons: i) each thread will only modify a specific cell in the vector; ii) one thread will not read other cells in the list.

Parallelism using OpenMP The most straightforward way to parallelize the algorithm is to use the OpenMP directive omp parallel for, which is the combination of omp parallel and omp for. The directive omp parallel invokes additional threads, and the original and new created threads will execute the code inside the scope marked by this directive. The directive omp for splits the work and assign them among the threads.

We can also select the strategy of scheduling among static, dynamic and guided. The static strategy assigns the iterations to the threads evenly; the dynamic assigns a given number to each thread, and the thread which finishes its work firstly will take the next group of tasks; guided is similar to dynamic, but the number of tasks in each group to be allocated will decrease exponentially.

We put the second phase of the raytracing minimization inside the omp parallel for scope using dynamic strategy, and the undetermined constraints will be assigned to the active threads.

Parallelism Using TBB TBB provides tbb::parallel_for function, which implements similar strategy to omp parallel for. It accepts a range of tasks, and the range will be partitioned into subranges. As it is shown in Algorithm 7, tbb::parallel_for calls the function operator() of a given class automatically with a divided subranges.

Parallel Parametric Linear Programming Solver

We have seen in Algorithm 2 that the algorithm of the PLP solver performs independent tasks, and thus the algorithm can be parallelized. Here we will present the parallel algorithm of the PLP solver.

Scheduling tasks

We use tbb::parallel_do to assign tasks to the threads, and use the provided feeder to add new tasks. Once one thread finishes its task, another task will be assigned.

In OpenMP version, we used the OpenMP tasking model to manage the tasks, in which a single thread creates tasks, and each task will be fetched by one thread. The tasking feature is introduced since OpenMP 3.0. Different from the omp parallel directive, in which the number of iterations is known in advance, OpenMP tasks can be dynamically created, and thus it can be used in while loops and recursive functions. must obtain the same basis. Thus it is possible to detect the repetition of computations by comparing the basis before obtaining the region. Once obtain an optimal basis from GLPK, we store it into a hash table. If the basis has already been encountered, the thread should stop and abandon the current task. As the cost of solving an LP problem is much less than the whole process of computing a region, we can nearly avoid the repeated computation.

Dealing with shared data

The optimal functions are stored in a shared list. The corresponding regions are also stored for checking weather a task point belongs to any known region and testing adjacency of two regions. The elements in the list may be accessed by multiple threads while other threads are appending new regions into the list.

In TBB version, the set of regions are stored in tbb::concurrent_vector, which is a concurrent data structure provided by TBB. Multiple threads can access the stored elements, grow the container and append new elements concurrently. In OpenMP version, we use our simple lock-free implementation 3based on an array with a large, statically defined maximal capacity, using atomic operations for increasing the current number of items.

There is a problem in both versions. The size of the container will be increased before the new item is appended, and thus not all the items in the list are ready for reading. A segmentation fault will occur if we try to fetch the new item when it is not completely appended yet. For instance, if we append a new element to the container containing n elements, whose indices are 0 to n -1, at a moment the size of the container is n + 1 while it only has n elements. Thus if we try to get the new element by its index n, the segmentation fault will occur.

To solve this problem, we use the mutex lock in C++ standard library or a lock-free implementation, in which we store the number of items which are ready for reading. If the number of ready items n ready is less than the size of the container n f ill , then the threads need to wait until all the items are ready for reading. The algorithm of concurrent appending is shown in Algorithm 9.

Our solution to against this problem is to find out the "smallest" basis among all the bases associated to the degenerated regions, meaning that the basic variables have smallest subscripts. Once this basis is found, we append it into a hash table containing the bases that have been explored. Thus when another thread obtains a task point inside a degenerated region corresponding to the same optimal function, it will also find out the same smallest basis, which is always the same one. Then it will figure out that the basis has been seen and stop the computation.

Algorithm 10 shows the algorithm to find out the ordered basis list, from which we can obtain the "smallest" basis. We firstly count the number of degenerated basic variables, i.e., the number of zeros appears on the right-hand side of the reconstructed constraints and denote this number by degNum. Then from the degenerated basic variables and the nonbasic variables select degNum variables that will construct the basis with the remaining basic variables.

Example 5.1. Given the basic variables λ 2 , λ 4 of which λ 2 = 0, and the nonbasic variables λ 1 , λ 3 , we need to choose one variable among λ 1 , λ 2 and λ 3 to construct the basis with λ 4 . The first choice is λ 1 as it has the smallest subscript. We can obtain an ordered list of candidates of bases: {(λ 1 , λ 4 ), (λ 2 , λ 4 ), (λ 3 , λ 4 )}. We firstly take (λ 1 , λ 4 ), and test if the corresponding region has non-empty interior. If the region is flat, we then take (λ 2 , λ 4 ) and so on. In TBB version we benefit from the provided interfaces and concurrent containers. One drawback of TBB is that in general cases it provides stack-like constructs, which means that the tasks are last-infirst-out 4 . The algorithm of projection with PPLP is essentially a graph search algorithm. In our case it is better to do a breadth-first search to avoid repeatedly extracting the same regions (although the recomputation has been stopped by testing the existence of the basis obtained from GLPK), whereas the last-in-first-out containers lead to depth-first search.

OpenMP is easy to use. We do not need to install any library, and only use the provided compiler directives. But we need to guarantee that the compiler we used supports OpenMP.

Sequential algorithm for minimization

We first compare the performance of three minimization methods: the minimization in NewPolka library of Apron1 , rational minimization algorithm using Farkas' lemma, and the raytracing method. We consider three parameters: the number of constraints (C), the number of variable (V ) and the redundancy ratio (R). For example, if we have 10 constraints where 3 are redundant, then R = 30%. The redundant constraints are created as combinations of irredundant constraints. All the irredundant constraints have 20% density, i.e., there are 20% coefficients of the constraints are 0.

Then we compare the raytracing minimization and the floating point minimization algorithm using Farkas' lemma for analyzing the advantage of raytracing minimization on the polyhedra having low redundancy ratio.

Number of constraints

In this experiment we used the polyhedra whose number of constraints ranges from 2 to 46. We can see in Figure 6.1 that when the number of constraints is small the performance of the three methods are similar, but the execution time of Apron increases significantly after that. This can be explained. As Apron uses double description of polyhedra, the cost depends on the number of generators. We have said that the number of generators increases exponentially when the polyhedra are hypercubes, meaning that a polyhedron in n dimension having 2 n generators. When C ′ = V + 1 the polyhedra are simplexes, where C ′ denotes the number of irredundant constraints. A simplex is a n-dimensional polyhedron whose has n + 1 generators, i.e. the number of generators increases linearly.

When C ′ > V + 1 the number of generators increases at a high growth rate. to be a simplex when C ′ and V tend to be equal, and the number of generators (see Table 6.1) will decrease if C ′ is large enough corresponding to V .

The execution time of rational minimization using Farkas' lemma increases steadily, and it performs better than Apron when the generator number is large (the middle in Figure 6.2). The line of raytracing minimization is nearly flat, which shows the advantage of floating-point arithmetic. Number of redundant constraints This experiment shows the effect of redundancy ratio, which is from 10% to 90%. In Figure 6.3 the cost of Apron is high when R < 50%, i.e., the number of irredundant constraints C ′ < 20, which is more than twice that of variables. In this case the polyhedra are hypercubes, and the cost of double description explodes. The number of generators is shown in Table 6.3. This shows the disadvantage of double description. The other methods seem to be less sensitive to the redundancy ratio compared to Apron. In higher dimension We have seen that the floating-point raytracing method has advantages with respect to the rational methods. Here (Figure 6.4) we focuses on the two floating-point methods: raytracing minimization and floating point minimization using Farkas' lemma. We did experiments on polyhedra which have more constraints in higher dimensions. All the experiments show the advantage of the raytracing method when the polyhedra contain more irredundant constraints than redundant ones. In higher dimension the raytracing method has more obvious advantages.

Results analysis

In general the raytracing method is more efficient than the rational ones. The constraint-only description is more efficient when the polyhedra are hypercubes. The raytracing method is more efficient than the floating point algorithm using Farkas' lemma when the ratio of redundancy is low. Therefore, it is better to use raytracing minimization for the polyhedra contain more irredundant constraints than redundant ones.

Parallel minimization

These experiments are performed on the server S 2 whose characteristics are described on Page 87. As it is shown in Figure 6.5 and Figure 6.6, when the size of the polyhedra is not large, i.e. C = 200 and V = 30 in our experiment, using more threads than 8 does not improve the performance, especially when the redundancy ratio is low, as there are not enough parallel tasks. Besides the libraries need to pay extra effort for arranging the needless threads. In Figure 6.7 and Figure 6.8, the size of polyhedra are increased to C = 600 and V = 50. Using a larger number of threads shows better performance. Figure 6.9 compares the performance of OpenMP and TBB version on a polyhedron generated with C = 600, V = 50 and R = 95%. We can see that the parallelisms using OpenMP and TBB have similar performance.

Projection via sequential PLP

In this section, we analyze the performance of our PLP solver on projection operators. We compare its performance with that of the NewPolka library of Apron.Since NewPolka does not exploit parallelism, we compare it to our library running with only one thread, and thus we consider our PLP solver as a sequential algorithm.

Experiments on randomly generated benchmarks

In the line charts, the execution time of NewPolka library of Apron is depicted in blue, and our solver (denoted by PLP) appears in red. To illustrate the performance benefits from the floating-point arithmetic, we turned off GLPK and always use the Simplex algorithm implemented in rational numbers (denoted by rational), whose execution time is shown by the orange lines 2 . It reveals that solving the LP problems in floating-point numbers and reconstructing the rational Simplex tableau leads to significant improvement of performance.

The generation of polyhedra is controlled by 4 parameters: number of constraints (C), number of variables (V ), projection ratio (PR) and density (D). The projection ratio is the proportion of eliminated variables: for example if we eliminate 6 variables out of 10, the projection ratio is 60%. The density represents the ratio of zero coefficients: if there are 2 zeros in 10 coefficients, the density is 20%. Through numerous experiments, we found that when the parameters C = 19,V = 8, PR = 62.5% and D = 37.5%, the execution time of PLP and Apron are similar, so we maintain three of them and make the other one varying to analyze the variation of performance.

Recall that in order to obtain the projection of a convex polyhedron P defined by constraints, Apron (and all libraries based on the same approach, including PPL) computes a generator description of P, projects it and then computes a minimized constraint description.

Projection ratio In Figure 6.10(a) we study the impact of the projection ratio. We can see that the execution time of PLP is almost the same for all the cases, whereas that of Apron changes significantly. Apron spends a lot of time on computing the generator description of each polyhedron. We plot the execution time of PLP (Figure 6.10(b)) and the number of regions (Figure 6.10(c)), which vary with the same trend. That means the cost of our approach depends on the number of regions to be explored.

The more variables are eliminated, the lower dimension the projected polyhedron has. Then the cost of Chernikova's algorithm for converting from the generators into the constraints will be less (Table ??). This explains why Apron is slow when the projection ratio is low, and becomes faster when the number of variables to be eliminated is larger.

Number of constraints

Fixing the other parameters, we increase the number of constraints from 12 to 30. The result is shown in Figure 6.11. We can see that Apron is faster than PLP when constraints are 2 The minimization is still computed in floating-point numbers. This experiment shows that in high dimensions our algorithm using PLP solver is more efficient.

Experiments on SV-COMP benchmarks

In this experiment we used the analyzer Pagai [START_REF] Henry | Pagai: A path sensitive static analyser[END_REF] and SV-COMP benchmarks [START_REF] Beyer | Automatic verification of c and java programs: Sv-comp 2019[END_REF]. We randomly selected some C programs from the category of Concurrency Safety, Software System and Reach Safety.

Then we used Pagai to analyze the program and extracted the polyhedra to be projected. We compute the projection on these polyhedra, and compare the results of our PLP algorithm on projection, NewPolka of Apron and ELINA (introduced in Chapter 1). In Table 6.4, we present the name of programs, the average time spent on projection in milliseconds, the average number of constraints (ACN), variables (AVN) and regions (AR).

PAGAI was designed to minimize the complexity of the polyhedra and thus it takes a small subset of the variables into the invariants; relationship with other variables are ignored. Thus in the benchmarks most polyhedra have small number of constraints and dimensions. As a result of our experiments on the random generated polyhedra, our approach does not have advantage in this case. As it is shown, our algorithm has advantage over Apron when the polyhedra contain more constraints and/or in higher dimension, e.g, polyhedra in ldv-linux-3.0-module-loop and ldv-linux-3.0-bluetooth, as we get rid of maintaining double description. ELINA is the most efficient on SV-COMP benchmarks.

Analysis

We conclude that our approach has remarkable advantage over Apron for projecting polyhedra having large number of constraints or/and variables; but it is not a good choice for solving problems with few constraints in low dimensions. Our sequential algorithm is less efficient than ELINA on SV-COMP benchmarks, as the number of constraints/dimensions of the polyhedra is small, in which case our approach does not have advantage.

Convex hull via sequential PLP

In this section we present several experiments on computing the convex hull of two polyhedra. We first perform the experiments on our randomly generated benchmarks, then on the SV-COMP benchmarks.

We use V to denote the number of variables and C1 and C2 to denote the number of constraints of the two polyhedra respectively. We can see in Figure 6.14 and Figure 6.15 that when the polyhedra having a small number of constraints and variables, the NewPolka library of Apron is more efficient than our approach using PLP. When the size of polyhedra is larger, for instance in Figure 6.16 and Figure 6.17 the algorithm of PLP solver has better performance. We obtain a similar result to the experiments on projection: our algorithm of PLP solver has advantage for the polyhedra having a large number of constraints and variables. Table 6.5 presents the results of experiments on the random selected programs from SV-COMP benchmarks. Our PLP algorithm is slower than both Apron and ELINA on these benchmarks, in which the polyhedra have a small number of constraints and variables.

Parallel PLP

In this section we report some experiments of our parallel algorithm of the PLP solver on projection. First we performed experiments on randomly generated benchmarks. Secondly we used the selected SV-COMP benchmarks.

6.5.1 Randomly generated benchmarks

Performance on various numbers of constraints and dimensions

In the first experiment, we consider polyhedra in very low dimensions. We project out 1 variable from polyhedra with 4 constraints and 2 variables. Figure 6.18 shows that the execution time increases slowly with the number of threads. As we sum up the running time of 10 polyhedra with the same parameters, the average running time of one polyhedron is only 0.35 to 0.7 milliseconds. Our algorithm takes time to construct the PLP and to preprocess the constraints (Gaussian elimination, eliminating implicit constraints, etc.) which are sequential. The parallelized approaches then solve the LP problems and reconstruct the rational results. When the LP problems are small, the parallel steps take a low proportion of the total execution time, and thus the parallelism cannot speed up the algorithm significantly. Besides, the average number of regions of the projected polyhedra is 2, and thus there are not enough tasks for all the threads. With the increase of the number of threads, there is extra cost for scheduling the tasks and managing the shared data. Thus the performance becomes worse when the number of threads gets larger. Now let us see the results of performance on the polyhedra having 25 constraints in 15 dimensions. We can see that in Figure 6.19, in both versions the performance basically remains the same when the number of threads is greater than 8. Our understanding is that when there are not enough tasks, increasing thread number cannot improve the performance.

The results on polyhedra with 100 constraints and 50 variables are shown in Figure 6.20. The two versions have similar performance. With 12 threads, the speedup factor reaches 4. The performance of parallelism is mainly limited by the total number of regions (or the number of tasks in the work list) and computations of degenerated regions. If the number of tasks is less than that of threads, there will be some threads waiting for tasks, and thus launching more threads will not improve performance. The degenerated regions, which are regions corresponding to the same optimal function, are computed by one thread in our approach (shown in Algorithm 4). When there are a great amount of degenerated regions, the one thread computing on them may take more time than other threads, i.e., this thread may complete its tasks much later than the others. In this case the total execution time depends on the thread computing the degenerated regions.

Generation graph of regions

Let us see the generation of regions for analyzing the impact factors of the parallelism. We perform the experiments on a polyhedron having 12 constraints in 10 dimension. The performance and speedup are shown in Figure 6.21. We can see that the performance improves while the number of threads increases from 1 to 6. Beyond that the execution time tends to stabilize. The execution time even slightly increases when the number of threads is greater than 12.

In the generation graphs in Figure 6.22 to Figure 6.33, a node represents one region. A region is generated from the region above it, and the two regions are connected by a solid or dotted line. The solid line represents the region below is computed by the task point generated by the corresponding region above. The degenerated regions are connected by dotted lines. Please note that the numbers labeling do not matter, and they are just assigned according to the generation order. What is of interest is the shape of the graph, which depends on the geometry of the regions (how many adjacent regions do they have), and on the selection of the next exploration point in the work list. The width of the graph captures the possibility to keep the threads busy. Figure 6.22 and Figure 6.23 illustrate the generation graph with 1 thread. It is likely that TBB uses a last-in-first-out strategy to manage the tasks, and OpenMP uses first-in-first-out strategy. In TBB version a region is computed with a point generated by the latest obtained region, and in OpenMP version the threads tend to take the oldest tasks, i.e., the ones generated by the region 0.

From 1 thread to 6 threads (Figure 6.22 to Figure 6.26, Figure 6.29 to Figure 6.31), it is likely that all the threads are busy. In Figure 6.27 and Figure 6.32 the region 0 generated 11 task points, and one of the 12 threads needs to wait for a task. In the case of using more threads(Figure 6.28, Figure 6.33), there will be more threads waiting for tasks. On one hand not all the threads are used when there are not enough tasks, and on the other hand TBB/OpenMP needs to manage the shared data, concurrent data structures, etc, which takes time. The speedup of both versions are likely to be limited by the computation of degenerated regions when they use 16 threads, which are shown in Figure 6.

SV-COMP benchmarks

We used the same selected SV-COMP benchmarks as for the experiments of the sequential algorithm for computing projection. In this experiment, we show the average execution time of 5 executions on the same polyhedron, i.e.

total running time number of polyhedron×5 . Most of the polyhedra in these benchmarks exhibit a small number of regions, as it is shown in Table 6.4, for most of them the performance becomes worse as the number of threads increases (some of them are slightly improved at first) (shown in Figure 6.34(a) to Figure 6.34(i)). These results are similar to our experiment on randomly generated polyhedra which obtained 2 projected regions on average (Figure 6.18). The parallelism on ldv-linux-3.0-bluetooth (Figure 6.34(j)) shows good performance, as the projected polyhedra on these programs have a significant number of regions.

Conclusion

We did a set of experiments on minimization, projection and convex hull. The raytracing minimization using floating point arithmetic has a big advantage over using rational numbers. When the redundancy ratio is low, raytracing minimization is more efficient. We use it in our PLP algorithm for minimizing the obtained regions, as they contain more irredundant constraints than redundant ones.

Our PLP solver improved the efficiency significantly by using floating point arithmetic instead of rational computations. Compared to NewPolka library which used double description of polyhedra, our approach are more efficient when solving problems in high dimensions. Our approach also has an advantage for solving big problems which contains a large number of constraints. Our algorithm is less efficient for computing small sized polyhedra.

For solving the problems which are not very small, the parallelism improved the performance obviously. But it is not a good strategy to use parallel algorithm for solving small problems, which means when the number of obtained regions is small, as there will be not enough tasks for all the threads and considerable time will be spent on scheduling the tasks and managing the shared data. Thus we conclude that it is a good strategy to use our PLP solver for computing the projection/convex hull of the polyhedra having a large number of constraints and/or variables. It is better to use the parallel algorithm in this case, as when the constraint number/variable number is large, the obtained polyhedron tends to have a large number of regions.

ii) transforming the problem into its dual form and apply our method to solve primal degeneracy to the dual problem.

Perturbing the objective function The perturbated objective function is shown in Equation 7.1. We can choose a small value for ε i , and then the LP problem with instantiated objective function can be solved by GLPK. But this will lead to numerical problem. Another possibility is to deal with ε ε ε as "infinitesimal" and compare the ratio in lexico-order, as what we did for perturbing the constraints of the PLP. But the disadvantage is that we cannot use GLPK to solve the LP problems and we need to add the perturbation terms to the rational matrix, as GLPK does not support the perturbation method.

There is a third method. When the dual degeneracy occurs, one or multiple coefficients of the nonbasic variables in the objective function equal to zero. We can pick up these variables, say x i , • • • x j , and optimize the obtained objective function furthermore in each direction of x i , • • • x j one by one. Each optimization is based on the result of the previous one. Then by solving several optimization problems, a unique optimal vertex will be found.

Towards a pure floating point algorithm of PLP solver

In this section we talk about the possibility of implementing a pure floating point PLP solver for solving projection and convex hull without impact on soundness.

For the situation where the exact solution is not required, it is possible to get rid of the rational approaches and obtain a sound floating point solution, i.e., the resulting polyhedron contains all the possible states that can be reached by the program in running time. We remove the rational processes from the flowchart in Figure 3.1 and obtain the flowchart of pure floating point algorithm, which is shown in Figure 7.1.

To obtain a sound result, the obtained polyhedron could lose some faces, but it should not contain any face that narrows the correct polyhedron. We have seen that in the raytracing minimization, some irredundant constraints could be removed. Thus the witness point of the corresponding constraint will never be added, and there will be a risk of missing regions, i.e. the obtained polyhedron may lose some faces. Even thought the obtained polyhedron is sound.

Two processes in the PLP solver which were computed in rational numbers need to be replaced by floating point arithmetic: checking adjacency and verifying the feasibility of the result obtained from GLPK. The criterion of checking adjacency is that: i) the regions which are not adjacent could be judged as adjacent; ii) the regions which are adjacent should never be judged as nonadjacent. If the nonadjacent regions are misjudged as adjacent, we will end up with missing the region between them, and the final result will be larger but sound. But if the adjacent regions are misjudged as nonadjacent, we will always try to add an extra task point between them, and the algorithm may loop forever.

The objective function of the PLP problem is a non-negative combination of a set of constraints, and we expected to obtain an over-approximated polyhedron. However, as what we already said, GLPK may obtain a solution which is in fact infeasible. In this case some of the variables have negative value, and thus the obtained polyhedron is not a non-negative combination anymore, which results in an unsound result. Therefore, we need to guarantee the feasibility of the solution from GLPK using a floating point checker.

Other operators in the polyhedra domain using floating point arithmetic

We saw that the results of minimization, projection and convex hull can be guaranteed to be sound using pure floating point computations. Now we talk about the inclusion and equality test.

By our raytracing method and PLP solver, we can guarantee the soundness of eliminating redundancy, meet, join and projection operators. Our approach does not implement the widening operator, so we will not talk about it. Let us see the other operators.

The inclusion operator can be implemented by minimization: considering two polyhedra P and P ′ , P ⊑ P ′ if ∀C ′ i ∈ P ′ , C ′ i is redundant with respect to P. But we cannot guarantee the soundness of the inclusion operator computed by our floating point minimization, as our strategy is to report an undetermined constraint as redundant to ensure the soundness of minimization. When we test if P is included in P ′ by our strategy, the approach may report P ⊑ P ′ when in fact P ̸ ⊑ P ′ . An analyzer using polyhedra would consider that P ′ is a sound approximation of the reachable states. Thus some states of P which are not included in P ′ could be reached at running time.

The equality test can be implemented by double inclusion test. Similarly, we may report two polyhedra are equal, which are in fact unequal, and thus it is unsound.

As a conclusion, we cannot guarantee the soundness of inclusion and equality test using floating point minimization with the current strategy. We need to explore new methods for these operators.

Example 2 . 1 .

 21 The Figure 2.1(a) shows an unbounded polyhedron, which has 4 faces and 3 of them are unbounded; the bounded polyhedron which contains 5 bounded faces is shown in Figure 2.1(b).

  (a) Unbounded polyhedra. (b) Bounded polyhedra.

Figure 2 . 1

 21 Figure 2.1

Definition 2 . 1 (

 21 Convex set). A set K ⊂ Q n is convex if and only if for each pair of distinct points a a a, b b b ∈ K, the closed segment with endpoints a a a and b b b is contained in K.

  (a) Constraint description. (b) Generator description.

Figure 2 . 2

 22 Figure 2.2

Table 2 . 2 - 2 . 4 .

 2224 floating point numbers format in IEEE 754 standard Example What will happen if we run the piece of code below? double a = 0.1 ; double b = 0.6 ; assert ( a * 6 == b ) ;

13 )

 13 Matrix representation of an LP problem Given a matrixA ∈ Q m×n and vectors b b b ∈ Q m , o o o ∈ Q n , theLP problem can also be written in matrix form, which is shown in Problem 2.14. maximize Z(λ λ λ ) = o o o T λ λ λ subject to Aλ λ λ = b b b and λ λ λ ≥ 0 (2.14)

  25 shows the PLP in matrix form. maximize Z(λ λ λ , x x x) = x x x T Oλ λ λ subject to Aλ λ λ = b b b and λ λ λ ≥ 0 (2.25)

Example 2 . 10 .

 210 Considering the same constraints with Example 2.6, the PLP problem in Problem 2.26 can be obtained by replacing the objective vector o o o T with x x x T O, where the matrix O

Figure 3 . 1 -

 31 Figure 3.1 -Flowchart of our PLP procedure

Figure 3 . 2 -

 32 Figure 3.2 -Redundant constraint

  ∃y y y = [y 1 , . . . , y m , y 0 ] ≥ 0 m ∑ j=1, j̸ =i a jk y j = a ik ∀k = {1, • • • , n} m ∑ j=1, j̸ =i b j y j + y 0 = b i (3.2) ∃y y y = [y 1 , . . . , y m , y 0 ] ≥ 0, My y y = d d d, where M

3 )

 3 Theorem 3.1 (Farkas' lemma). Let A ∈ R m×n and b b b ∈ R m . Then exactly one of the following two statements is true: • There exists an y y y ∈ R n such that Ay y y = b b b and y y y ≥ 0.

Figure 3 . 3 -Example 3 . 3 .

 3333 Figure 3.3 -The point w 1 is a witness of the irredundancy of C 1 : It satisfies all constraints but C 1 . Therefore, removing C 1 would change the polyhedron because w 1 would then belong to it.

Figure 3 . 4 -

 34 Figure 3.4 -The process of checking the redundancy status of the constraint C 1 using the second phase of raytracing minimization.

Figure 3 . 5 -Theorem 3 . 2 .

 3532 Figure 3.5 -How to obtain an interior point in different cases: In Figure (a) and Figure (c) the central points are found; for the case of Figure (b) the constraints are shifted for obtaining an interior point.

( 3 . 8 ) 3 . 4 .

 3834 Example Let use consider the polyhedron shown in Figure3.6:

  Figure 3.7(a) we shift the non-strict inequality C ′ 1 : x 2 > 2 and obtain C ′′ 1 . By solving the satisfiability problem {x 2 ≥ 2 + ε, x 1 + x 2 ≤ 7, -x 1 + x 2 ≤ 0} we get the witness point w w w. If we compute the optimum in the direction x 2 with constraints {-x 1 + x 2 ≤ 0, x 1 + x 2 ≤ 7}, as it is shown in Figure 3.7(b), we obtain a witness point w w w ′ . (a) The ray hit C 2 first and found p p p ′ . (b) The ray launched from p p p to p p p ′ hit C 3 first.

Figure 3 . 6 (

 36 Figure 3.6

Figure 3 . 7

 37 Figure 3.7

Algorithm 1 :

 1 Raytracing minimization algorithm.Input: poly F : the polyhedron to be minimized Output: the indices of the irredundant constraints and the witness point Function Minimize(poly F ) // try to obtain a central point of the polyhedron

Figure 3 . 8 -

 38 Figure 3.8 -When the gray area is extremely small, the ray may misjudge the first hit constraint as C 2 instead of C 1 .

Figure 3 . 9 -

 39 Figure 3.9 -Projecting a 3D polyhedron to 2D.

  Figure 3.10(b) there are two more regions which are in yellow and gray. The yellow region corresponds to a redundant constraint of the projected polyhedron, and the gray region corresponds to the trivial constraint 0 ≤ 1.

  (a) The regions with normalization of Problem 3.9.(b) The regions with the constraint ∑ λ i = 1.

Figure 3 . 10

 310 Figure 3.10

Figure 3 .

 3 Figure 3.11 -The convex hull of two polyhedra in 2D.

Definition 3 . 3 (

 33 Augmented matrix). Given the system of equalities Ax x x = b b b, its augmented matrix is [A|b b b].

Example 3 . 7 .

 37 Considering the matrix in reduced form.

  which is equivalent to the substitution

  , which is equivalent to the original PLP in Problem 3.11. minimize λ 0 (-

  )

  we apply this substitution to the objective matrix O to get the matrix of the new objective function O ′ : O ′ = O -ΘM. As O B = ΘM B , by applying O -ΘM the columns corresponding to the basic variables become 0.

Figure 3 . 12 -

 312 Figure 3.12 -An example of missing regions

FFigure 3 . 13 -

 313 Figure 3.13 -Limitation of syntactic adjacency checking

Figure 3 . 14 - 2 :

 3142 Figure 3.14 -Add point between regions

Figure 3 . 15 -

 315 Figure 3.15 -The expected projected polyhedron by projecting out x 1 of P is P ′ .

Figure 4 .

 4 Figure 4.1 -Degeneracy

4

 4

  (a) Primal degeneracy: one optimum, one optimal vertex with 4 possible bases. (b) Dual degeneracy: one optimum, 4 possible vertices and one basis for each of them.

Figure 4 . 2 -

 42 Figure 4.2 -Optimal vertices.

Example 4 . 3 .

 43 Considering again the example of Problem 4.2, the perturbated problem is shown in Problem 4.3.

  maximize

Figure 4 . 3 -

 43 Figure 4.3 -Split vertex.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . row j . . . . . . . . . . . a jq • • • • • • 1 • • • 0 • • • 0 b j . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . row k . . . . . . . . . . . a kq • • • • • • 0 • • • 1 • • • 0 b k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4 .

 4 5(a) or in Figure 4.5(b).

Figure 4 . 4 -

 44 Figure 4.4 -A region in 3D.

Figure 4 . 5 -

 45 Figure 4.5 -Example of overlapping regions. The regions obtained from bases B 1 and B 2 cover the whole space corresponding to their face without overlapping; similar for B 3 and B 4 . But we may obtain the bases B 1 and B 4 , which results in overlaps and an incomplete coverage of the space.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . row j . . . . . . . . . . . a jq • • • • • • 1 • • • 0 b j . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Theorem 4 . 4 .

 44 The region obtained by doing one pivoting in the Simplex tableau must be adjacent to the original region.

Figure 4 . 6 -

 46 Figure 4.6 -An impossible situation: R 1 , R 2 , R 3 have their adjacent regions, but R 4 is missing.

Figure 4 . 7 -

 47 Figure 4.7 -The region R i has two adjacent regions R j and R k , whose optimal function Z j and Z k are not equal.

Algorithm 8 :

 8 Checking bases (basicIndices, nonbasicIndices) = GlpkSolveLp(w F , pl p F ) if basicIndices is infeasible then (basicIndices, nonbasicIndices) = RationalSimplex(w F , pl p Q ) // test if the basis has been obtained for avoiding recomputation atomic existed = AddBasis(basicIndices) if existed then return else • • •

5. 3 . 4

 34 Advantages and drawbacks of two versions 5.3.4.1 TBB

Algorithm 10 :

 10 Obtain the basic variables who has smallest subscript. Input: basicIdx: the indices of basic variables nonbasicIdx: the indices of nonbasic variables Output: a list contains the bases, which are sorted in lexico order of their subscript. Function GetBasisList(basicIdx, nonbasicIdx) idxList = none remained = none Insert(idxList,nonbasicIdx) degNum = 0 // find all the basic variables which equal to 0 foreach basic variable v do if v == 0 then Insert(idxList, GetIdx()v)) degNum = degNum + 1 else Insert(remained, GetIdx(v)) /* now idxList contains the indices of nonbasic variables and that of the basic variables which equal to 0 */ Sort(idxList) /* get all possible combinations of degNum variables from idxList, and the combinations are sorted by their indices in lexico-order */ selectedList = ChooseFrom(idxList, degNum) /* construct a set of possible bases, which are sorted by their indices in lexico-order */ newBasisList = GetBasis(remained, selectedList) return newBasisList

  Figure 6.1(b) is a zoomed figure of Figure 6.1(a), and the variation of the number of generators is shown in Table 6.1.

Figure 6 . 1 -

 61 Figure 6.1 -Variation of constraints. C=[2,46], V =12, R=50%.

Figure 6 . 2 -

 62 Figure 6.2 -Variation of variables. C=40, V =[2,20], R=50%

Figure 6 . 3 -

 63 Figure 6.3 -Variation of redundancy. C=40, V =10, R=[10%,90%].

Figure 6 .

 6 5 to Figure 6.8 shows the performance of the raytracing minimization with different number(a) C=200, V =30, R=[5%,95%] (b) C=600, V =50, R=[5%,95%] (c) C=200, V =150, R=[5%,95%] (d) C=60, V =200, R=[5%,95%]

Figure 6 . 4 -

 64 Figure 6.4 -Performance of floating-point algorithms for minimization in higher dimension

Figure 6 . 5 -

 65 Figure 6.5 -Parallel raytracing minimization with OpenMP, C=200, V =30, R=[5%,95%]

Figure 6 . 6 -Figure 6 . 7 -

 6667 Figure 6.6 -Parallel raytracing minimization with TBB, C=200, V =30, R=[5%,95%]

Figure 6 . 8 -

 68 Figure 6.8 -Parallel raytracing minimization with TBB, C=600, V =50, R=[5%,95%]

Figure 6 . 9 -

 69 Figure 6.9 -Parallel raytracing minimization with different threads, C=600, V =50, R=95%

Figure 6 . 10 -

 610 Figure 6.10 -C=19, V =8, D=37.5%, PR=[25%,87.5%]

Figure 6 . 11 -

 611 Figure 6.11 -C=[12,30], V =8, D=37.5%, PR=62.5%

Figure 6 . 12 -

 612 Figure 6.12 -C=19, V =[3,15], D=37.5%, PR=62.5%

Figure 6 . 13 -

 613 Figure 6.13 -C=19, V =8, D=[12.5%,75%], PR=62.5%

Figure 6 . 14 -

 614 Figure 6.14 -CN1=4, CN2=[2,7], V =4.Figure 6.15 -CN1=4, CN2=5, V =[2,7].

Figure 6 . 15 -

 615 Figure 6.14 -CN1=4, CN2=[2,7], V =4.Figure 6.15 -CN1=4, CN2=5, V =[2,7].

Figure 6 . 16 -

 616 Figure 6.16 -CN1=15, CN2=[12,19], V =6.Figure 6.17 -CN1=20, CN2=15, V =[2,6].

Figure 6 . 17 -

 617 Figure 6.16 -CN1=15, CN2=[12,19], V =6.Figure 6.17 -CN1=20, CN2=15, V =[2,6].

Figure 6 . 18 -

 618 Figure 6.18 -Variation of redundancy, C=4, V =2, D=0%, eliminating 1 variable. The average number of regions is 2.

Figure 6 . 19 -

 619 Figure 6.19 -C=25, V =15, D=20%, eliminating 1 variable. The average number of regions is 172.5.

Figure 6 . 20 -

 620 Figure 6.20 -C=100, V =50, D=20%, eliminating 1 variable. Average region number is 2928.4.

  28 and Figure 6.33.

Figure 6 . 21 -

 621 Figure 6.21 -C=12, V =10, D=20%, eliminating 1 variable.

Figure 6 . 22 -

 622 Figure 6.22 -C=12, V =10, D=20%, running with 1 thread, TBB.

41 Figure 6 . 23 -

 41623 Figure 6.23 -C=12, V =10, D=20%, running with 1 thread, OpenMP.

Figure 6 .

 6 Figure 6.24 -C=12, V =10, D=20%, running with 2 threads, TBB.

Figure 6 . 25 -

 625 Figure 6.25 -C=12, V =10, D=20%, running with 4 threads, TBB.

  0

Figure 6 . 33 -

 633 Figure 6.33 -C=12, V =10, D=20%, running with 16 threads, OpenMP.

Figure 6

 6 Figure 6.34

. 1 )( 7

 17 x jb i + ε i )λ i (7Dual problem Given an LP problem in Problem 7.2, the corresponding dual problem is shown in Problem 7.3. The primal problem does not have dual degeneracy if its dual problem does not have primal degeneracy. Hence we can add perturbation terms to the right-hand side of the constraints of the dual problem for avoiding primal degeneracy, and thus the original problem can get rid of dual degeneracy. maximize Z(λ λ λ ) = λ j ≤ b i ∀i ∈ {1, • • • m} and λ λ λ ≥ 0 y i ≥ o i ∀i ∈ {1, • • • m} and y y y ≥ 0 (7.3)

Table 2 .

 2 7 -Simplex tableau example.

	Example 2.11. Considering the same example with Problem 2.26, By substituting the parameters x x x
	with the point (-1, 0), we obtain the LP problem in Example 2.6. As the LP problem is optimized
	with the basic variables λ 3 and λ 4 , the PLP problem also reaches optimum with the same basis.
	By substituting λ 2 and λ 3 we obtain the final objective function Z ′ (λ λ λ , x x x) = (-3x 1 -17 3 x 2 -25 3 )λ 1 +
	(6x 1 + 2x 2 + 4)λ 4 + (8x 1 + 19x 2 + 17), and the optimal solution is (R 2 , Z * 2 (x x x)).

  Equation 3.1 can be expressed in another way, which is shown in Equation 3.2. By writing it into matrix form, we have Equation 3.3. If the system of linear inequalities has a solution, the first statement of Farkas' lemma is true; otherwise the second statement is true.

			3.1) : a constraint C i is
	redundant if and only if it is a linear combination of other constraints. In other words, C i : a a a i x x x ≥ b i is
	redundant if and only if it fulfills Equation 3.1.
		m
	∃y y y = [y 1 , . . . , y m , y 0 ] ≥ 0, a a a i x x x -b i =	∑	y j (a a a j x x x -b j ) + y 0	(3.1)
		j=1, j̸ =i

  cannot find a central point, compute any point inside the polyhedron p F = GetInteriorPoint(poly F ) if p F == none then report polyhedron has empty interior terminate * if the floating point minimization cannot determine the redundancy status, set the constraint as redundant for guaranteeing soundness */ Considering the example in Figure3.8, in which the polyhedron is colored in blue, and p p p is the interior point. We want to check if C 1 is redundant. The ray r r r launched from p p p to p p p ′ hit C 1 first, and then hit C 2 .

	if cannot determine then
	SetAsRedundant(idx)
	else
	/* if the irredundancy witness point exists, the corresponding constraint is
	irredundant; otherwise it is redundant */
	if found irredundancy witness point then
	SetAsIrredundant(idx)
	else
	SetAsRedundant(idx)
	return indices of irredundant constraints with their witness point
	computation.

// launch rays from the interior point to each boundary rays F = LaunchRays(poly F , p F ) /* compute the distance from the interior point to each boundary along the direction of ray F , and set the constraint whose boundary was first met as irredundant */ foreach ray F in rays F do constraintIdx = FirstHitConstraint(poly F , ray F , p F ) SetAsIrredundant(constraintIdx) foreach constraint idx in undetermined constraints do /

  13) Solving Problem 3.13 with GLPK leads to choosing λ 1 , r 2 as basic variables. Then we can restate Problem 3.11 into an equivalent PLP with the basis λ 1 , r 2 : ( 40 41 x 2 -7 41 ) + λ 2 (x 1 -800 41 x 2 + 960 41 ) + λ 3 (x 1 -765 41 x 2 -3x 3 + 918 41 ) + λ 0 (-40 41 x 2 + 48 41

	minimize	r 1
	subject to
	and

  .1. Let M denote the matrix of constraints and O . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

		λ 1	• • •	λ m	= constant
	. . .
	row j	a j1	• • •	a jm	b j
	. . .
	objective ∑ n k=1 o 1k

. .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	row j	. . . . . . . . . . .	1 • • • • • •	1 a jq

  This because the polyhedron tends

	constraint number	2	4	6	8	10	12	14	16
	generator number 130	130	130	130	130	130	130	130
	constraint number	18	20	22	24	26	28	30	32
	generator number 130	130	130	130	474	1252	2689 4836
	constraint number	34	36	38	40	42	44	46
	generator number 8999 15725 22472 35947 50000 72448 97582
	Table 6.1 -Number of generators. C=[2,46], V=12, R=50%.

Table 6 . 3

 63

	redundancy ratio	10%	20%	30%	40%	50%	60% 70% 80% 90%
	generator number 383823 240358 135260 61426 24435 6334 827 110 110

-Number of generators. C=40, V =10, R=[10%,90%]%

Table 6 .

 6 [START_REF]OpenMP Application Programming Interface[END_REF] -Performance of PLP algorithm on projection on SV-COMP benchmarks.

	Program	Apron (ms) ELINA (ms) PLP (ms) ACN AVN	AR
	pthread-complex-buffer	0.29	0.18	0.32	3.25	3.06	2
	ldv-linux-3.0-module-loop	0.57	0.22	0.37	3.16 16.19	2
	ssh-clnt-01	0.29	0.18	0.32	3.53	2.61	1.98
	ldv-consumption-firewire	0.45	0.28	0.7	7.19	6.14	3.2
	busybox-1.22.0-head3	0.75	0.38	1.31	10.94 6.14	5.5
	ldv-linux-3.0-usb-input	0.27	0.16	0.29	3	2	2
	bitvector-gcd	0.33	0.19	0.73	5	3	4
	array-example-sorting	0.33	0.2	0.63	4.78	3.67	2.89
	ldv-linux-3.0-bluetooth	255.66	2.46	12.46	20.62 17.66 18.86
	ssh-srvr-01	0.43	0.24	1.2	5.91	4.68	4.61

Table 6 .

 6 [START_REF] Beyer | Automatic verification of c and java programs: Sv-comp 2019[END_REF] -Performance PLP algorithm on convex hull on SV-COMP benchmarks.

	Program	Apron (ms) ELINA (ms) PLP (ms)	ACN	AVN AR
	array-memsafety	0.28	0.18	1.98	2.63, 2.63	2.5	6.75
	busybox-1.22.0	0.34	0.21	0.73	3.5, 2.5	7	3.5
	bitvector	0.31	0.19	0.5	2.0, 1.5	4.5	1.5
	loops	0.28	0.19	0.9	2.88, 2.88	3.0	4.88
	ldv-linux-3.0-atm	0.43	0.26	1.36	4.5, 4.5	7.5	5.5
	loop-invgen	0.38	0.24	9.89	5.8, 5.9	3.4	21.5

  Figure 6.26 -C=12, V =10, D=20%, running with 6 threads, TBB.Figure 6.31 -C=12, V =10, D=20%, running with 6 threads, OpenMP.Figure 6.32 -C=12, V =10, D=20%, running with 12 threads, OpenMP.

																											0
							0										1					3	4		5	6	8	11	15	14	16
		1	2	3	4	5	6	9	11	14	15	17		2	7	9	12	13	17		20	18	22	36	19	23	28	26	30	29
	7	12	10	8	18	13	19	20	23	26	27	31		10	21	25	27	32				31			33	35
		25	22	21	30	16	39	33	34	36	41			24	34
																							0
		37	28	32		24								37
			1						2				3			4				5			6		7	9		8	10	11
			38			29								38
	12	14	17	18		13	16	19	20		21	22	23	25	26		34	27	28		29	30	33		32	35	15	36	37
						35								39
						40								40													24
														41													31
	Figure 6.27 -C=12, V =10, D=20%, run-
	ning with 12 threads, TBB.						Figure 6.28 -C=12, V =10, D=20%, running with 16 39
														threads, TBB.									40
																											41
							0
			1	3	4	5	6	7	8	9	10	13
																										0
		15	17						11	22
																2				1			3	4	5	6	7	8	9	10	11
									12	24
													13	14	16	17	18	20	21	23	24	26	28	29	32	31		12	41
									14	26
																								34	35	36		15
									16	28
																								39		38		19
									18	30
																										40		22
									20	32
																												25
										34	35
																												27	30
										41	37
																												33
											39
																												37
													Figure 6.30 -C=12, V =10, D=20%, running with 4 0
		Figure 6.29 -C=12, V =10, D=20%, 1 2	threads, OpenMP. 3 4			5			7	8		9	14	15
		running with 2 threads, OpenMP.					1		2	3		6		4		5
	6	10		11	12	13		16	28		20	18	26	31		19	21	22		23	24		27	34		29	30	32
																9		7	10	11		8		12
	25					35	36
																								0
																19		15	16	13		14		18
			2						4						3					5			1	6		7	8	9	10	12
																27		20	25	17		21		23
	11	14		16	18	19	20		21	22	23	25	26		27	29		30	31	32	35	37	36	13	38	39	40
	41																	26	32	22		29		30		17
																											24
																		33	36	24		34		37
																											28
																		38	40	28		41
																											33
																					31
																											34
																					35
																					39

Now distributed as part of Apron: http://apron.cri.ensmp.fr/library/

Contributors: Alexis Fouilhé, Alexandre Maréchal, Sylvain Boulmé, Hang Yu, Michaël Périn, David Monniaux.

http://elina.ethz.ch/

There is another form of PLP, in which the parameters appear on the right-hand side of the constraints.

3 (2x 1x 2 -1) +[START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF] 3 , meaning that C 4 is redundant. In this case the first statement of Farkas' lemma is true.

https://www.gnu.org/software/glpk/

To adapt to the description in our approach, we define the reduced echelon from differently from that in the textbook.

We use Flint, which provides exact rational scalar, vector and matrix computations, including solving of linear systems. http://www.flintlib.org/

https://www.threadingbuildingblocks.org/

The OpenMP version is implemented by Camile Coti and David Monniaux.

The lock-free data structure is provided by David Monniaux.

If a thread has enough tasks to do, it will do depth-first execution; otherwise it will do temporary breadth-first execution. More details are presented in https://software.intel.com/en-us/node/506103.

https://github.com/antoinemine/apron

https://www.threadingbuildingblocks.org/

Acknowledgements

First I would like to express my special thanks to my supervisors Dr. David Monniaux and Dr. Michaël Périn, who spare no effort to help me and inspire me during the three years. They are always so patient with my buggy code and the reports having annoying structures. They not only teach me the knowledge, but also show me the method to present the work and the way to do research. I also want to present my great thank to Dr. Laurent Mounier, who helped me so much with my first attempt to teach. Then I want to show my special thank to my families, who always support me at any time, and my two cats, who always sit beside my laptop and try to type some sentences when I was writing my thesis in the night. I would like to thank all my friends and colleague in Verimag for their help and kindness.

Summary of the chapter

In this chapter we presented our sequential algorithm of PLP solver with two minimization methods. We mainly used the floating point raytracing minimization. When the floating point computation cannot determine the redundancy status of any constraint, the rational minimization implemented using Farkas' lemma will be invoked. Because of the rounding errors of the floating point arithmetic, we did some adaptation to the raytracing minimization for guaranteeing termination of the algorithm and soundness of the result, i.e., the resulting polyhedron is an over-approximation of the correct one.

We illustrated the polyhedral projection and convex hull geometrically, and also explained the way to compute them using PLP. We used GLPK to solve the LP problems in floating point numbers. To reconcile our approach with GLPK, we need to simplify the equality constraints of PLP into inequalities. In each step of the PLP solver, we must consider the impact of floating point arithmetic. Due to the floating point computations, the results may be unsound or imprecise. The result is unsound if GLPK obtain a reports an infeasible solution as a feasible one. In other words, GLPK obtain an infeasible basis but returns it as a feasible one, which means that the real feasible basis has not been found. Thus we provided a checker for testing the feasibility of the result obtained from GLPK. The checker for guaranteeing the precision will be presented in Section 4.4. The focus of our PLP solver is on computing the polyhedra with inequalities. As we explained, the equalities need to be dealt with separately.

Chapter 4

Degeneracy and Overlapping Regions

Degeneracy [START_REF] George | Linear programming 1: introduction[END_REF][START_REF] George | Linear programming 2: theory and extensions[END_REF][START_REF] Chvatal | Linear programming[END_REF][START_REF] Gal | Postoptimal Analyses, Parametric Programming, and Related Topics: degeneracy, multicriteria decision making, redundancy[END_REF] is a difficulty of our approach. It is a general issue in LP which may lead to the cycling of the Simplex algorithm, and there are mature solutions to deal with it. In PLP it leads to overlapping regions. In this chapter we present the fundamental of degeneracy, the impact on LP and PLP, and the methods to handle it.

Base on the work in [START_REF] Jones | Lexicographic perturbation for multiparametric linear programming with applications to control[END_REF] we implement an algorithm to deal with the degeneracy in our PLP solver, and thus the overlapping regions can be avoided. Once the overlapping regions are avoided, we can check the precision of the resulting polyhedron. If there are missing regions, it means that the floating point computation is not precise enough, and we use a rational procedure to compute these regions.

Introduction to degeneracy

There are two kinds of degeneracy. The primal degeneracy [START_REF] George | Linear programming 1: introduction[END_REF][START_REF] Chvatal | Linear programming[END_REF] means that the unique optimum can be expressed by several different bases, and the dual degeneracy [START_REF] George | Linear programming 2: theory and extensions[END_REF][START_REF] Gal | Postoptimal Analyses, Parametric Programming, and Related Topics: degeneracy, multicriteria decision making, redundancy[END_REF] happens when there exist more than one optimums. Now let us consider the example shown in Figure 4.1. The pyramid in 3 dimensions is made of 5 constraints and 4 of them intersect at the top point.

Primal degeneracy

Let us compute the optimum in the direction toward the top of the pyramid shown in Figure 4.1: the optimal solution is given by the top vertex.

Once the system is rewritten as equality constraints by introducing slack variables. The number of basic variables, defined in terms of the others, corresponds to the number of equalities. The original variables are nonbasic variables, and thus the number of non-basic variables equals the number of dimensions of the constraints m i=1 a a a i x x x ≤ b i , which is 3 in this example. A vertex in 3 dimensions is specified by 3 faces. That means 3 of the constraints associated to the nonbasic variables intersect at the optimal vertex, and the non-basic variables reach their lower bound (for minimization problems), which is 0 in our case. But in this example, a fourth constraint also intersects at the same vertex, which Algorithm 4: Algorithm to avoid overlapping regions.

Input: w F : the task point pl p Q×F : the PLP problem to be solved Output: degeneracy regions correspond to the same optimal solution Function DiscoverNewRegion(w F , pl p Q×F ) // get the basis which leads the LP problem to the optimum glpkBasis = GlpkSolveLp(w F , pl p F ) basicIdx = FeasibilityChecker(pl p Q , glpkBasis) /* if the adjacent region has not been found, compute the missing region using a rational procedure */ if adj == false then

Summary of the chapter

In this chapter we illustrated the primal and dual degeneracy with examples. For LP problems the primal degeneracy is more noticed than dual degeneracy as the former may lead to cycling. There are various methods of avoiding cycling, such as Blind's rule and the perturbation method.

For PLP the degeneracy leads to overlapping regions. We proved the absence of dual degeneracy in our PLP problems for computing projection and convex hull, and thus we only need to deal with the primal degeneracy. Based on the work in [START_REF] Jones | Lexicographic perturbation for multiparametric linear programming with applications to control[END_REF] we adapted the lexicographic perturbation to our PLP solver. After the perturbation we can always obtain a unique basis with any task point in each region, and thus overlapping regions are avoided.

We proved that, without overlapping regions, if all the regions have their adjacent regions, there is no missing region. If any adjacent region is missing, we use a rational procedure to compute this region. The adjacency checker guarantees the precision of our result.

Chapter 5 Parallelism

In this chapter we first present some preliminaries of our parallel algorithm. Then we show the parallel algorithm of the PLP solver with a second level parallelism for the raytracing minimization. In the parallel raytracing minimization, each thread computes the redundancy status of one constraint. The parallelism is simple to implement as the computations of the redundancy status of the constraints are independent. In the parallel PLP solver, each thread computes one optimal function with the corresponding region. It is more complex to parallelize this algorithm as the number of regions is unknown in advance, i.e., the tasks are dynamically added. Besides there is data shared by multiple threads. To implement the parallelism we need to modify our sequential algorithm shown in Algorithm 2 and Algorithm 4 at several steps. The work in this chapter was presented in our paper [START_REF] Coti | Parallel parametric linear programming solving, and application to polyhedral computations[END_REF].

Preliminaries

We will explain race conditions, which may lead the parallel algorithm to an undesirable situation, and the techniques to avoid it. Then TBB and OpenMP will be introduced, using which we implemented the parallelism.

Race conditions and thread safe

A race condition occurs when the result of an algorithm depends on the sequence of the accessing of threads. In other words, if the sequence of multiple threads is not controlled, because of the race condition, the result may be incorrect.

Let us see a simple example. Assuming we have two pieces of code in Algorithm 6. The value of val2 could be 5, if the thread A increases val after the thread B reading it, or 6 otherwise.

Thread safe means that the data types or the piece of code can be executed by multiple threads without race conditions.

Mutex and lock-free algorithm

Mutex (Mutual exclusion) is a strategy to prevent race conditions. If a piece of code is blocked by a mutex, only one (or a given number) thread can access it. By locking the shared data, the race conditions can be avoided.

Lock-free data structures can be accessed by multiple threads safely without using mutex locks. Lock-free strategy could provide better performance than using mutex.

TBB and OpenMP

We implemented the parallel algorithm using Intel's Thread Building Blocks (TBB) 1 , and OpenMP tasking 2 [START_REF]OpenMP Application Programming Interface[END_REF] respectively, which can be selected during compilation time.

Intel Threading Building Blocks (TBB) is a C++ template library for task-based parallel programming. It provides concurrent data structures and task scheduling algorithms. To use TBB, we just need to specify the tasks without any details about threads, and the TBB library will map the tasks onto threads efficiently. Besides, TBB also supports nested parallelism.

OpenMP (Open Multi-Processing) is an application program interface (API) that is used for shared memory multiprocessing programming in C/C++ and Fortran on multi-platform. It includes both compiler directives and runtime library routines. Parameters such as the number of threads to be used can be set both programmatically and through environment variables.

Parallel raytracing minimization

Now we explain the parallelism on raytracing minimization. As the first phase of the algorithm only contains a matrix multiplication and selecting the minimum value which is the minimum distance from the interior point to the constraint, the execution time can be negligible with respect to the second phase. Thus we only parallelized the second phase. The parallelism is easy to be implemented as the rays are launched independently.

The only shared data is the list containing the redundancy status of the constraints, which is stored in std::vector<Status>, where Status is a user-defined structure. Any constraint has one of the three

Dealing with multiple threads

As we explained, we do not need to take care of the assignment of tasks, as it will be scheduled by TBB or OpenMP. The only matter we need to pay attention to is avoiding multiple threads computing on the same problem, which will result in the inefficiency of the algorithm. Consider that when the thread A is computing on the region R 1 , the thread B completed the computation on R 2 and generated a witness point w w w which is inside R 1 . Then the thread C can fetch w w w and compute on the same region as the thread A, as R 1 has not been generated.

It is impossible to completely avoid multiple threads computing on the same region, but we can reduce the possibility of its occurrence by Algorithm 8. If two threads compute on the same region, they Algorithm 9: Current appending Input: (newOptimum Q ,newRegion Q×F ): the new optimal function with its region that will be appended (optimums Q ,regions Q×F ): the list of optimal functions with their regions Function CurrentAppend

while n ready < i do spinning or using a condition variable // n ready == i here atomic { n ready = i+1 }

Updating algorithm of preventing degeneracy

In sequential version of Algorithm 4, we simply start from any basis and find out all the degenerated regions which correspond to the same optimal function. But this is not enough for avoid overlapping regions in parallel algorithm. Consider the situation in Figure 5.1: thread 1 took the task point w w w 1 , and it is computing on the region R 1 ; meanwhile thread 2 obtained the task point w w w 2 , which is not covered by any known region as the computation of R 1 has not completed. Thus thread 2 obtains the basis corresponding to the region R 2 by solving the LP problem. As the obtained basis has not been found, thread 2 will continue the computation and obtain R 2 , which overlaps with R 1 . Thus we need to update Algorithm 4 for avoiding this situation, and the updated algorithm is shown in Algorithm 11.

Please note that in the sequential version this problem does not appear because the single thread needs to complete the computation on R 1 before fetching the task point w w w 2 . Then it will report that w w w 2 is covered by R 1 . // find the smallest basis in lexico-order

if newRegion Q has non-empty interior then basicIdx = b break // start from the smallest basis, compute all the degenerated regions size = GetSize(basicIdx)

// the remaining part is the same as the sequential algorithm

Summary of the chapter

We parallelized our algorithms of raytracing minimization and PLP solver using TBB and OpenMP separately. Thanks to TBB and OpenMP the tasks are scheduled automatically, and we do not need to assign them to the threads by hand. One advantage of TBB is that it provides concurrent data structures, and thus the shared data can be accessed and appended concurrently. In OpenMP version we provided a simple implementation of lock-free data structure. In both version we need to pay attention to the shared data stored in the concurrent vector (for TBB version) or the lock-free data structure (for OpenMP version) when we access the items by their indices, since appending items and increasing the size of the container is not atomic.

The parallel PLP solver is different from the sequential algorithm in some details. To guarantee a unique basis for each region, we found out the "smallest" possible basis for avoiding two threads obtaining different bases with the task points inside the same region, and thus the overlapping regions were avoided.

Chapter 6 Experiments

In this section we will present some experiments on raytracing minimization and the PLP solver for computing the projection and convex hull. We first analyze the advantages and disadvantages of our sequential algorithms, and then we analyze the speedup (the ratio of the runtime of the sequential algorithm for solving a problem to the time taken by the parallel algorithm to solve the same problem) of our parallelism by running it on various of threads.

We used three libraries in our implementation:

-Eigen 3.3.2 for floating-point vector and matrix operations; -FLINT 2.5.2 for rational arithmetic, vector and matrix operations; -GLPK 4.6.4 for solving linear programs in floating-point.

The experiments of minimization and some experiments on the projection and convex hull are performed using our randomly generated benchmarks, as we need a single varying parameter for analyzing the pros and cons. These benchmarks contain randomly generated polyhedra in which the coefficients of constraints range from -50 to 50. In each experiment, we use 10 polyhedra generated with the same parameters. To smooth out experimental noise, we perform the experiments on each polyhedron for 5 times resulting in 50 executions for each set of parameters. Then we sum up the running time of 50 executions, and divide the total time by 5 for obtaining the average execution time of 10 polyhedra. The second part of experiments on projection (and convex hull) used SV-COMP benchmarks [START_REF] Beyer | Automatic verification of c and java programs: Sv-comp 2019[END_REF], by which we will see the performance of our algorithm on real-world programs. We compare the performance on these benchmarks with ELINA library [START_REF] Singh | Fast polyhedra abstract domain[END_REF].

The experiments with sequential algorithms are carried out on a machine with Intel(R) Core(TM) i5-6200U CPU. And the parallel experiments are performed on two servers: S

Future Work

In this chapter we talk about the future work. We will first explain some possible improvements and supplements to the current version of the PLP solver. Then we will show that it is possible to implement a pure floating point PLP solver for computing projection and convex hull, if the rational result is not required. Without rational procedures the soundness of the result can be guaranteed, i.e., the resulting polyhedron is over-approximate, but the precision could be lower. For other operators, the pure floating point solver cannot obtain a sound result, and we need new algorithms to compute these operators using pure floating point arithmetic.

Improving the current algorithm of the PLP solver

In this section we explain the possible improvements on our PLP solver. First the method of checking adjacency could be more efficient. Second the PLP solver could be adapted for solving more general PLP problems which do not have the normalization constraint.

The approach for checking adjacency

We currently store the regions that have been explored into an unstructured array; checking whether an optimization direction is covered by an existing region is done by a linear search. This could be improved in two ways: i) regions corresponding to the same optimum (primal degeneracy) could be merged into a single region; ii) regions could be stored in a structure allowing fast search. For instance, we could use a binary tree where each node is labeled with a hyperplane, and each path from the root corresponds to a conjunction of half-spaces; then each region is stored only in the paths such that the associated half-spaces intersects the region.

Dual degeneracy in general PLP problems

For more general cases, i.e., PLP problems without the normalization constraint, we need to deal with the dual degeneracy. We have two strategies: i) applying perturbation method to the objective function; Verified polyhedra library, VPL, 8 Index Abstract VPL (Verified Polyhedra Library) is an abstract polyhedra domain using constraint-only description. All main operators boiled down to polyhedral projection, which can be computed using Fourier-Motzkin elimination [START_REF] Maréchal | Scalable minimizing-operators on polyhedra via parametric linear programming[END_REF]. This method generates many redundant constraints which must be eliminated at a high cost. A novel algorithm was implemented in VPL for computing the polyhedral projection using parametric linear programming (PLP) [START_REF] Gal | Postoptimal Analyses, Parametric Programming, and Related Topics: degeneracy, multicriteria decision making, redundancy[END_REF], which can replace Fourier-Motzkin elimination. This PLP solver can also be used for computing the join operator (convex hull). Our work focuses on improving the efficiency of the algorithm of PLP solver.

In prior work, PLP [START_REF] Jones | Lexicographic perturbation for multiparametric linear programming with applications to control[END_REF][START_REF] Maréchal | New Algorithmics for Polyhedral Calculus via Parametric Linear Programming[END_REF] was done in arbitrary precision rational arithmetic. In this work, we present an approach where most of the computations are performed in floating-point arithmetic, and then exact rational results are reconstructed. The result obtained from our approach is guaranteed to be sound. We also propose a workaround for a difficulty called degeneracy, which plagued previous attempts at using PLP for computations on polyhedra: in general the (parametric) linear programming problems are degenerated, resulting in redundant computations and duplicates in geometric descriptions. The algorithm of the PLP solver is intrinsically parallelable, however it was developed in VPL with OCaml, which does not well support parallelism programming. In our approach, which is implemented with C++, we proposed a task-based scheme for parallelizing it. Two parallel implementations have been realized using two different parallel mechanisms: Intel's Thread Building Blocks (TBB) 1 and OpenMP tasks [START_REF]OpenMP Application Programming Interface[END_REF].

Résumé

VPL (Verified Polyhedra Library) est un domaine de polyhèdres abstraits utilisant une description uniquement par contrainte. Tous les opérateurs principaux se résumaient à une projection polyédrique pouvant être calculée à l'aide de l'élimination de Fourier-Motzkin. Cette méthode génère de nombreuses contraintes redondantes qui doivent être éliminées à un coût élevé. Un nouvel algorithme a été mis en oeuvre dans VPL pour calculer la projection polyédrique à l'aide de la programmation linéaire paramétrique (PLP), qui peut remplacer l'élimination de Fourier-Motzkin. Ce solveur PLP peut également être utilisé pour calculer l'opérateur de jointure (coque convexe). Notre travail est axé sur l'amélioration de l'efficacité de l'algorithme du solveur PLP.

Dans des travaux antérieurs, le PLP était effectué en arithmétique rationnelle de précision arbitraire. Dans ce travail, nous présentons une approche dans laquelle la plupart des calculs sont effectués en arithmétique en virgule flottante, puis les résultats rationnels exacts sont reconstruits. Le résultat obtenu grâce à notre approche est assuré d'être solide. Nous proposons également une solution de contournement à une difficulté appelée dégénérescence, qui entachait les tentatives précédentes d'utilisation de PLP pour les calculs sur polyèdres : en général, les problèmes de programmation linéaire (paramétrique) sont