An algorithm for constructing the skeleton graph of degenerate systems of linear inequalities

Derive the quantitative predictions of constraint-based models require of conversion algorithms to enumerate and construct the skeleton graph conformed by the extreme points of the feasible region, where all constraints in the model are fulfilled. The conversion is problematic when the system of linear constraints is degenerate. This paper describes a conversion algorithm that combines the best of two methods: the incremental slicing of cones that defeats degeneracy and pivoting for a swift traversal of the set of extreme points. An extensive computational practice uncovers two complementary classes of conversion problems. The two classes are distinguished by a practical measure of complexity that involves the input and output sizes. Detailed characterizations of the complexity classes and the corresponding performances of the algorithm are presented. For the benefit of implementors, a simple example illustrates the stages of the exposition.


Introduction
Mathematical modeling in areas of science such as the physics of quantum nonlocality [1][2][3][4][5][6][7] and systems biology [8][9][10][11], frequently take the form of a system of linear inequalities in some Euclidean space. Every inequality in the system defines a half-space and the intersection of all half-spaces constitutes the polyhedral feasible region of the model. In physics, the no-signaling approach to quantum nonlocality [2][3][4]7] leads to degenerate systems of linear constraints for the correlations between the parties playing in a correlational set up [4,6,7]. In systems biology, the structural analysis approach for the mathematical modeling of a biological complex applies constraint-based methods that take the form of degenerate linear systems of inequalities [8][9][10][11].
In order to transform the model constraints into quantitative predictions, conversion methods are necessary to derive the relational structure (a skeleton graph or a network) that is conformed by the set of extreme points of the feasible region. In the no-signaling approach to quantum nonlocality, the extreme points of the polytope of correlations provide the local operations and the elementary non-local Popescu-Rohrlich channels [2][3][4] which are necessary (and in any cases sufficient) to simulate any no-signaling correlational scenario [5]. Quite similarly, the "elementary constituents", "archetypes" or "modes" in biology [12], are provided by the extreme points of the feasibility region. In both cases, the resulting output descriptions display a wide range of structural complexities, characterized by measurements as graph entropy [13,14] and graph similarity [15].
The most effective conversion methods available arise from combinatorial geometry [16][17][18][19]. However, when the half-space descriptions of the extreme points of the feasible polytope are degenerate a combinatorial explosion is produced that is the cause of stubborn difficulties for their enumeration [19][20][21][22][23].
The paper introduces an algorithm to convert a half-space description into the skeleton graph conformed by the set of extreme points of the feasible region. Using a combination of incremental [18,19] and pivotal [24,25] methods, an algorithm with a good performance to resolve degeneracy and to complete the traversal of the set of extreme points is presented. The effects of degeneracy are studied computationally for a very large number of half-space descriptions, organized into four families according to the degree of degeneracy of the input half-space descriptions and on the complexity of the output graphs.
The standard formulation of a constraint-based model is the system of linear inequalities Matrix A is real and of size m × n. The model constraints determine the entries of A and the entries of vector b 2 R m . Each one of the m constraint inequalities in Eq (1) defines a halfspace of the euclidean space of dimension n, E n . Every row of matrix A is the constraint vector a i 2 R n , with index i 2 I ¼ f0; . . . ; m À 1g. Vector a i defines the i-th feasible half-space H i ¼ fx 2 R n : a i Á x b i g. The intersection of all half-spaces i 2 I constitutes the set of feasible values of x 2 R n , which conform the convex polytope P ¼ fx 2 R n : Ax ⩽ bg. We assume P is bounded and of full affine dimension n, for which is necessary that m > n. The description of P that is provided by Eq (1) is known as a half-space description, or H-description.
However, what is physically meaningful is the combinatorial structure that is encoded in the skeleton graph G(P) of P, known as the V-description of P. The conversion of the H-description Eq (1) into the graph G(P) is accomplished when the set of vertices V = {p 2 P: p is extreme} and the set of edges E & V × V have been determined. Then, the skeleton graph G(P) = (V, E) discloses the organizational structure that is implicit in the set of linear constraints Eq (1).
Whether p is a vertex of G(P) is decided by the non-negativity of its slack vector s(p) = b − Ap and the rank of its set of active hyperplanes ZðpÞ ¼ fi 2 I : a i Á p ¼ b i g. First, a point p is in the feasible region P whenever the slack vector s(p) is non-negative. Then, p 2 P is an extreme point if, and only if, rank ðZðpÞÞ ¼ dimðspanfa i : i 2 ZðpÞgÞ ¼ n. An extreme point p is regular (or non-degenerate) if it has an active set ZðpÞ of cardianlity jZðpÞj ¼ n. The active set of a regular extreme point p, ZðpÞ, is a basis. Otherwise, when jZðpÞj ¼ n þ s and σ ! 1, the extreme point p is σ-degenerate.
The combinatorial triviality of regular extreme points does not present any difficulty to the simplex pivoting rules [16,24,25]. However, when pivoting around a σ-degenerate extreme point p, the method faces up to nþs n À Á potential bases in the active set ZðpÞ, so the exhaustive search of the neighboring points of p requires the examination of nþs n À Á simplex tableaus and a factible pivote has to be looked for by testing the n entries of m − n − σ rows of every tableau. The misery then is that a σ-degenerate vertex demands the simplex method to do a search among a multiplicity of alternatives, just for finding the neighboring points of p. The amount of searches Eq (2) may be-per vertex!-much larger than the total number of vertices of the complete skeleton graph G(P). Besides, when doing the search, the simplex method may get trapped in an endless cycle or just get stalled [22]. Several approaches have been tried out to overcome such deficiencies [20][21][22][23].
The practice of the double description method has proved its efficiency [18,19] in the determination of the extreme rays of highly degenerate polyhedral cones [26]. However, the method is not efficient for the construction of the complete skeleton graph of large degenerate systems of linear inequalities, mainly due to the very large number of tentative vertices that are produced at intermediate stages of the conversion procedure [26]. The majority of intermediate vertices are discarded at the end. To overcome this situation, we have designed a swift and compact pivoting method to determine the neighbors of extreme points, by taking as the input the extreme rays of a cone. In this way we have combined in Algorithm 5 the best of two methods: the incremental slicing of cones to defeat degeneracy [18,19], and pivoting around extreme points for a swift traversal of extreme points [24].
The incremental procedure goes slicing a cone, starting with a regular cone that is broader than and includes the target cone. The preparation of the base cone to be the input of the incremental procedure is explained in Section 1: a basis is chosen from the active set ZðpÞ and a standard algebraic method produces the extreme rays of the corresponding regular cone. The half-spaces in the active set ZðpÞ that are not part of the basis set are inserted by the incremental procedure, one-by-one, until they are exhausted and the target cone has been sculpted. The half-space insertion procedure is explained in Section 2. The explanation includes (i) the alternative combinatorial or algebraic test necessary to identify the 2-face cuts during the slicing procedure and (ii) a recording strategy that helps the algorithm to reduce the number of tests.
The extreme rays that are produced by the incremental procedure provide the scanning directions for the pivoting rule that is followed to determine the set of neighboring points of p. The pivoting rule is developed in Section 3. In Section 4 the incremental and pivoting methods are combined in Algorithm 5, which converts the system of linear inequalities into the skeleton graph conformed by the set of extreme points of the feasible region.
The computational practice in Section 5 affords understanding about the effects that the degeneracy present in the input systems has on the performance of Algorithm 5. The very large number of input systems employed in the practice of Section 5 is organized in four families that offer a controlled and distributed sampling of the complexity spectrum of the conversion problem. The amount adopted to estimate the complexity combines the average degeneracy hσi that is present in the input system and the average connectivity κ of the output graph G(P).
The family with the lowest complexity consists of regular (non degenerate) polytopes with H-descriptions produced at random [7]. The other three are one-parameter families. The family with the highest complexity consists of Birkhoff polytopes [27,28]. The other two families, with intermediate complexities, consist of no-signaling polytopes [7]. Section 5 details the characterization of the four families, produced by applying Algorithm 5.
The computational practice of Section 5 distinguishes two classes of conversion problems. A first class consists of systems of linear inequalities that have a combined complexity which becomes smaller as a function of the input size. The systems in this class convert into skeleton graphs with a number of vertices that grows faster than their vertex-connectivity, as a function of the input size. For these conversion problems (I) the CPU time consumed by Algorithm 5 is mostly applied to complete the traversal of extreme points and not to resolve degeneracy, (II) the algebraic test for 2-face cuts in Algorithm 5 is faster than the combinatorial test and (III) the incremental procedure is highly sensitive (i) to the choice of the input basis set and (ii) to the insertion order of the cutting half-spaces.
The second class of conversion problems distinguished by the computational practice in Section 5 has a combined complexity that does not decrease as a function of the input size. As the input size of the systems of linear inequalities in this class is increased, the number of vertices of their skeleton graphs increases and the vertex-connectivity does not decrease. For these conversion problems (I) the CPU time consumed by Algorithm 5 is mostly applied to resolve the degeneracy present at the H-description and not to complete the traversal of extreme points, (II) the combinatorial test for 2-face cuts in Algorithm 5 is faster than the algebraic test and (III) the incremental procedure is not sensitive (i) to the choice of the input basis set and (ii) neither to the insertion order of the cutting half-spaces.
The two complementary classes described above are detailed in Sections 5 and 6. For the benefit of the implementor we make use of a simple, but rich enough, no-signaling constraint H-description [7] as example to illustrate our exposition.

