Maximum Consensus Floating Point Solutions for Infeasible Low-Dimensional Linear Programs with Convex Hull as the Intermediate Representation

This paper proposes a novel method to efficiently solve infeasible low-dimensional linear programs (LDLPs) with billions of constraints and a small number of unknown variables, where all the constraints cannot be satisfied simultaneously. We focus on infeasible linear programs generated in the RLibm project for creating correctly rounded math libraries. Specifically, we are interested in generating a floating point solution that satisfies the maximum number of constraints. None of the existing methods can solve such large linear programs while producing floating point solutions. We observe that the convex hull can serve as an intermediate representation (IR) for solving infeasible LDLPs using the geometric duality between linear programs and convex hulls. Specifically, some of the constraints that correspond to points on the convex hull are precisely those constraints that make the linear program infeasible. Our key idea is to split the entire set of constraints into two subsets using the convex hull IR: (a) a set X of feasible constraints and (b) a superset V of infeasible constraints. Using the special structure of the RLibm constraints and the presence of a method to check whether a system is feasible or not, we identify a superset of infeasible constraints by computing the convex hull in 2-dimensions. Subsequently, we identify the key constraints (i.e., basis constraints) in the set of feasible constraints X and use them to create a new linear program whose solution identifies the maximum set of constraints satisfiable in V while satisfying all the constraints in X. This new solver enabled us to improve the performance of the resulting RLibm polynomials while solving the corresponding linear programs significantly faster.


INTRODUCTION
Linear programs (LPs) arise in many domains, such as robotics [13,47], databases [10], machine learning [41,48,53,61], computer graphics [32], etc.They are also widely used in the programming languages community for analyzing vulnerabilities in C source code [26], designing correctly rounded math libraries [2,3,37,38,40], repair of deep neural networks [52,55], Presburger arithmetic for polyhedral compilation [46], and many other problems.A linear program is called feasible if all the constraints can be satis ed simultaneously, otherwise, it is called infeasible.Given the widespread use of linear programs, numerous seminal algorithms and methods have satis es all constraints in the rst set (i.e., X) and the slack variables allow some of the additional constraints in the second set to be satis ed.In contrast to our approach, the prior maximum consensus formulation [34] cannot be directly applied to the entire set of constraints because it destroys the low-dimensional nature of the linear program.The resulting linear program that has billions of constraints without the low-dimensional nature cannot be solved by any existing solver.
Convex hull as the intermediate representation.To split the set of constraints into the two sets described above, our idea is to use the convex hull as the intermediate representation (IR).Speci cally, we exploit the geometric duality [14,30] between linear programs and convex hulls and observe that some of the constraints that correspond to points on the convex hull are precisely those constraints that make the low-dimensional linear program infeasible.We compute the convex hull V of the entire system of constraints.Then, we split A into V and X = A − V. We use the RLibm's adaptation of the Clarkson's method to check if X is a feasible LDLP.If so, we have split the entire set of constraints A into a feasible LDLP X and a superset of violated constraints V. Otherwise, we iteratively compute the convex hull of the set X and the new set of violated constraints V ′ .If the set X − V ′ is a feasible LDLP, then we have found our partition and V ∪ V ′ is the superset of the infeasible constraints.The above procedure continues until we end up with a feasible LDLP.This procedure terminates as A is nite.
Exactly computing the convex hull of points that are embedded in the -dimensional Euclidean space R has an asymptotic run-time complexity of ( ) [22], which makes it intractable when is in the order of billions.To address this issue, we observe that the RLibm constraints have a special structure and utilize it to make the problem tractable.The RLibm linear constraint for an input is of the form, ≤ 0 + 1 + 2 2 + 3 3 + .. + ≤ ℎ, where 's are unknowns and /ℎ represents lower/upper bounds of the rounding interval (see Section 2.1).Using geometric duality transformations (see Section 2.3), these constraints represent the points, ( , 2 , 3 , ..., , ) or ( , 2 , 3 , ..., , ℎ), whose convex hull is being computed in our approach.We observe that these points are intrinsically 2-dimensional because the points ( , 2 , 3 , ..., ) are parameterized by a single parameter and represent a 1-dimensional curve in high dimensions.Further, (or ℎ) indicates the freedom to satisfy the constraint.Hence, we propose to project points in R to a 2-dimensional space (i.e., R 2 ) and compute the convex hull in R 2 .The complexity of computing the convex hull of points in R 2 with our projections is ( log ), where is the number of points on the convex hull [33].This makes our approach extremely fast in practice and can still handle billions of constraints.Theoretically, we have to compute the convex hull with 2-D projections iteratively to identify the superset of infeasible constraints in A. Empirically, we found that computing the convex hull with a 2-D projection identi ed a small superset of the infeasible constraints in A in a single iteration.
Our prototype can solve infeasible linear programs with billions of constraints signi cantly faster while producing the best solutions.Using the solutions from our solver, we have also improved the performance of correctly rounded RLibm's math libraries.

BACKGROUND ON LINEAR PROGRAMS IN THE RLIBM PROJECT
We provide background on our RLibm project, the linear programs generated with the RLibm approach, and the notion of geometric duality between convex hulls and linear programs that will be necessary to understand our solution for infeasible linear programs.