Example (Outset)
The conversion problem consists of the two-party correlations that are feasible for a no-signaling box [4] with a binary input per party and asymmetric in its outputs, producing one out of 3 and of 2 possible outcomes per party respectively. The no-signaling and non-negativity constraints on correlations [7] produce the system of linear inequalities in Table 1. The feasible polytope P NS is the intersection of m = 24 half-spaces in an Euclidean space of dimension n = 14.
Using the constraint vectors in Table 1 one verifies that the origin p = 0 is an extreme point of P NS since its active set ZðpÞ ¼ f0; . . . ; 19g has rank ðZðpÞÞ ¼ 14 and the slack vector at p = 0, s(0) = b − A p = b, is not negative. This extreme point is degenerate with s ¼ jZðpÞj À n ¼ 6 and for the simplex method it represents a multiplicity Eq (2) of μ = 2,170,560 search options. This huge value of μ is to be compared with the 6 cutting halfspaces that the double description method needs to insert, one at a time.

Regular cones as the base case
The cone described by the active set ZðpÞ of the extreme point p is the set K ZðpÞ ¼ fx 2 R n : a i Á x b i ; i 2 ZðpÞg. The cone D ZðpÞ :¼ fx 2 R n : a i Á x 0; i 2 ZðpÞg is the translation of cone K ZðpÞ to the origin, Table 1. Parameters (A, b) of the heuristic example.
a 23 ¼x 5 þx 6 þx 7 þx 8 þx 9 þx 10 þx 11 Àx 12 Àx 13 Since rank ðZðpÞÞ ¼ n, both cones K ZðpÞ and D ZðpÞ are peaked, with apices located at p and the origin, respectively. The polyhedral cone K ZðpÞ fits the feasible region P and the 1-faces (or extreme rays) of K ZðpÞ provide the directions to scan for the neighbors of p. Then, and in view of Eq (3), the first step towards the skeleton graph G(P) is to convert the half-space description of cone D ZðpÞ into its set of extreme rays X ZðpÞ , such that D ; .
The determination of the set of extreme rays of a degenerate cone is the subject matter of the next section. Meanwhile, the extreme rays of a regular cone may be obtained from its halfspace description by methods of linear algebra. A regular cone D B is the intersection of the half-spaces of a basis B ⊊ ZðpÞ. The set X B of extreme rays of D B is given in the following.
The set X B in Lemma 1.1 is the negative of the biorthogonal companion of B. Given that B ⊊ ZðpÞ we have that D B ' D ZðpÞ . The set of extreme rays of cone K B is the set X B þ p.
For a degenerate vertex p with active set ZðpÞ the double description method, discussed in the next section, produces the set X ZðpÞ of extreme rays of cone D ZðpÞ & D B by starting with the set of rays provided by Lemma 1.1 for a basis B ⊊ ZðpÞ.