The RLibm Project
The RLibm project proposes a new method to build correctly rounded math libraries by generating polynomial approximations that approximate the correctly rounded result rather than the real value  of an elementary function [2,3,37,38,40].Typically, the implementation of a correctly rounded function for a 32-bit oat representation uses the 64-bit double precision representation internally.Given a correctly rounded result for a 32-bit oat input, there is an interval of values in the double precision representation around the correctly rounded result such that any value in that interval rounds to the correctly rounded value (see Figure 1(A)), which is called the rounding interval.The rounding interval for an input is represented as [ , ℎ ], where is the lower bound and ℎ is the upper bound.When the goal is to generate a polynomial of degree − 1 with terms, the rounding interval speci es the linear constraint shown in Figure 1(B) on the result of the polynomial for input .The RLibm project generates a system of linear inequalities corresponding to all 32-bit inputs and their corresponding rounding intervals.Then, the goal is to identify the coe cients (i.e., 's) of a polynomial of a particular degree that satis es these inequalities.The RLibm project uses an LP solver to solve this system of inequalities.
Before generating such polynomial approximations, it is necessary to reduce the original input from the domain of a 32-bit oat representation (i.e., [2 −150 , 2 128 )] to a small domain (e.g., [0, 1]).This step is called range reduction.The original input is range reduced to ′ .Subsequently, the polynomial that we wish to generate approximates the result for ′ , which is then output compensated to compute the nal output for .The RLibm project uses speci c range reduction methods that ensure that the reduced inputs are positive, which we leverage in our approach.
The RLibm project has also proposed a method to generate a single polynomial approximation that can produce correctly rounded results for all FP representations up to -bits [40].The key idea is to generate a polynomial that produces correctly rounded results for the ( + 2)-bit representation (which has 2 additional precision bits compared to the -bit representation) using the round-to-odd rounding mode.When that result is double-rounded to any representation with less than -bits, it produces correct results for the target representation.The system of linear constraints generated to produce correctly rounded results for the ( + 2)-bit representation with the round-to-odd mode for all inputs also makes them infeasible in many cases.

Solving Low-Dimensional Linear Programs
The trajectory of our RLibm project has been shaped by the ability to solve linear programs with billions of constraints.Initially the RLibm project used piecewise polynomials because there were no publicly available solvers that could handle such a large number of constraints.Another challenging issue is to generate oating-point (FP) solutions because a real-valued solution from an LP solver when rounded to the FP representation may not satisfy all the constraints.
Given a reduced input and its rounding interval [ , ℎ ], the linear constraint generated to create a polynomial of degree − 1 with terms is shown in Figure 1(B).We will canonicalize the above constraint in the standard LP format that has just one inequality in each constraint.Putting all the unknown polynomial coe cients in a column vector , the constraint in Figure 1(B) can be re-written in the following equivalent form: where the row vector = [1, , 2 . . ., −1 ] and +1 = [−1, − 1 , − 2 . . ., − −1 ].There are two inequalities for each RLibm constraint after canonicalization.The number of unknown variables (i.e., 's) is also known as the dimension of the linear system.
Low-dimensional linear program.When the dimension of is much smaller than the number of constraints, then the linear program is called a low-dimensional linear program (LDLP).There is a large amount of redundancy in LDLPs, in the sense that a few key constraints determine the solution for all the other constraints [19].
Solving feasible low-dimensional linear programs.The seminal results from Meggido [42] and Clarkson [20] describe algorithms to solve large feasible LDLPs.If the LDLP is feasible, then the small dimension property can be used to design a fast randomized algorithm [20] for computing the solution with real values.The RLibm project uses Clarkson's method [20] to solve feasible LDLPs with oating-point solutions [2].The essence of this algorithm is to employ weighted random sampling [23] to choose a small subset of 6 2 constraints S from the entire set of constraints A. It associates a weight with each constraint.Initially the weight is 1 for all constraints.The algorithm is as follows: • Step 1. Sample S constraints from A with weighted random sampling where |S| = 6 2 , where is the number of terms in the polynomial.• Step 2: Solve the sample S optimally using an LP solver to compute the sample solution x ★ .
• Step 3: Exit if x ★ satis es all constraints in A, otherwise check how many constraints of A are not satis ed by x ★ .If the sum of the weights of the violated constraints is less than or equal to 1/(3 − 1) fraction of the combined sum of the weights of constraints satis ed by x ★ , then discard this sample and go to Step 1. Otherwise, double the weights of all violated constraints and then go to Step 1.
These three steps are repeated until a solution x ★ satis es all constraints in A. Clarkson [20] provides seminal results on the convergence of this procedure.Speci cally, with probability at least 1/2, x ★ can only violate 1/3 constraints in A. In other words, the above procedure converges to a solution that satis es all the constraints in A in 6 log iterations in expectation.
To solve feasible LDLPs with the above algorithm, the RLibm project uses the SOPLEX solver to generate the solution for the sample.SOPLEX solves the sample with real values (i.e., with GMP rational numbers).When the sample solution is rounded to oating point coe cients, some constraints in the sample may not be satis ed.In such cases, the RLibm project restricts the intervals to check if the sample can still be solved.The need to generate a oating point solution can also turn the feasible LDLP into an infeasible LDLP.
What happens with infeasible LDLPs?When the LDLP is infeasible, Clarkson's algorithm [20] adapted by the RLibm project is not guaranteed to converge to a solution.Running more iterations of the above algorithm with an infeasible LDLP does not help in reducing the number of constraints violated by the sample solution.The above algorithm uses weighted random sampling and doubles the weight for the violated constraints.Although this doubling of weights for the violated constraints increases the probability of them being selected in the next sample, the entire sample becomes infeasible.In many LDLPs, it fails to generate an acceptable solution that satis es the maximum number of constraints.In fact, the state-of-the-art approach is to run this algorithm for a su ciently long time and choose the best result that has been generated so far.

Geometric Duality and Linear Programming
Our proposed method to solve infeasible LDLPs uses ideas from computational geometry, speci cally the geometric duality transformations between points and lines.We provide a brief background on it.We refer the reader to excellent books on this subject for detailed information [11,30].
Geometric duality transformations.For simplicity of exposition, consider the 2-D case rst.A geometric duality : R 2 → R 2 is a transformation that maps a set of lines to a set of points and vice versa, in a one-to-one manner.Let = ( , ) be a 2-D point, then the dual of is denoted by ★ and is de ned as the line = − .The dual of a line = + is de ned as the point ( , − ).Duality transformations have several nice properties that can be useful for analyzing geometric structures.For example, duality transforms are incidence preserving, i.e., a point lies on a line if and only if the dual point ★ lies on the dual line For example, three lines 1 , 2 and 3 are concurrent in the primal space if their dual points ★ 1 , ★ 2 and ★ 3 are collinear in the dual space.More generally, three lines 1 , 2 and 3 are approximately concurrent in the primal space if their corresponding dual points ★ 1 , ★ 2 and ★ 3 are approximately collinear in the dual space, as shown in Figure 2. The latter problem can be e ciently solved using linear regression.A popular result in computational geometry, which our proposed approach builds upon, is that the lower envelope of an arrangement of lines (i.e., the lowest partition of R 2 that is produced by the arrangement) corresponds to the upper convex hull of the dual points [12], as shown in Figure 4(B).In > 2 dimensions, geometric duality : R → R maps a hyperplane ℎ to a point ℎ ★ and vice versa, as de ned below: Visualizing Linear Programs through the lens of Geometric Duality.For visualizing constraints in an LDLP from a geometric viewpoint, consider the constraint below: Fig. 3.The shaded region denotes the solution space for an LDLP that is bounded from above by upper halfspaces (red) and from below by lower halfspaces (blue).
Equation (2) represents a halfspace in R , i.e., a partition of R into two pieces by a hyperplane that only includes one piece.If the coefcients 1 , . . ., are all positive, then equation (2) describes an upper halfspace, i.e., it is bounded from above, as shown in Figure 4(A)(b).Similarly, if the coe cients 1 , . . ., are all negative, then equation (2) describes a lower halfspace, i.e., it is bounded from below, as shown in Figure 4(A)(a).Thus, the solution space for a system of linear constraints, which is the intersection of all these halfspaces, can   Fig. 5. Our approach to generate an FP solution for an infeasible LDLP that satisfies the maximum number of constraints.The first step is to split the entire set of constraints (A) into two sets X 1 and V 1 where X 1 is a feasible LDLP.Using the convex hull as the IR, we compute the superset of the infeasible constraints,V 1 , using geometric duality.Second, we use the basis of the solution for X 1 (i.e., B 1 ) and a new linear program with slack variables that satisfies the maximum number of constraints in V 1 to get a new feasible LDLP X and a small set V. Third, we perform an iterative algorithm with weighted random sampling from X with the maximum consensus formulation for constraints in V to produce an FP solution that satisfies the maximum number of constraints in A.
be viewed as being bounded from above by upper halfspaces, and being bounded from below by lower halfspaces, as shown in Figure 3. Constraints that border the solution space are precisely the set of key constraints that determine the solution for all the other constraints.This visualization also presents an opportunity to e ciently compute the set of constraints that border the solution space.Speci cally, if we only consider the set of upper halfspaces, then their lower envelope corresponds to the set of constraints that border the solution space (see Figure 3).This lower envelope can be e ciently computed using the upper convex hull of the dual points [12], as shown in Figure 4(B).In a similar fashion, the set of lower halfspaces that border the solution space can be e ciently identi ed by computing the upper envelope of the lower halfspaces, which corresponds to the lower convex hull of the dual points.
Identifying infeasible constraints in an LDLP.For an infeasible LDLP, the shaded region shown in Figure 3 does not exist, or in other words, the intersection of the lower envelope of the upper halfspaces and the upper envelope of the lower halfspaces is empty.In this case, identifying the two envelopes can still provide crucial information regarding the set of constraints that make the LDLP infeasible.Speci cally, constraints that border the two envelopes are the most con icting constraints, and removing them is the rst step towards making the LDLP feasible.If the LDLP is still infeasible, then this process can be repeated until the LDLP becomes feasible (see Figure 6).

AN OVERVIEW OF OUR APPROACH
Given an infeasible low-dimensional linear program with a large number of constraints (i.e., billions of them) but with a small number of unknown variables, our goal is to nd a oating point solution that satis es the maximum number of constraints.Our approach is motivated by such infeasible LDLPs generated in the RLibm project to generate polynomial approximations with oating point coe cients to produce correctly rounded math libraries.This section provides an overview of our approach before going into the details in the subsequent sections (as illustrated in Figure 5).
Existing methods for infeasible LDLPs.Since we have an infeasible linear program, there does not exist a single solution that satis es all the constraints.Our work is inspired by the maximum consensus formulations in the computer vision community where they create another feasible linear program by using slack variables whose solution satis es the maximum number of constraints in the original linear program [34].This approach creates four additional constraints for each original linear constraint and destroys the low-dimensional nature of the linear program.Hence, it works only with a small number of constraints typically in the order of thousands of constraints.When we tried that formulation with our infeasible LDLPs, it could not produce a solution even with real values.The huge number of constraints requires us to explore linear or near-linear time algorithms (i.e., even an algorithm with a time complexity of ( 2 ) makes it impractical for our purposes) while preserving the low-dimensional nature of the linear program.
Insights.Our high-level idea to solve infeasible LDLPs with a large number of constraints is to partition the set of constraints into two sets: a set of constraints that is feasible and a small superset of infeasible constraints that makes the original set infeasible.This partitioning enables us to intelligently combine existing solvers for feasible LDLPs from the RLibm project [2] and the maximum consensus formulation from the vision community [34] to produce maximum consensus solutions for infeasible LDLPs with approximately billions of constraints.Figure 5 illustrates our approach to solve a large infeasible LDLP in the RLibm context.The research challenges that we have to address are the following: (1) how do we perform the above partitioning quickly?(2) how do we combine the maximum consensus formulation that handles small number of constraints with the solver for feasible LDLPs such that we can generate maximum consensus solutions for a system with billions of constraints?
3.1 Convex Hull as the IR to Create a Feasible LDLP One of the key ideas of this paper is to use the convex hull as the intermediate representation (IR) to perform the partitioning of the constraints.Geometric duality transforms, which we described in Section 2, allow us to represent a linear constraint as a point in the dual space.Then, nding the feasible region (i.e., a solution) for all the constraints can be done by identifying the convex hull in the dual space.Thus, we use the convex hull as the IR to solve large infeasible LDLPs.It also enables us to perform further optimizations and approximations that would not have been feasible otherwise, similar to how compiler IRs enable canonicalization, optimization, and approximation while compiling high-level languages.
We represent every linear constraint as a point in the dual space.Subsequently, we identify the upper hull of the points in the dual space corresponding to the constraints that specify the lower bound and identify the lower hull of the points in the dual space corresponding to the constraints that specify the upper bound.Our insight is that the constraints that make the LDLP infeasible must lie either on the upper or the lower hull in the dual space (see Figure 6).
Given the original set of constraints A, we compute the lower hull L and the upper hull U, and their union i.e., V 1 = L ∪ U. Subsequently, we compute the set X 1 = A − V 1 .We check if the set X 1 is a feasible LDLP using the Clarkson method adapted to the RLibm context [2].If so, we Fig. 6.Intuition on why computing the lower hull and upper hull iteratively enables the identification of infeasible constraints.Each line represents a constraint.The bordered lines represent constraints on the convex hull of the dual points.A er we iteratively remove the constraints on the convex hull, we end up with a feasible LDLP.
have successfully partitioned A into a feasible LDLP X 1 and a superset of infeasible constraints V 1 .Otherwise, we compute the lower hull and upper hull of X 1 and check if removing them from X 1 makes it feasible.We repeat this procedure iteratively until we end up with a feasible LDLP.The union of all the lower hull and upper hull constraints identi ed in the intermediate steps constitutes the superset of the infeasible constraints.Given that the set A is nite, the process terminates either when we identify a feasible LDLP or when the size of the superset of infeasible constraints identi ed exceeds some threshold (e.g., 20,000).Figure 6 pictorially shows how removing the lower and upper hull constraints identi es the infeasible constraints and eventually results in a feasible LDLP.
Approximation using the convex hull IR.Computing the convex hull with constraints and in dimensions has time complexity ( ), which makes the problem intractable when we have billions of constraints in the RLibm context.We observe that RLibm constraints have a special structure; they are intrinsically 2-dimensional because a given constraint depends on a single input (which is positive) and the bounds or ℎ (i.e., 0 + 1 + 2 2 + .. + ≤ ℎ ).Hence, we propose to project points from -dimensions to 2-dimensions and compute the convex hull in 2-dimensions.We subsequently identify the feasible LDLP using the computed 2-dimensional hulls with an iterative process that we described above.We empirically show that computing the convex hull in 2-dimensions in a single iteration is su cient to identify a superset of infeasible constraints in the RLibm context.Our use of the convex hull as the IR and approximation over it by projecting to 2-dimensions, makes our method extremely fast in practice with time complexity ( log ), where is the number of vertices on the convex hull.

Maximum Consensus using the Basis of the Feasible LDLP
After computing the feasible LDLP X 1 and a superset of the infeasible constraints V 1 using the convex hull, we now want to identify the maximum number of constraints in V 1 that can be satis ed while also satisfying all the constraints in X 1 .Our high level strategy is to create a new linear program with slack variables that satis es all the constraints in X 1 and the maximum number of constraints in V 1 inspired by prior maximum consensus formulations [34].However, we cannot directly use the maximum consensus formulation because we have more than a billion constraints.
Using the maximum consensus formulation for all the billion constraints will violate the "small dimension" property (small number of variables and a large number of constraints).General linear programs with billions of constraints cannot be solved by any modern LP solver.
To address this issue, we leverage the insight that any feasible LDLP has a small set of "key" basis constraints and any solution (with real values) that satis es these key constraints will satisfy all other constraints in the LDLP [19].We use the feasible LDLP solver from the RLibm project to identify (an overapproximation of) the basis constraints B 1 , which has 6 2 constraints where is the number of unknowns in the LDLP (i.e., dimension of the LDLP).
Next, we generate a new LP formulation that uses slack variables for the constraints in V 1 and the original constraints for the basis constraints of B 1 .The use of slack variables makes the new LP feasible.Because the number of variables in the new LP increases the number of unknowns, the new linear program is not a low-dimensional linear program.Given that the superset of the infeasible constraints computed using the convex hull has already reduced V 1 to a few thousand constraints and we use the basis of X 1 that just has 6 2 constraints, the new feasible LP in total has a few thousand constraints, which can be easily solved by modern LP solvers.Section 5 provides detailed explanation to create the maximum consensus solution using B 1 and V 1 .
When we solve the above new linear program using any existing LP solver (e.g., SOPLEX), it produces solutions with real values.Our goal is to produce oating point maximum consensus solutions.When we round the real-valued solution to oating point values, it may violate some of the constraints in B 1 or V 1 .In those cases, we attempt to solve the system by strengthening the constraint to account for rounding errors (i.e., decreasing the upper bound or increasing the lower bound).It forces the LP solver to nd a solution for a much stricter constraint.When we are able to nd a solution, we have identi ed a nearly maximum consensus solution.The middle part of Figure 5 illustrates these steps.

Combining Clarkson's Method with the Maximum Consensus Approach
The solution obtained from the above step is the nearly-maximum consensus solution (not the maximum consensus solution) due to FP rounding errors.The challenging issue in generating FP solutions arise when the solution from the above step satis es B 1 and the maximum number of constraints in V 1 but violates another constraint in X 1 .To identify the maximum consensus solution, we split the entire set of constraints into two sets: X consists of constraints that are satis ed by the nearly maximum consensus solution and V consists of constraints violated by the nearly maximum consensus solution.
We iteratively combine Clarkson's method for feasible LDLPs in the RLibm project [2] (i.e., for X) and the maximum consensus formulation for violated constraints (i.e., V).The goal is to identify an FP solution that satis es all the constraints in X and the maximum number of constraints in V. We could have directly used this iterative method combining Clarkson's method and the maximum consensus formulation on the feasible set X 1 and the violated set V 1 identi ed using the convex hull.However, each iteration of the Clarkson's method would have been signi cantly slower because there are thousands of violated constraints in V 1 .The nearly maximum consensus solution ideally reduces the cardinality of violated constraints from thousands to a handful.
We assign weights to each constraint in X.We sample 6 2 constraints from X and use all the constraints from V with the maximum consensus formulation and ask an o -the-shelf LP solver for a solution.The use of slack variables for the constraints in V enables us to create a feasible linear program for the sample.We validate if the sample solution satis es all the constraints in X.If they are violated by the sample solution, then the weights of those constraints in X are doubled.When all the constraints in X are satis ed, we have an FP solution that satis es the maximum number of constraints.The rightmost part of Figure 5 shows the iterative loop.

FINDING THE SUPERSET OF INFEASIBLE CONSTRAINTS WITH THE CONVEX HULL
We make a case for using the convex hull as the intermediate representation to generate maximum consensus solutions for infeasible low-dimensional linear programs.We achieve this by using geometric duality transformations to represent a constraint as a point in the dual space [30].Looking at the convex hull of the points in the dual space can help in the identi cation of feasible regions in the case of feasible LDLPs [51].As explained in Section 2, we observe that by visualizing constraints as halfspaces, the lower envelope of all upper halfspaces (resp.upper envelope of all lower halfspaces) includes all constraints that make the LDLP infeasible.The upper and lower envelopes can be e ciently computed by using geometric duality and computing the lower and upper hull of the dual points.Hence, the convex hull serves as a good intermediate representation for computing maximum consensus solutions for infeasible linear programs.
Linear constraints to points in the dual space.The rst step in our approach is to represent each constraint as a point in the dual space.Given a linear constraint from the RLibm project of the form ≤ 1 + 2 1 + 3 2 + 4 3 ≤ ℎ , we rst canonicalize it as two constraints: The points in the dual space for these constraints are (1, 1 , 2 , 3 , ℎ ) and (−1, − 1 , − 2 , − 3 , − ), respectively.Since the rst coordinate is the same for all the dual points (1 for dual points of constraints that specify the upper bound and −1 for dual points of constraints that specify the lower bound), we can ignore the rst coordinate, making the mapping from R to R .
Special structure of the RLibm constraints and computing the convex hull with 2-D projections.Computing the convex hull in -dimensions has complexity ( ), which becomes intractable when is in order of billions.To address this issue, we leverage the special structure of the RLibm constraints and compute the convex hull in 2-dimensions.We project the point set P that is embedded in a high-dimensional Euclidean space R to a lower-dimensional Euclidean space R 2 .The RLibm constraints are of the form, for example, ( , 2 , 3 , ..., , ) or ( , 2 , 3 , ..., , ℎ) where and ℎ are real numbers.These points are intrinsically 2-dimensional because the points ( , 2 , 3 , ..., ) are parameterized by a single parameter and represent a 1-dimensional curve in high dimensions.
As we show next in Theorem 1, the points ( , 2 , 3 , ..., ) also lie on a convex manifold.Moreover, they are monotonically increasing.Hence, if we ignore the last coordinate ( , ℎ), then all the constraints are satis able.Intuitively, and ℎ represent the freedom provided by the RLibm approach for computing the correctly rounded result.If this freedom is large, then it makes it easier for a single polynomial approximation to exist.Thus, it is only because of the last coordinate (or ℎ) that some constraints become infeasible.Our approach of computing the 2-dimensional hull is based on the intrinsic 2-dimensionality of the point set described above.Technically, the constraint with the highest coordinate in any direction should be below the constraint with the lowest ℎ coordinate in that direction for the LDLP to become feasible.This observation is our main motivation for using the 2-D convex hull for nding infeasible constraints (i.e., with the extreme and ℎ coordinates).
We now explain the details behind this projection by rst proving the following theorem: Theorem 1.All points of the form ( , ), where 1 ≤ < , lie on a convex curve.
Proof.The points ( , ) lie on the graph of the function = / .The curvature of this function can be computed as: In the context of the RLibm project, is the reduced input, which is always positive [38,39].Furthermore, since > ≥ 1, it follows that the curvature is always positive.The positivity of the curvature implies that the function = / is convex.□ The convex hull of a point set P is the smallest convex set containing P. However, if P was sampled from a convex curve, then by de nition, all points in P lie on the convex hull.An argument similar to that made in Theorem 1 can also be made for all triplets, quadruplets, etc., of powers of , implying that they lie on a convex manifold.However, we skip this proof for simplicity of exposition as it follows directly from the structure of cyclic polytopes [1].
The above argument leads us to individually consider the pairs ( , ℎ ), ( 2 , ℎ ), . . ., ( , ℎ ) when computing the convex hull for the dual points corresponding to constraints that specify the upper bound.In a similar fashion, we only consider the pairs (− , ), (− 2 , − ), . . ., (− , − ) when computing the convex hull for the dual points corresponding to constraints that specify the lower bound.Intuitively, this means that it is the lower and upper bounds and ℎ that ultimately decide whether the linear program is feasible or infeasible.Indeed, as the gap between and ℎ increases for all inputs, so does the likelihood of a possible solution for the linear program.
We provide a proof for the general case (i.e., without the special structure of the RLibm constraints) that points identi ed using the convex hull of the 2-D projection maps to the convex hull of the dual points ( , 2 , . . ., , ℎ ) or ( , 2 , . . ., , ) in -dimensions.For this purpose, we recall that a point lies on the convex hull, if and only if there exists a hyperplane that passes through and all the other points lie only on one side of , i.e., in one halfspace de ned by [11].
Theorem 2. Consider a point set P in -dimensions and let Q be a 2-D point set that is computed by projecting all points in P to a 2-D space by choosing two of the coordinates.Then, any point on the convex hull of Q maps to a point on the convex hull of P.
Proof.Without loss of generality, assume that the point set Q was computed by projecting all the -dimensional points in P along the rst two coordinates.Let ≡ ( 0 , 1 ) be a point on the convex hull of .By de nition of the convex hull, there is a hyperplane that passes through and all other points in Q lie only on one side of .Let ≡ ( 0 , 1 ) be the normal vector for .The above de nition implies that 0 • 0 + 1 • 1 ≤ 0 • 0 + 1 • 1 for all other points ≡ ( 0 , 1 ) in Q.If we consider the normal vector ≡ ( 0 , 1 , 0, . . ., 0) in -dimensions, then it passes through the point ′ ≡ ( 0 , 1 , 2 , . . ., −1 ) which was projected to .Moreover, for all points ′ ≡ ( 0 , 1 , 2 , . . ., −1 ) in P, we still have the inequality 0 from the de nition of the 2-D convex hull, i.e., the point ′ lies on the convex hull of P. □ In summary, without the structure of the RLibm constraints, points identi ed via the 2-D hull lie on the -D hull, but not vice versa.The speci c structure of the RLibm constraints (i.e., they are inherently 2-dimensional) enables us to identify the infeasible constraints using the 2-D hull.
Identifying a superset of infeasible constraints by computing the convex hull with 2-D projections.To identify the superset of infeasible constraints, our approach is to use an iterative procedure, which we described in Section 3.1.We compute the convex hull with the 2-D projection and check if the remaining constraints after removing the constraints on the 2-D hull are feasible using the RLibm's feasible LDLP solver.This process repeats until we end up with a feasible LDLP or the total number of constraints accumulated by computing the 2-D hull iteratively exceeds the user-de ned threshold.Crucially, our approach always reports a superset of infeasible constraints by construction because we have a de nitive way to check whether a system is a feasible LDLP using RLibm's solver.We empirically show in Section 7 that the convex hull computed with the 2-D projection computes a superset of infeasible constraints in a single iteration in the RLibm context.
Computing the convex hull in 2-dimensions with our projections can be done in ( log ) time using the Kirkpatrick-Seidel algorithm [33] where is the number of vertices on the convex hull.This computation is performed with high precision arithmetic using the CGAL library [56].With this approach, we are able to split the set of constraints A into a feasible set of constraints X 1 and a superset of infeasible constraints V 1 .

COMPUTING NEARLY MAXIMUM CONSENSUS SOLUTIONS
Using the convex hull of the constraints in the dual space, we have partitioned the entire set of constraints (i.e., A) into a set of feasible constraints (i.e., X 1 ) and a superset of infeasible constraints (i.e., V 1 ).The set V 1 can have a few thousand constraints.Now, the goal is to identify a solution that satis es all the constraints in the rst set and the maximum number of constraints in the second set.
Identifying the basis of the feasible LDLP X 1 .Given that the feasible set is an LDLP, there exists a small set of basis constraints (i.e., B 1 in Figure 5) whose solution satis es all constraints in the feasible LDLP X 1 [19].For the feasible set of constraints, we can use the solver from the RLibm project for feasible low-dimensional linear programs [2] to identify the basis constraints.Inspired by maximum consensus formulations from the computer vision community [34], we design a new linear program whose solution in real values satis es all the constraints in B 1 and the maximum number of constraints in the set V 1 .We describe our adaptation of the maximum consensus formulation in Section 5.1.The new linear program for determining maximum consensus introduces two new variables for each original constraint.Hence, it cannot be applied to all the constraints in A as modern LP solvers cannot solve regular LP problems (i.e., number of variables are proportional to the number of constraints) with more than a few thousand constraints.
Transitioning from real valued solutions to FP solutions.The solution that we obtain by solving the basis constraints and the new LP formulation for the constraints in V 1 is with real values.When we round the real valued solution to a oating point representation, it may fail to satisfy some constraints in B 1 .In such cases, we reduce the bound (i.e., ) of the constraint by 1 ULP (units in the last place) in double precision (i.e., the representation used for the implementation of the math library) for constraints of the form ( − ≤ 0).This is analogous to shrinking of the intervals in the RLibm project.By doing so, we ask the solver to solve a much stricter constraint.We repeat this process until we are able to solve the set B 1 ∪ V 1 .This process terminates because B 1 is a feasible LDLP.When the solution with oating point values satis es all constraints in B 1 , we have a nearly maximum consensus solution.
The solution may still not be the maximum consensus solution for the entire set of constraints A because the rounding error introduced by rounding the real value to a FP solution may violate a few constraints in A. The important point to note is that the number of violated constraints will signi cantly shrink from a few thousands to a handful by identifying the nearly maximum consensus solution.Subsequently, we employ an iterative algorithm that combines maximum consensus with Clarkson's method [20] to nd the maximum consensus solution, which we describe in Section 6.

A New Linear
Program for Maximum Consensus Among V 1 and B 1 Given the set of basis constraints B 1 and an overapproximation of violated constraints V 1 , the set B 1 ∪ V 1 is an infeasible LDLP.We want to identify a solution that satis es the maximum number of constraints in V 1 while satisfying all constraints in B 1 .Mathematically, This formulation, as stated above by itself, is not a linear program.The objective function is counting over discrete sets.Our goal in this section is to create an LP formulation.
Slack variables to move to a feasible LP and indicator variables to count the number of violated constraints.We want to create a new feasible LP from an infeasible LP.There are four key ideas in creating this formulation.First, we add slack variables to each constraint.Each original constraint − ≤ 0 now becomes − ≤ , where ≥ 0 is a slack variable.The slack variable provides a bit more freedom to satisfy the constraint.We want the slack variable to be non-zero only when the constraint is violated, which we accomplish by adding additional constraints that we describe next.When the ℎ constraint is violated in the original LP, we have − > 0. With the slack variables, we have − ≤ .Hence, will be non-zero when the constraint is violated.
Second, we add indicator variables for each constraint to indicate whether the constraint is violated or not.With the use of indicator variables, nding a solution that satis es the maximum number of constraints boils down to minimizing the sum of the indicator variables (see equation ( 5)).
Third, we want the indicator variable to be 1 when the constraint is violated.Since the slack variable is non-zero when the constraint is violated, we force the indicator variable to be 1 with the constraint (1 − ) = 0.
Finally, we want to force the indicator variable to be 0 and the slack variable to be 0 when the ℎ constraint is satis ed.We accomplish this with the constraint ( − + ) = 0.When the original LP constraint is satis ed, − ≤ 0. Hence, − + ≥ 0. Given that ≥ 0, we have − + ≥ 0, which forces to be zero.The constraint (1 − ) = 0 forces to be zero when is 0. The constraint ( − + ) = 0 forces − + = 0 when the ℎ constraint is violated because is 1 in that case.Further, it also ensures that the slack variable captures the amount by which the constraint + ≤ 0 is o in the original LP.Summarizing these ideas, the formulation to identify the maximum set of constraints in V 1 while satisfying all constraints in B 1 is as follows: Note that all constraints from the basis of X 1 (i.e., B 1 ) are used unmodi ed from the original LP (see equation ( 7)).Slack and indicator variables are only introduced for the constraints in V 1 .
Removing non-linearity with iterative predictor and corrector steps.Equations ( 5)-( 11) are still non-linear due to the constraints in equations ( 8) and (9).We use the idea from prior work [34] of removing this non-linearity by moving the non-linear constraints in equations (8) and ( 9) into the objective function.These new constraints are controlled by the penalty parameter , which forces as many constraints to be satis ed as possible. min The above system has a non-linear objective function.The idea is to solve it in an iterative fashion with a predictor-corrector loop.In the predictor step, we x the values of , which makes them constant and the system linear.The resulting LP can be solved by an LP solver to obtain the values of and .Subsequently in the corrector step, we use the values for and from the predictor step and make them constant and create a linear program and solve them for ′ .
The linear program being solved by the predictor step after xing the variables with constants is shown in equation ( 13), which has been simpli ed with constant propagation for the variables and some algebraic simpli cation.
After obtaining the solution to equation ( 13), the corrector step xes the variables and and solves for the following set of constraints in equation ( 14) to determine the 's.
Interestingly, this equation can be solved in closed form in a way that also ensures that ∈ {0, 1} .Speci cally, if 1 − ( − ) ≤ 0, set = 1, else set = 0.In the subsequent iterations, the predictor uses these values of .The parameter, when set to a large value, will force the 's to be 1, unless the constraint is satis ed.The predictor will try to satisfy as many violated constraints as possible in the next iteration.We alternate between the predictor and the corrector steps until the objective function of the predictor step in equation ( 13) does not change across iterations.

Illustration of Maximum Consensus with B 1 and V 1
We illustrate our above formulation to nd the nearly maximum consensus solution with a small infeasible LP with three basis and two violated RLibm constraints.Assume that the polynomial ( ) that we are generating from the RLibm project has the form ( ) = 0 + 1 + 2 2 .Let us say the basis constraints (i.e., B 1 ) corresponding to the inputs 1 , 2 , and 3 are: Let's say the violated constraints, V 1 , corresponding to inputs 4 and 5 are: After canonicalization, each original basis constraint from equation ( 15) becomes, Similarly, the violated constraints from equation ( 16) after canonicalization become, Our goal is to nd a nearly maximum consensus solution that satis es all the basis constraints and the maximum number of the violated constraints.Obviously, in this particular scenario, it is possible to exhaustively test whether the solution is maximal by inserting each of the violated constraints one by one and checking if the linear program is feasible.However, in practical scenarios, this is not a feasible approach when there are several billions of constraints.
We introduce slack and indicator variables for each of the constraints in the set of violated constraints.We use 4 and 4ℎ to represent slack variables for the violated constraint in equation (18).Similarly, 5 and 5ℎ represent slack variables used for the violated constraint in equation (19).Similarly, 4 , 4ℎ , 5 , and 5ℎ represent indicator variables corresponding to the violated constraints in equations ( 18) and (19), respectively.
Initially, we can set indicator variables (i.e., 's) to either 0 or 1 for the predictor step when we try to nd the 's and .The following linear equation from the predictor step (i.e., equation ( 13)) attempts to nd and satisfy the maximum number of violated constraints while satisfying all the basis constraints: Note that basis constraints are used without any slack or indicator variables in equation (20).The introduction of the slack variables 4 , 4ℎ , 5 , 5ℎ ensures that the constraints [ 4 , ℎ 4 ] and [ 5 , ℎ 5 ] may not be exactly satis ed by the computed solution at either the lower or the upper bound, while ensuring that the linear program de ned in equation ( 20) still computes a non-trivial solution, i.e., some of the slack variables may be non-zero.This is in stark contrast to prior solvers for feasible LPs [2], which will simply fail to nd a solution if all the constraints cannot be exactly satis ed.
The corrector stage uses the solution from the LP problem above and computes 4 , 4ℎ , 5 , and 5ℎ using equation (14).We set to a large pre-determined value to provide higher penalty to violated constraints.The process is now repeated with the corrected values until the objective function in equation ( 13) does not change or the user-speci ed number of constraints are satis ed.

ITERATIVELY FINDING THE MAXIMUM CONSENSUS FP SOLUTION
The new linear program for maximum consensus with the basis and an overapproximation of the violated constraints is solved with real values and the solution is rounded to oating point values.The oating point solution satis es all the basis constraints and the maximum number of violated constraints.However, due to rounding errors, the oating point solution may violate some constraints in X 1 , the feasible LP.There is also some bias introduced by the xed basis B X 1 because the basis computation was decoupled from the set of violated constraints V 1 .Hence, the solution from the formulation presented in Section 5 is a solution that is close to the maximum consensus solution but not the maximum consensus solution itself.Using the above solution, we partition the entire set of constraints into two sets: a set of constraints X that is satis ed by the nearly maximum consensus solution and a set of violated constraints V.The set of violated constraints V is signi cantly smaller than the overapproximation of the violated constraints that we started with (i.e., V 1 ).Next, we propose an iterative combination of the maximum consensus formulation with Clarkson's algorithm [20] for feasible LDLPs.Clarkson's method tries to remove the bias introduced by any xed basis by iterating over many di erent bases for X, in an e ort to nd the basis that is best suited for satisfying the maximum number of constraints in V.Each iteration with the sets X and V is signi cantly faster in comparison to running our iterative algorithm with X 1 and V 1 , which were generated by our approximation using the convex hull.
Here are the steps of our iterative algorithm, which uses Clarkson's method for feasible LDLPs from the RLibm project and combines it with the maximum consensus formulation that we described in Section 5.1.Given the feasible LDLP X and a small set of violated constraints V, our goal is to generate a solution with oating point values that satis es the maximum number of constraints in V while satisfying all the constraints in X.
• Step 1: We maintain weights with each constraint from X.We initialize the weights to 1.
• Step 2: Now we sample 6 2 constraints from X using weighted random sampling, where is the number of unknowns.We use the original LP for the sampled constraints.We generate the maximum consensus LP formulation for the constraints in V. We ask the LP solver to solve the sampled constraints from X along with the maximum consensus constraints of the form in equation (13) for those in V.This modi ed linear program only has + 2|V | variables, which is much smaller than the total number of constraints |A|, thereby maintaining a low-dimensional structure.• Step 3: We check how many constraints of X are not satis ed by the sample solution.If the sum of the weights of the constraints in X that are violated by the sample solution is more than 1/(3 − 1) fraction of the combined sum of the weights of constraints in X satis ed by the sample solution, then we discard this sample and go to Step 2. Otherwise, we double the weights of the constraints violated by the sample solution in X and go back to Step 2.
• Step 4: The algorithm terminates when the sample solution satis es all constraints in X and the number of violated constraints does not change for a xed number of iterations.

EXPERIMENTAL EVALUATION
We describe our prototype [4] for solving infeasible LDLPs, our experimental methodology for evaluating its ability to solve linear programs, and the performance of the resulting math libraries with the maximum consensus solution from our solver.Prototype and methodology.Our solver for infeasible LDLPs takes linear programs as input and produces oating point solutions in double precision for all the unknown variables.It is written in C++ and handles linear programs with several billion constraints.The prototype uses the SOPLEX solver, which uses exact rational arithmetic to solve small linear programs internally.To evaluate our infeasible LDLP solver, we use linear programs generated from the RLibm project for generating correctly rounded elementary functions.A single polynomial approximation for an elementary function from the RLibm project produces correctly rounded results for all rounding modes and multiple representations (up to 32-bits).When we create a feasible subset of the original infeasible LDLP, we use the solver for feasible linear programs from the RLibm project [2], which implements Clarkson's method [20] to produce oating point solutions.We use the computational geometry package, CGAL [56], to compute the convex hull using our 2-D projections.
We experimented with the following LP solvers to solve the linear programs generated in the RLibm project: Gurobi, CPLEX, and SOPLEX.All these solvers failed to solve any of the the LPs from the RLibm project.For the experimental evaluation, we compare our solver with RLibm's solver for feasible linear programs [2] and our implementation of the maximum consensus solution from the computer vision community [34].When given an infeasible linear program, RLibm's solver does not terminate.We modi ed it to print out the number of constraints satis ed and violated by the best solution after every iteration.We ran the RLibm's solver for 500 iterations and use the best solution produced for our evaluation.We also implemented the maximum consensus solution for small programs from the vision community for the entire LDLP [34] because there are no publicly available versions.We use the wall clock time to measure the time taken to produce a solution.
We incorporated the resulting oating point solutions generated by our solver for each elementary function into RLibm's math libraries.We validated that the resulting functions produce correctly rounded results for all representations and for all inputs.During this process, we also improved the implementation of range reduction and output compensation code in the RLibm project.We evaluate the performance bene ts from these changes along with those from the maximum consensus solutions.We have committed the resulting functions to the RLibm git repository.
Checking the validity of the resulting functions can be completed within a minute.Subsequently, we measure the performance of the resulting functions.We conducted our experiments on a 2.10GHz Intel Xeon(R) Silver 4310 server with 256GB of RAM running Ubuntu 22.04.3LTS that has both Intel turbo boost and hyper-threading disabled to minimize perturbation.We compiled our resulting math libraries with -march=native -O3 ags.For measuring performance, we use hardware performance counters using the perf tool to count the number of cycles taken to compute the result for each input.We aggregate these counts for all inputs to compute the total time taken for each elementary function.
Ability to solve infeasible linear programs and produce maximum consensus solutions.Table 1 compares the ability of the three solvers: our solver, the solver from the RLibm project, and the MCS solver from the computer vision community to solve linear programs for generating correctly rounded implementations for 16 functions.The maximum consensus solver from the vision community [34] failed to produce a solution for any of the functions.The SOPLEX solver used to solve the MCS LP formulation times out without producing a solution.This is because Table 1.This table reports whether the solvers produce the maximum consensus solution (i.e., ✓) or not (i.e., ✗).We report the type of LP and the number of RLibm constraints (i.e., ≤ ( ) ≤ ℎ ) for a particular function.We compare our solver to the RLibm project's solver [2] and our implementation of the maximum consensus solver from the computer vision community [34].We report the speedup in generating a solution with our solver in comparison to RLibm's solver when it generates an acceptable solution with less than 30 violated constraints.that formulation destroys the low-dimensional nature of the LP problem and no solver can solve LPs with billions of unknowns.In contrast, our solver is able to generate the best known maximum consensus solution for both infeasible and feasible LDLPs corresponding to these functions.
For feasible LDLPs, there exists a solution that satis es all the constraints.RLibm's solver can solve such feasible LDLPs in similar time as our solver (i.e., for log2f, cospif, sinpif, and coshf).In the case of an infeasible LDLP, there does not exist a single solution that satis es all the constraints.Our goal is to identify a solution that satis es the maximum number of constraints (i.e., only a few violated constraints).RLibm's solver cannot generate the solution produced by our solver within the speci ed number of iterations.In contrast, our solver is able generate a better solution in signi cantly faster time (174× faster for tanf).Overall, our solver is able to generate a better solution than RLibm's solver in at least 10× shorter time for infeasible LDLPs.
Analysis of the solutions from the various stages of our solver.For each function, Table 2 reports the number of canonicalized constraints, which is 2× the number of RLibm constraints.The number of constraints for each function is a function of range reduction.For trigonometric functions, one needs to use 1 with a large number of bits to produce correctly rounded results [43], which causes each original input to result in a unique reduced input (and a constraint) after range reduction.Hence, trigonometric functions have more than 2 billion canonicalized constraints.
The third, fourth, and fth columns of Table 2 report the number of extreme points on the convex hull that approximates the number of infeasible constraints (i.e., |V 1 | from Figure 5), the number of constraints violated by the solution generated by using the basis and the new LP for the nearly maximum consensus FP solution (i.e., |V | from Figure 5), and the number of constraints violated by the solution with the iterative loop for maximum consensus, respectively.The use of convex hull to over-approximate the infeasible constraints reduces the number of constraints from billions to a few thousands (i.e., 11,773 in the worst case with exp10f), which allows us to solve the new LP for maximum consensus with SOPLEX.Fig. 7. Speedup of the resulting correctly rounded functions generated using our approach.The first bar reports the speedup just due to the reduction in the number of special cases with our maximum consensus solution (MCS).The second bar reports the speedup due to reduction in the special cases along with other optimizations to range reduction and output compensation when compared to functions in the RLibm project.
For feasible LDLPs, our nearly-maximum consensus formulation produces the nal solution.For infeasible LDLPs, the solution reported by our nearly-maximum consensus solution is almost as good as the best solution that we generate.The only exception being the log10f function where our iterative loop for maximum consensus is needed to account for both the rounding errors resulting from rounding an exact rational solution to an FP solution and reducing the bias caused by a xed basis, which was computed independently of the set of violated constraints V 1 .The number of violated constraints reduces from 39 to 8 with our iterative maximum consensus loop.
The sixth column of Table 2 reports the number of violated constraints in the best solution generated by the RLibm solver.The RLibm solver did not generate the best solution reported by our solver.Finally, the seventh column of Table 2 reports the time taken to produce the best solution with our solver.The time taken increases with the increase in the number of constraints.We are able to generate the best solution in approximately 36 minutes in the worst case (for tanf), which is two orders of magnitude faster than RLibm's solver (see seventh column of Table 1).
Improvements in the performance of the resulting elementary functions.Figure 7 reports the speedup of the resulting elementary functions as a result of: (1) maximum consensus solutions generated by our solver and (2) other improvements to range reduction and output compensation along with the maximum consensus solutions.Each constraint that is violated by our solution is added as a branch condition.We use __builtin_expect intrinsics to provide hints to the compiler to optimize the common case.On average, our implementations with maximum consensus solutions are 8% faster than the corresponding RLibm functions with the same range reduction and output compensation code.The feasible LDLPs do not see any improvement in performance.The performance improvement is signi cant for infeasible LDLPs with a large number of special cases.The logf function has 24% speedup because it has a signi cant reduction in the number of constraints violated between the RLibm solution (i.e., 24) and our solution (i.e., 10).When our improvements to range reduction and output compensation are combined with maximum consensus solutions, we improved the performance of the resulting implementations by 20% on average when compared to previous RLibm functions.

RELATED WORK
We describe closely related work on solving linear programs and on computing the convex hull.
Solving Linear Programs.Detailed exposition on various approaches to solve linear programs ranging from the simplex algorithm to interior-point and ellipsoid methods can be found in the textbooks [59].As a result of algorithmic and engineering advances, modern LP solvers can easily solve several thousand constraints and unknowns.Meggido [42] and Clarkson [20] made seminal advances by observing that certain classes of "feasible" linear programs that have low-dimensions can be solved much more e ectively with randomized algorithms.Based on the advances of Clarkson, Seidel [51] observed that computing the convex hull can help in solving low-dimensional feasible linear programs with a simpler analysis.When these algorithms are applied to infeasible linear programs, they will not terminate and do not generate the maximum consensus solution.Further, all these algorithms work with real values and may not produce oating point solutions.In our prior work in the RLibm project, we used the ideas from Clarkson's method and designed a solver for feasible LDLPs that produce oating point solutions with the solve-and-re ne loop [2].The key idea in RLibm's adaptation of Clarkon's method is to solve the sample with real values and tighten the constraints when the sample solution with real values rounded to oating point values does not satisfy the sample.When given an infeasible LDLP, RLibm's solver fails to terminate and nd the maximum consensus solution similar to Clarkson's method.We use our previous RLibm solver for feasible LDLPs internally once we create a feasible subset of the infeasible LDLP.
Irreducible infeasible subsets.A subset C of constraints that itself is infeasible, but any proper subset of C is feasible, is called an irreducible infeasible subset (IIS).The concept of an IIS was rst introduced for general optimization problems and later introduced for linear programs [29,58].The rst practical methods for computing them were presented in the seminal work of Chinneck and Dravnieks [18].Since then, after identifying an IIS in a linear program, various methods have been proposed for removing constraints until the system becomes feasible [17,18,54,60].
Maximum consensus formulations.The concept of maximum feasible subsets (MAXFS) of constraints were rst introduced by Amaldi et al. [6].Some equivalent formulations for MAXFS include the minimum unsatis ed linear relation problem [5] and the minimum cardinality IIS set-covering problem [16].All these formulations are known to be NP-hard [49].Greenberg et al. [28] have proposed a mixed-integer linear program computing the MAXFS for linear constraints.The vision community has explored similar maximum consensus solutions for many model-tting problems on real-world data, which is contaminated by noise and outliers.Similar to the LDLP formulation of RLibm [2], there are several "hypothesize-and-verify" methods in computer vision as well, predominantly RANSAC [25] and its variants.They do not provide good solutions for infeasible linear programs.Hence, globally optimal algorithms for maximum consensus, such as branch-and-bound search [35,62], tree search [15], or exhaustive search [24,45] have been explored.However, they work only for small problem sizes.Our formulation is inspired by the prior MCS work [34] which decomposed the maximum consensus problem into two separate linear programs using a penalty formulation.However, this MCS formulation [34] destroys the low-dimensional property, which is crucial for solving LPs with billions of constraints.Hence, it does not solve any infeasible or feasible system as we show in our evaluation in Section 7.
Computing the convex hull.Exactly computing the convex hull of points that are embedded in -dimensional Euclidean space has an asymptotic run-time complexity of ( ) [22].The fastest known algorithm for computing the convex hull has a complexity of ( 2 + | | log ) [50], where | | is the number of faces on the convex hull.This is still prohibitively expensive when is in the order of billions of constraints.The complexity of exactly computing the convex hull was recently recognized by Müller et al. [44] in the context of formal veri cation of deep neural networks where they propose an approximation algorithm for computing the convex hull.
For identifying a superset of the infeasible constraints, we only need to identify the points that lie on the convex hull and do not require knowledge of the connectivity structure (i.e., the topology of the convex hull).Some researchers have recognized this problem to be easier than computing the convex hull and termed it as the frame problem [21].A good discussion of the di erences in complexity of the convex hull and the frame problem is provided in [31].The method proposed by Dulá et al. [21] requires a complex initialization procedure and is serial in nature, where a linear program is needed to be solved per iteration.This approach is very slow for problem sizes with billions of constraints.Thus, we designed a di erent approach that can very quickly identify a superset of the infeasible constraints by computing the convex hull in the dual space.
Correctly rounded math libraries.The seminal book [43] by Muller provides a detailed survey on generating correctly rounded math libraries.The standard approach to develop correctly rounded math libraries prior to the RLibm project was to create polynomial approximations that minimizes the maximum error across all inputs using the Remez algorithm.In contrast to minimax methods, the RLibm project makes a case for directly approximating the correctly rounded result and creating a rounding interval, which naturally leads to an LDLP.The RLibm project has shown that the freedom available to the polynomial generator is much larger than minimax methods.Hence, solving large LDLPs with billions of constraints with few violated constraints is crucial for performance of the resulting functions.Given that infeasible LDLPs are common with the RLibm project, our solver makes the resulting math libraries much faster by producing maximum consensus solutions.

CONCLUSION
We propose a new method for solving infeasible low-dimensional linear programs with billions of constraints to generate oating point solutions that satisfy the maximum number of constraints.Our key idea is to create a superset of infeasible constraints by computing the convex hull and subsequently create a new linear program with slack variables whose solution satis es the maximum number of constraints in the original LDLP.In the context of LDLPs generated in the RLibm project, our method not only produces the best solution but also does it signi cantly faster.The resulting solutions from our solver along with other optimizations to range reduction and output compensation helped us improve the performance of RLibm's math libraries by 20% on average.Our approach to solve large infeasible LDLPs with oating point solutions will likely be useful in various domains such as neural network veri cation, robotics, and computer vision, which we plan to explore in the future.

Fig. 1 .
Fig. 1. (A) The rounding interval of a correctly rounded result v2.(B) The linear constraint generated to produce a correctly rounded result for input where the interval [ , ℎ ] is 1 ULP (units in the last place).
Interpreting constraints as half spaces (B) Mapping the lower envelope of half spaces to the upper convex hull

Fig. 4 .
Fig. 4. We use geometric duality transforms to map halfspaces to points in a dual space.Then we can compute the lower envelope of the halfspaces by computing the upper convex hull of the points in the dual space.
★. Duality transforms are also order preserving, i.e., a point lies above a line if and only if the dual point ★ lies above the dual line ★ .Geometric duality allows one to reason about lines by converting the problem into another equivalent problem about points.

Table 2 .
This table reports the total number of canonical constraints, the number of extreme points on the convex hull, violated constraints in the nearly maximum consensus solution and our final maximum consensus solution, the number of violated constraints in the RLibm solution, and the time taken to solve the infeasible LDLP with our approach.