Example (The base cone)
The extreme point p = 0 of the no-signaling polytope P NS has degeneracy σ = 6. The basis B ¼ f0; . . . ; 13g ⊊ ZðpÞ, through lemma 1.1, provides us with the set of rays Use the constraint vectors in Table 1 to verify the membership conditions of set Eq (4).

Incremental slicing of cones
The incremental procedure to generate the set of extreme rays X ZðpÞ of the cone at a degenerate extreme point p, which is described by the set of active planes ZðpÞ, starts with the approximate set X B , provided by lemma 1.1 for a basis B ⊊ ZðpÞ. The set B 0 ≔ ZðpÞnB is not empty by degeneracy. Then, the half-spaces remaining in B 0 are introduced one at a time. The insertion process produces a non-increasing chain of cones by eliminating the current extreme rays which are not in the feasible half-space introduced and adding the new extreme rays that are created by the half-space that has been added. When all hyperplanes in the degeneracy set B 0 have been inserted, the set of extreme rays X ZðpÞ , describing cone D ZðpÞ , is produced. Assume the insertion procedure has gone adding half-spaces from B 0 as far as to generate coneðX J 0 Þ & D B , for some J 0 ⊊ ZðpÞ. Assume the procedure is to advance one step farther by adding the half-space H k , for some k remaining in ZðpÞnJ 0 . Then, by working on the current set of rays X J 0 and the vector a k that is associated to H k , the procedure will produce the set of extreme rays X J that results from the insertion of the half-space H k . This one-plane insertion procedure defines the function The procedure represented by F begins with the partition of the current set of rays X J 0 into three sets. A first set collects the rays that are within the feasible half-space H k , the set X þ ≔ fr 2 X J 0 : a k Á r < 0g. A second set is X À ≔ fr 2 X J 0 : a k Á r > 0g, which is the set of rays outside the feasible half-space H k . The third set is X 0 ≔ fr 2 X J 0 : a k Á r ¼ 0g, which collects the rays lying on the hyperplane @H k . The rays in X + [X 0 remain extreme for the next cone D J . We have X J X þ [ X 0 . Rays in X − become unfeasible, but they are necessary to complete the set of extreme rays of cone D J .
New extreme rays in X J are the intersections of the hyperplane @H k with 2-faces of the current cone D J 0 . In order to decide whether a pair of extreme rays ρ and ρ 0 constitute a 2-face of D J 0 or not, let J 0 jr ¼ fi 2 J 0 : a i Á r ¼ 0g and let J 0 jðr; r 0 Þ ¼ J 0 jr \ J 0 jr 0 be the joint active subset of the pair (ρ, ρ 0 ) in D J 0 . Central to the incremental slicing procedure [7,19] is the following. Lemma 2.1 (Tests of colaminarity [19]). For some J ⊊ ZðpÞ with rankðJ Þ ¼ n, let ρ and ρ 0 be extreme rays of D J . The following statements about the pair (ρ, ρ 0 ) are equivalent.
When statement (a) holds we say that the pair of extreme rays (ρ, ρ 0 ) is colaminar in D J and denote the relation by ρ * ρ 0 .
The completion of X J is achieved by incorporating all the intersections that the current hyperplane @H k makes with the 2-faces of D J 0 that are framed by pairs of rays (ρ 0 , ρ)2X − × X + . When the case is that ρ * ρ 0 , the ray is extreme for the sliced cone The collection of all such rays, completes the set of extreme rays of D J , The test of colaminarity in Eq (7), corresponding to line 7 of the pseudo-code for function F in Algorithm 1, may proceed in one of two standard ways. Either by applying the algebraic test (b) in lemma 2.1, which applies methods of linear algebra, or the combinatorial test (c), which runs over the set X J 0 nfr; r 0 g [19,29].

Example (A first plane insertion)
The extreme point p = 0 of P NS has active set ZðpÞ ¼ f0; . . . ; 19g. The set X B of extreme rays of the base cone D B , for basis B ¼ f0; . . . ; 13g, was determined in Eq (5). There remains B 0 ¼ f14; . . . ; 19g half-spaces to be inserted. By inserting H 14 first, the current set X B is partitioned into the following subsets, For the base cone D B the test of colaminarity for the pairs in X − × X + may be skipped, because in a regular cone, as is the case for D B , every pair of extreme rays is colaminar (Lemma 2.2 below). The set of new rays Y 14 that is obtained by applying Formula (6) to every pair in X − × X + is shown in Table 2. The new cone in the chain D J 0 & D B , with J 0 ¼ f0; . . . ; 13; 14g, has gained six new rays and lost the two rays in X − , which became unfeasible.
The set X ZðpÞ of extreme rays of the degenerate cone D ZðpÞ is produced by the iteration of the insertion function F in Algorithm 1. The half-spaces in the degeneracy set B 0 ¼ ZðpÞnB are introduced one by one until the set is exhausted. This incremental slicing of cones constitutes function X in Algorithm 2. The base case to start the iteration, line 1 of Algorithm 2, is the set of extreme rays X B of a basis B ⊊ ZðpÞ, as is given by lemma 1.1.
Algorithm 2 Incremental slicing with standard insertion function F.

Example (Standard incremental slicing)
For the extreme point p = 0 of P NS the standard slicing function X is applied to basis B ¼ f0; . . . ; 19g, with complement B 0 ¼ ZðpÞnB ¼ f14; . . . ; 19g. The first half-space H 14 was inserted in section 2.1 already.
When the next three half-spaces in {15, 16, 17} are inserted, all pairs in X − × X + pass the test of colaminarity, producing each a new extreme ray for the subsequent cone in the chain. Consequently, the number of intermediate extreme rays goes up. It is not so for the insertion of the last two half-spaces, H 18 and H 19 . The 2-faces of the current cone that are slashed by the insertion of H 18 are identified by the test of colaminarity. This time the subset of colaminar pairs in X − × X + , shown in Table 3, is rather sparse. The last two insertions, H 18 and H 19 , make the number of intermediate rays come down. Upon conclusion, the standard function X in Algorithm 2 returns a set of 54 extreme rays for the cone D Zðp¼0Þ , which is almost four times the 14 extreme rays at a non-degenerate extreme point in dimension 14.
The computational practice has shown [19,26] that the typical behaviour of the number of intermediary rays during the slicing process is to go up and then come down. At the step corresponding to Table 3 of Example 2.2, the elements of X − listed along the left column of Table 3 become unfeasible but they give rise to as many new rays as colaminarity symbols * are entered in Table 3. Although a great deal of the intermediary rays computed Rays produced by the intersection of H 14 with every pair in X − × X + that is colaminar, ρ 0 * ρ. by function X in Algorithm 2 do not survive as extreme rays of the target cone D ZðpÞ , the standar slicing function F has to apply the test of colaminarity (either (b) or (c) in lemma 2.1) to every pair in the current set X − × X + and at every step in the chain leading to X ZðpÞ . The standard incremental slicing procedure, Algorithms 2 and 1, gets bogged down in applying an excess of colaminarity tests.
Unfortunately, we do not have a method to avoid the tests of colaminarity during the slicing procedure. What we have is a method to reduce the number of times the test is applied. The method is simple but improves considerably the CPU time of the slicing procedure. The idea is to combine a record of pairs of rays we know are colaminar, as to avoid a re-testing of pairs, with a necessary condition of colaminarity. If the necessary condition is not fulfilled, the pair is rejected and that is it. Otherwise, the pair is searched in the record. If it is not found, then the colaminarity test is applied to the pair. After a positive test, a new ray is produced and the pair is included in the record of colaminar pairs. This improvement is implemented in Algorithm 3.

Algorithm 3
Incremental cone slicing with 2-face recording and rejection test.
if jJ jðr 0 ; rÞj < n À 2: The record of pairs of colaminar rays in Algorithm 3 is initialized in line 2 to X B Â X B because the following.   The improved Algorithm 3 rejects a pair as colaminar when the following condition does not hold.
The minimum condition in Lemma 2.3 is tested in lines 7 and 8 of Algorithm 3. The incremental method we have described for the enumeration of extreme rays of cones is not effective when extended for the conversion of the complete system of linear inequalities. A weakness of the method is that the number of intermediary vertices can grow exponentially as compared to the number of true vertices. That is why we follow a pivoting method instead.

Getting the neighbouring extreme points
The rays in the set X ZðpÞ (computed by Algorithm 3) define the directions to scan around p, looking for its neighboring extreme points. Assuming the feasible region P is compact, the positive linear span of every r 2 X ZðpÞ necessarily intersects one of the hyperplanes in I nZðpÞ: the set of hyperplanes N p ðrÞ ¼ ft 2 InZðpÞ : a t Á r > 0g is not empty. Then, for t 2 N p ðrÞ the hyperplane @H t is pierced by ρ at the point q t = p + ρλ t , with The point q t with the smallest (positive) λ t is the neighboring extreme point of p along ρ. Lemma 3.1 (Pivoting around p [7]). For every r 2 X ZðpÞ let N p ðrÞ be as defined above. The set of neighboring extreme points of p is The unordered pair {p, q}, for every q 2 V p , is an element of the set of edges of p, denoted by E p . Lemma 3.1 is the core of the pivoting function P that returns the set V p of neighbors and the set E p of edges of p. The function P is defined by the pseudo-code listed in Algorithm 4.

Assembling the skeleton graph
The traversal of extreme points starts from a known extreme point p of P, which becomes the first vertex of the skeleton graph G(P). The set of vertices V is initialized as {p}. The connectivity of p in G(P) is determined by function P in Algorithm 4 by providing us with two sets: the set of neighboring vertices V p and the set of edges E p . The set of edges E of the skeleton is initialized as E p . All the neighboring vertices of p are awaiting to be scanned. They take a place in the queue Q, initialized with the set V p . With this provision, an exhaustive search proceeds by applying repeatedly the following three steps, finishing when there are no vertices awaiting in the queue.
1. Pick the next vertex v that is available in Q and remove v from Q. When the procedure stops, the pair (V, E) is the graph G(P). The procedure defines function G, which is shown in Algorithm 5.

Algorithm 5 Assembling the skeleton graph G(P).
Pop out next v 2 Q and update Q = Q \ {v} 5 Choose a basis B v and let

Example
Function G in Algorithm 5 starts at the extreme point p = 0 to produce the skeleton graph G(P NS ). The output graph has 108 vertices and 1,548 edges, which represent a 27% of the edges of the complete graph K 108 . Two thirds (72) of the 108 vertices are sparsely connected, having degree 16. The other 36 vertices are densely connected, having degree 54 (connecting with a half of all the vertices in the graph). To appreciate how densely connected is the output graph G(P NS ), compare the 108 vertices it has against the upper bound 16,016 for regular polytopes described by the same number of half-spaces and in the same Euclidean space [30].

Computational practice
Algorithm 5 converts systems of linear inequalities into the skeleton graph G(P) of the feasible polytope P. The combinatorial complications introduced by degeneracy are explored computationally on a wealth of systems of linear inequalities, organized in four families of varying complexity.
In our scheme the input systems to Algorithm 5 are sized by the product variable z = nmhZi, which includes (a) the dimension (or number of constrained variables) n, (b) the number m of half-spaces (or constraints) and (c) the average number hZi = n + hσi of active hyperplanes of the extreme points of the feasible region P. The size of the output graph G(P) is measured by (a) the number of vertices |V| of the graph and (b) its connectivity κ = hXi/|V|, where hXi is the average number of edges attached to a vertex of G(P). The number of vertices |V| weights the bulk of the graph and the connectivity κ < 1 is a simple measure of the complexity of the graph's topology. The two output variables, κ and |V|, and the input size z provide a practical characterization of the complexity of a conversion problem, including degeneracy and its consequences.
Algorithm 5 translates the input size into some form of complexity of the output graph. Extreme simplicity is reached by the family of non-degenerate half-space descriptions that converts into regular graphs. The regular systems we are including in our practice have m constraint hyperplanes tangent to the (n − 1)-sphere in the Euclidean space E n . The points of tangency were generated at random and equally distributed on the positive ortant of the sphere. The family of regular systems labeled k0 in Fig 1 (represented by hollow diamonds) covers the ranges 7 n 20 and 16 m 64. The random systems so produced are regular with probability 1. For any value of the input size z, the non-degenerate systems of linear inequalities produce the bulkiest, Fig 1B, and most sparsely connected, Fig 1A, output graphs. The trend followed in Fig 1A by the non-degenrate k0 family implies that degenerate half-space descriptions have combined output-input size values κz 2 ! 2 × 10 4 .
The zones of intermediate complexity in Fig 1A and 1B are populated by bipartite no-signaling system [7], represented by circular dots. The star is the tripartite no-signaling binary system [6]. The bipartite system with the smallest input size z corresponds to Bell's experimental setup [1] and the hollow dot is the heuristic example in Table 1. The upper sequence of nosignaling systems in Fig 1A (red online) |referred to as the k2 family| emerges from Bell's setup by increasing the number v of outcomes for one of the parties only. The output graph has |V| = 2v 2 (2 + (v − 1) 2 ) [7].
The lower sequence of no-signaling systems in Fig 1A (green online) |referred to as the k1 family| emerges from Bell's setup by increasing the number of input options available to one of the parties. The degenerate systems in the k1 family produce the bulkiest but most sparsely connected class of no-signaling graphs.
Each sequence of no-signaling systems in Fig 1 (red and green online) constitutes a oneparameter family. The blue dots dispersed between the two border sequences correspond to binary systems differing from Bell's setup in at least two parameters. Fig 1 confirms that the An algorithm for degenerate systems of linear inequalities no-signaling polytopes constitute a good example of highly degenerate half-space descriptions that produce output graphs with a large volume |V| and a moderately dense connectivity κ.
Birkhoff systems of half-spaces produce densely connected skeleton graphs, which appear in the upper part of the output-input map of Fig 1A. They are referred to as the k3 family of ℓ × ℓ doubly-stochastic matrices [p ij ] and are represented by square dots in Fig 1 (orange  online). The polyhedral region P ℓ of doubly-stochastic matrices is known as Birkhoff's polytope. It is a notable polytope in various branches of mathematics [27,28]. A polytope P ℓ is the feasible region, in Euclidean space of dimension n = (ℓ − 1) 2 , of the non-negativity constraints p ij ! 0, subjected to the normalization conditions P ' i¼1 p ij ¼ P ' i¼1 p ji ¼ 1. The graph G(P ℓ ) has |V| = ℓ! vertices [27]. The output-input maps in Fig 1 show that Birkhoff's polytopes have a highly degenerate half-space description and for high values of the input size z their output graphs become bulky while keeping a dense connectedness (the diameter of G(P ℓ ) is 2 for every ℓ [27]). The one-parameter k3 family of Birkhoff's polytopes constitute our most complex exemplar. Fig 2A shows the practical estimate of complexity of the output graph G(P) that is afforded by volume |V| and connectivity κ. Each family follows a well defined "complexification path". The arrows in the figures indicate the direction the input size z of the half-space descriptions becomes greater. The general trend, as the input size gets bigger, is an increase of the volume |V| at the cost of loosing the connectivity κ: degeneracy at the input is converted into connectivity at different rates per family. This fact is better appreciated in Fig 2B by using the combined output-input size κz, instead of simply κ. The output graphs of the k1 family of systems exhibit a considerable increase in volume but degeneracy is not producing a dense connectivity. Fig 2B shows that family k2 exhibits the slowest growing rate of the volume of the graph and the combined output-input size κz is maintained around the value κz * 2.1 × 10 3 . This constant value of κz presents itself as a boundary between families.
The time Algorithm 5 takes to output the graph of a regular system is orders of magnitude longer than the time needed to produce the graph of a degenerate system of the same input size z. For large values of z the CPU time for regular systems is Oðz 9:6 Þ, shown as a solid segment in Fig 3A. In Fig 3A Algorithm 5 shows the best performance for the k2-family of no-signaling systems. The k2 family reaches in our exploration the biggest values of the input size z * 10 6 and the CPU time fits the law Oðz 3 Þ exactly. At z * 10 4 , Algorithm 5 consumes a CPU time to produce the skeleton graph of a k2-family degenerate system that is five orders of magnitud shorter than the CPU time needed for a regular system. The CPU times for the k1 and k3 families go, for large values of z in Fig 3A, as Oðz 7:2 Þ and Oðz 5 Þ, respectively.
Differences in performance of Algorithm 5 are better appreciated in terms of the combined output-input variable κz. Fig 3B shows a clear distinction between the four families of systems. Families k0 and k1 produce the most sparsely connected graphs and exhibit a similar growing rate of the CPU times as the combined variable κz gets smaller. Algorithm 5 shows the opposite behavior on systems of family k3: the CPU time grows with κz. In conoclusion, Algorithm 5 has a better performance on degenerate families with densely connected graphs |families k2 and k3.
The conversion algorithm involves two procedures that manage complex inputs. One is in line 6 of Algorithm 5 that resolves degeneracy to determine (a) the extreme rays of cones and (b) the neighboring extreme points. The other procedure is the exhaustive traversal of extreme points. The traversal procedure is Algorithm 5 itself, excluding the time employed by line 6. Next test decides whether the CPU time consumed by Algorithm 5 is either employed in the resolution of degeneracy or in traversing the set of extreme points.
The fraction of total CPU time employed to traverse the set of extreme points is shown in Fig 4A as  Degeneracy is the source of a combinatorial explosion of the search universe of simplexbased methods [25]. In contrast, degeneracy furnishes the incremental procedure with a huge set of options for the input ðB; B 0 Þ in lines 5 and 6 of Algorithm 5. The alternatives consist in choosing a basis B & ZðpÞ, an order for the n elements in B and an order for the σ elements in B 0 ¼ ZðpÞnB.
In our computational exploration we found that the incremental procedure is not sensitive to the order given to the elements of any chosen basis B. While for highly degenerate systems with less than moderately connected output graphs the choice of a basis B and the order adopted for B 0 are critical. The example system in Table 1 (in dimension n = 14) will be used to illustrates the situation.
The extreme point p = 0 of the system in Table 1 has jZðpÞj ¼ 20 and σ = 6. At this moderate level of degeneracy there is a set of 6144 bases, which represent a 16% of all possible  Fig  5B is the distribution for a bad basis having T =T ¼ 2:10. The difference between a good and a bad choice of the combination basis-insertion-order may represent in this example a factor greater than 2 in the CPU time.
A similar treatment for the three families of degenerate systems produces the distribution spreads T =T shown in Fig 6, plotted as a function of the combined variable κz. Degenerate systems that convert into densely connected graphs appear to have narrow spreads, T =T % 1:2 and T =T < 2:5 for the k3 and k2 families, respectively. The choice of an insertion order is not an issue for them. An algorithm for degenerate systems of linear inequalities Degenerate systems that convert into bulky but sparsely-connected graphs show in Fig 6 a much higher sensitivity to the election of basis and insertion order. The spread T =T of the dispersion of CPU times may be as wide as 24 for the cases used in our exploration and the spread keeps growing with z. Choosing a suitable combination of basis B and insertion order of B 0 is an issue for the optimization of the execution times of degenerate systems that convert into sparsely connected graphs.
The incremental slicing of cones may apply either a combinatorial or an algebraic test for 2-face cuttings. The results of a comparison of performance of the algebraic (T ) and the combinatorial (T) tests are shown in Fig 7. The combinatorial test shows a better performance on systems that convert into densely connected graphs, the k2-family has the ratio T=T $ 0:32. For systems that convert into sparsely connected graphs the algebraic test is the best option, specially for big input sizes. For the k3-family the combinatorial test is better for moderately bulky graphs, but when they get bigger the algebraic test turns out to be the better choice. Fig 8 compares the execution times of Algorithm 5 with the times employed by our own implementation of the pure-simplex based Balinski's algoritm [25]. Only systems that were

Concluding remarks
The algorithm presented converts degenerate systems of linear inequalities into their skeleton graphs. Algorithm 5 applies an improved version of the incremental method for the enumeration of extreme rays to defeat degeneracy and a simple pivoting rule for a swift traversal of the set of extreme points.
The results obtained by Algorithm 5 in the computational practice of Section 5 characterize conversion problems in two classes, which are distinguished by the combined output-input measure of complexity κz. Class A is constituted by systems of linear inequalities that convert degeneracy into densely connected graphs. Systems that convert degeneracy into bulky but The combinatorial test is the best option for systems that convert into densely-connected graphs while the algebraic test is appropriate for systems converting into sparselyconnected graphs. https://doi.org/10.1371/journal.pone.0175819.g007 An algorithm for degenerate systems of linear inequalities sparsely connected graphs conform the second class B. The almost stationary value of κz shown by the k2 family of no-signaling polytopes suggests that class B systems have values of κz ≲ 10 3 .
The computational practice in Section 5 showed that Algorithm 5 performs better on systems of class A. For these systems, the incremental method is neither sensitive to the choice of the input basis nor to the insertion order of the cutting half-spaces and there is a clear evidence that favors the combinatorial over the algebraic 2-face test.
On the contrary, for class B systems the performance of the incremental method is highly sensitive to the large number of combinatorial options introduced by degeneracy. The selection of an "optimal" input basis and an "optimal" insertion order remains a problem for the systems in class B. No solution to this problem is found in the literature, only recommendations to achieve some technical easiness, such as keeping a fixed insertion order and the remark that any dynamical reordering based on explorations consumes a longer CPU time [19,29,31]. The recommendation for class B derived from the computational An algorithm for degenerate systems of linear inequalities practice of Section 5 is to explore a few extreme points in advance as to find a good order of the constraint vectors for the sample of points and then adopt the finding for the whole procedure.
Concerning the 2-face test, the algebraic one is undoubtedly the best option for class B systems. In practice, it is advisable to run a competition of the combinatorial and algebraic tests on a small sample of extreme points and use the winning test to proceed with the full conversion.
In the extensive practice of Section 5, Algorithm 5 showed a superior performance respect the simplex-based pivoting algorithm by Balinski [25]. Our own implementation of Balinski's algorithm got stalled (did not finish) when trying to convert most of the systems studied in Section 5.
In this paper we have used the combined output-input variable κz as the measure of complexity for the output graphs. However, there exist other measures for characterizing such complexity, as the graph entropy [13,14] and the graph similarity [15]. Then, a future work that involves such measures to characterize the output graphs produced by degenerate systems of linear restrictions would be interesting.
The algorithm introduced in this paper aims to be a useful tool in applied problems that require a conversion mechanism for an appropriate interpretation. We are aware that a theoretical analysis of its computational complexity is required. However, for an analysis of computational complexity to be of some significance, a precise and definite characterization of the input descriptions is required. The problem is that, unlike regular input systems, degenerate graphs cannot be described simply by their dimension and number of half-spaces, as shown in Section 5. We did not find in the literature a standard characterization for the input descriptions as to typify degenerate systems. The no-signaling and Birkhoff's half-space descriptions are candidates to consider, but first one of them should be characterized as the specific standard for degeneracy. On the other hand, the effect of implementations that optimize the performance of the algorithm (such as the 2-face recording and the rejection test) must be differentiated and evaluated too. The above considerations tells us that, even necessary, a significant analysis of computational complexity is a very ambitious endeavor that lays beyond the scope of this article.