Jean-Pierre Dussault
email: pierre.dussault@usherbrooke.ca

Jean Charles Gilbert

Baptiste Plaquevent-Jourdain
email: baptiste.plaquevent-jourdain@usherbrooke.ca

J.-P Dussault

J Ch

Gilbert

On the B-differential of the componentwise minimum of two affine vector functions

Keywords: Mathematics Subject Classification (2020) 05A18, 05C40, 26A24, 26A27, 46N10, 47A50, 47A63, 49J52, 49N15, 52C35, 65Y20, 65K15, 90C33, 90C46

This paper focuses on the description and computation of the B-differential of the componentwise minimum of two affine vector functions. This issue arises in the reformulation of the linear complementarity problem with the Min C-function. The question has many equivalent formulations and we identify some of them in linear algebra, convex analysis and discrete geometry. These formulations are used to state some properties of the B-differential, like its symmetry, condition for its completeness, its connectivity, bounds on its cardinality, etc. The set to specify has a finite number of elements, which may grow exponentially with the range space dimension of the functions, so that its description is most often algorithmic. We first present an incremental-recursive approach avoiding to solve any optimization subproblem, unlike several previous approaches. It is based on the notion of matroid circuit and the related introduced concept of stem vector. Next, we propose modifications, adapted to the problem at stake, of an algorithm introduced by Rada and Černý in 2018 to determine the cells of an arrangement in the space of hyperplanes having a point in common. Measured in CPU time on the considered test-problems, the mean acceleration ratios of the proposed algorithms, with respect to the one of Rada and Černý, are in the range 15..31, and this speed-up can exceed 100, depending on the problem and the approach. Keywords B-differential • bipartition of a finite set • C-differential • complementarity problem • componentwise minimum of functions • connectivity • dual approach • Gordan's alternative • hyperplane arrangement • matroid circuit • pointed cone • Schläfli's bound • stem vector • strict linear inequalities • Winder's formula.

Introduction

Let E and F be two real vector spaces of finite dimensions n := dim E and m := dim F. The B-differential at x ∈ E of a function H : E → F is the set denoted and defined by ∂ B H(x) := {J ∈ L(E, F) : H ′ (x k ) → J for a sequence {x k } ⊆ D H converging to x}, (1.1) where L(E, F) is the set of linear (continuous) maps from E to F and D H is the set of points at which H is (Fréchet) differentiable (its derivative at x is denoted by H ′ (x)). Recall that a locally Lipschitz continuous function is differentiable almost everywhere in the sense of the Lebesgue measure (Rademacher's theorem [START_REF] Rademacher | Über partielle und totale differenzierbarkeit[END_REF]) and this property has the consequence that the B-differential of a locally Lipschitz function is nonempty and bounded everywhere [START_REF] Frank | Optimization and Nonsmooth Analysis[END_REF]. The B-differential is an intermediate set used to define the C-differential (C for Clarke [START_REF] Frank | Optimization and Nonsmooth Analysis[END_REF]) of H at x, which is denoted and defined by

∂ C H(x) := co ∂ B H(x), (1.2)
where co S denotes the convex hull of a set S [START_REF] Rockafellar | Convex Analysis[END_REF][START_REF] Hiriart-Urruty | Fundamentals of Convex Analysis[END_REF][START_REF] Michael Borwein | Convex Analysis and Nonlinear Optimization -Theory and Examples[END_REF]. Both intervene in the specification of conditions ensuring the local convergence of the semismooth Newton algorithm [START_REF] Qi | Convergence analysis of some algorithms for solving nonsmooth equations[END_REF][START_REF] Qi | A nonsmooth version of Newton's method[END_REF][START_REF] Marek | On a new exponential iterative method for solving nonsmooth equations[END_REF], which can be a motivation for being interested in that concept.

In this paper, we focus on the description of the B-differential of H at x when H : R n → R m is the componentwise minimum of two affine functions x → Ax + a and x → Bx + b, where A, B ∈ R m×n and a, b ∈ R m . Hence, H is defined at x by H(x) = min(Ax + a, Bx + b), (1.3) where the minimum operator "min" acts componentwise (for two vectors u, v ∈ R m and i ∈ [1 : m] := {1, . . . , m}: [min(u, v)] i := min(u i , v i )). This function is usually nonsmooth. A motivation to look at the B-differential of that function H comes from the fact that, when m = n and H is given by (1.3), as explained below, the equation H(x) = 0 This system expresses the fact that a point x ∈ R n is sought such that Ax + a 0, Bx + b 0 and (Ax + a) T (Bx + b) = 0 (the superscript " T " is used here and below to denote vector or matrix transposition). Problem (1.5) is a special case of the so-called (extended ) vertical LCP, which uses more than two matrices and vectors in its formulation [22,[START_REF] Sznajder | The generalized order linear complementarity problem[END_REF][START_REF] Zhang | Global error bounds for the extended vertical LCP[END_REF]. In the standard LCP, A is the identity matrix and a = 0 [START_REF] Katta | Linear Complementarity, Linear and Nonlinear Programming[END_REF]23].

The reformulation (1.4) of (1.5) is based on the fact that, for two real numbers α and β, min(α, β) = 0 if and only if α 0, β 0 and αβ = 0 [1,[START_REF] Pang | A B-differentiable equation-based, globally and locally quadratically convergent algorithm for nonlinear programs, complementarity and variational inequality problems[END_REF]. This reformulation serves as the basis for a number of solving methods and investigations [1,[START_REF] Kojima | Extension of Newton and quasi-Newton methods to systems of PC 1 equations[END_REF][START_REF] Pang | Newton's method for B-differentiable equations[END_REF][START_REF] Pang | A B-differentiable equation-based, globally and locally quadratically convergent algorithm for nonlinear programs, complementarity and variational inequality problems[END_REF][START_REF] Pang | Complementarity problems[END_REF][START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems (two volumes)[END_REF][START_REF] Ibtihel | Nonconvergence of the plain Newton-min algorithm for linear complementarity problems with a P -matrix[END_REF]9,[START_REF] Alexey | Newton-Type Methods for Optimization and Variational Problems[END_REF][START_REF] Ibtihel | An algorithmic characterization of P-matricity II: adjustments, refinements, and validation[END_REF][START_REF] Dussault | A lower bound on the iterative complexity of the Harker and Pang globalization technique of the Newton-min algorithm for solving the linear complementarity problem[END_REF]27,28]. If

(1.5) stands alone, it is appropriate to have m = n, but (1.5) may be part of a system with other constraints to satisfy [START_REF] Marchand | Fully coupled generalised hybrid-mixed finite element approximation of two-phase two-component flow in porous media. Part II: numerical scheme and numerical results[END_REF][START_REF] Marchand | Fully coupled generalised hybrid-mixed finite element approximation of two-phase two-component flow in porous media. Part I: formulation and properties of the mathematical model[END_REF][START_REF] Ibtihel Ben Gharbia | Gas phase appearance and disappearance as a problem with complementarity constraints[END_REF], in which case m n. In the computation of the B-differential of the Min function (1.3), m and n may be unrelated. Note that there are many other ways of reformulating problem (1.5) as a nonsmooth system of equations. It is frequent to use the Fischer function, whose B-differential is computed in [START_REF] Facchinei | A new merit function for nonlinear complementarity problems and a related algorithm[END_REF]. The function H in (1.3) has been less studied and used than the Fischer function, although it has various advantages: it is piecewise affine (but has more nondifferentiability kinks), the local convergence of a semi-smooth Newton algorithm using it can be established under weaker assumptions and may be finitely locally convergent for linear complementarity problems [36; § 9.2].

Occasionally, we shall refer to the nonlinear version of the above problem, in which a function H : E → R m is defined at x ∈ E by H(x) := min(F (x), G(x)), (1.6) where F and G : E → R m are two functions and the "min" operator still acts componentwise. The equation H(x) = 0 is then a reformulation of the complementarity problem "0 F (x) ⊥ G(x) 0".

As a first general remark, let us quote the fact that the B-differential of H cannot be deduced from the knowledge of the B-differential of its scalar components

H i : x ∈ E → H i (x) ∈ R, for i ∈ [1 : m],
which is trivial in the present context. Indeed, it is known that [20; proposition 2.6.2(e)]

∂ B H(x) ⊆ ∂ × B H(x) := ∂ B H 1 (x) × • • • × ∂ B Hm(x),
(1.7) but equality in this inclusion may not hold (see [36; § 7.1.15], counter-example 2. 3 and almost all the examples and test-cases below). Therefore, all the components of H must be taken into account simultaneously.

The B-differential of H at x is a finite set, made of Jacobians whose ith row is A i, : or B i, : (proposition 2.2). Consequently, its cardinality can be exponential in m and it occurs that its full mathematical description is a tricky task, essentially when there are many indices i for which (Ax + a) i = (Bx + b) i and A i, : = B i, : , a situation that makes H nondifferentiable (lemma 2.1).

Then, a rich panorama of configurations appears, which is barely glimpsed in this contribution.

The paper starts with a background section (section 2), which recalls a basic property of the minimum of two functions (lemma 2.1) and gives us a first perception of the structure of the Bdifferential of the function H, in particular its finite nature (proposition 2.2). A useful technical lemma is also presented (lemma 2.6).

In section 3, it is shown that the problem of computing ∂ B H(x) has a rich panel of equivalent formulations, related to various areas of mathematics. We have quoted two forms of the problem in linear algebra, which are dual to each other (section 3.2), two equivalent problems in convex analysis (section 3.3) and a last equivalent problem, which arises in computational discrete geometry and deals with the arrangement of hyperplanes having the origin in common (section 3.4).

Section 4 gives some properties of the B-differential of H, recalls Winder's formula of its cardinality, provides some lower and upper bounds on this one, proves necessary and sufficient conditions so that two extreme configurations occur and highlights two links between the B-differential and C-differential.

Section 5 presents algorithms for computing one (section 5.1) or all (section 5.2) the Jacobians of ∂ B H(x). In the latter case, the algorithms construct a tree incrementally and recursively (section 5.2.1), as proposed by Rada and Černý [62]. On the one hand (section 5.2.2), an algorithm based on the notion of matroid circuit of the matrix V expressing the "derivative gap" is proposed; it has the nice feature of requiring no linear optimization problem (LOP) to solve. On the other hand (section 5.2.4), various modifications of the algorithm of Rada and Černý [62] are proposed with the goal of decreasing the number of LOPs to solve. Numerical experiments are reported (section 5.2. [START_REF] Baldi | Deep learning in biomedical data science[END_REF], showing that the proposed algorithms significantly improve the performance of the Rada and Černý method, with mean (resp. median) acceleration ratios in the range 15..31 (resp. 5..20), measured by the computing time. This speed-up exceeds 100, for some algorithms and test-problems.

This paper is an abridged version of the more detailed report [29].

Notation. We denote by |S| the number of elements of a set S (i.e., its cardinality). The power set of a set S is denoted by P(S). The set of bipartitions (I, J) of a set K is denoted by B(K): I ∪ J = K and I ∩ J = ∅. The sets of nonzero natural and real numbers are denoted by N * and R * , respectively. The sign of a real number is the multifunction sgn : R ⊸ R defined by sgn(t) = {1} if t > 0, sgn(t) = {-1} if t < 0 and sgn(0) = [-1, 1]. We note R n + := {x ∈ R n : x 0} and R n ++ := {x ∈ R n : x > 0} (strict inequalities must also be understood componentwise; hence x > 0 means x i > 0 for all indices i). For a subset S of a vector space, we denote by vect(S) the subspace spanned by S. The vector of all one's, in a real space whose dimension is given by the context, is denoted by e. The Hadamard product of u and v ∈ R n is the vector u q v ∈ R n whose ith component is u i v i . The range space of an m × n matrix A is denoted by R(A), its null space by N (A), its rank is rank(A) := dim R(A) and its nullity is null(A) := dim N (A) = nrank(A) by the rank-nullity theorem. The ith row (resp. column) of A is denoted by A i, : (resp. A : ,i ). Transposition operates after a row/column selection: A T i, : is a short notation for the column vector (A i, : ) T and A T : ,i is a short notation for the row vector (A : ,i ) T . For a vector α, Diag(α) is the square diagonal matrix with the α i 's on its diagonal.

Background

Recall that F : E → F is said to be (Fréchet) differentiable at x if F (x + d) = F (x) + Ld + o( d ) for some L ∈ L(E, F), in which case one denotes by F ′ (x) = L the derivative of F at x. We say below that F is continuously differentiable at x if it is differentiable near x (like in [START_REF] Frank | Optimization and Nonsmooth Analysis[END_REF], "near" means here and below "in a neighborhood of" in the topological sense) and if its derivative is continuous at x.

The next famous lemma recalls a necessary and sufficient condition guaranteeing the differentiability of the minimum of two scalar functions (see [60; 1993, final remarks (1)], [78; 2011, theorem 2.1] and [29]).

Lemma 2.1 (differentiability of the Min function) Let f and g : E → R be two functions and h : E → R be defined by h(•) := min(f(•), g(•)). Suppose that f and g are differentiable at a point

x ∈ E. 1) If f (x) < g(x), then h is differentiable at x and h ′ (x) = f ′ (x). 2) If f (x) > g(x), then h is differentiable at x and h ′ (x) = g ′ (x). 3) If f (x) = g(x), then h is differentiable at x if and only if f ′ (x) = g ′ (x). In this case, h ′ (x) = f ′ (x) = g ′ (x).
The previous lemma shows the relevance of the following index sets:

A(x) := {i ∈ [1 : m] : (Ax + a) i < (Bx + b) i },
(2.1a)

B(x) := {i ∈ [1 : m] : (Ax + a) i > (Bx + b) i },
(2.1b)

E (x) := {i ∈ [1 : m] : (Ax + a) i = (Bx + b) i },
(2.1c)

E = (x) := {i ∈ E (x) : A i, : = B i, : },
(2.1d)

E = (x) := {i ∈ E (x) : A i, : = B i, : }.
(2.1e)

To simplify the presentation, we assume in the sequel that

E = (x) = [1 : p],
(2.2)

for some p ∈ [0 : m] (p = 0 if and only if E = (x) = ∅).
The next proposition describes the superset

∂ × B H(x) of ∂ B H(x)
given in the right-hand side of (1.7) (see [47; 1998, § 2] in a somehow different context, [25; 2000, before (8)] and [29] for a meticulous proof). This Cartesian product actually reads

∂ × B H(x) := {J ∈ L(E, R m ) : J i, : = A i, : , if i ∈ A(x), J i, : = B i, : , if i ∈ B(x), J i, : = A i, : = B i, : , if i ∈ E = (x), J i, : ∈ {A i, : , B i, : }, if i ∈ E = (x)}. (2.3) Proposition 2.2 (superset of ∂ B H(x)) One has ∂ B H(x) ⊆ ∂ × B H(x) = ∂ B H 1 (x) × • • • × ∂ B Hm(x). In particular, |∂ B H(x)| 2 p .
The following counter-example shows that one can have ∂ B H(x) = ∂ × B H(x) and highlights the interest of the B-differential for the convergence of the semismooth Newton algorithm on (1.4).

Counter-examples 2.3 Let n = 2, m = 2, A = ( -1 1 -1 -1 ), B = ( 1 1 1 -1 ) and a = b = ( 0 0 ). One has A(0) = B(0) = ∅, E (0) = E = (0) = {1, 2}, ∂ B H(0) = {A, B}, while ∂ × B H(0) = {A, B, ( -1 1 1 -1 ), ( 1 1
-1 -1 )}. This example also shows that all the Jacobians of ∂ B H(0) can be nonsingular, while the Jacobian

( 1 1 -1 -1 ) of ∂ × B H(0)
is singular and the central Jacobian (4.10), namely

1 2 (A + B) = ( 0 1 0 -1 ) ∈ ∂ C H(0)
, is also singular. Therefore, in this case, H is strongly BD-regular at 0 in the sense of [60; p. 233] and the conditions ensuring the local convergence of the semismooth Newton algorithm are satisfied [60; theorem 3.1].

⊓ ⊔

The previous proposition shows that ∂ B H(x) is a finite set. It also naturally leads to the next definition.

Definition 2.4 (complete B-differential) We say that the B-differential of H at x ∈ R n is complete if ∂ B H(x) = ∂ × B H(x) or, equivalently, if |∂ B H(x)| = 2 p . ⊓ ⊔ Definitions 2.5 (symmetry in ∂ B H(x)) For x ∈ E, we say that the Jacobian J ∈ ∂ × B H(x) is symmetric to the Jacobian J ∈ ∂ × B H(x) if Ji, : = A i, : if i ∈ E = (x) and J i, : = B i, : , B i, : if i ∈ E = (x) and J i, : = A i, : . The B-differential ∂ B H(x) itself is said to be symmetric if each Jacobian J ∈ ∂ B H(x) has its symmetric Jacobian J in ∂ B H(x).

⊓ ⊔

We shall use several times the following lemma, which, for the sake of generality, is written in a slightly more abstract formalism than the one we need below (one could take the Euclidean scalar product for

: i ∈ [1 : p]}| = p is dense in E. Proof Denote by Ξ the set of vectors ξ ∈ E such that |{ ξ, v i : i ∈ [1 : p]}| = p (i.e., { ξ, v i : i ∈ [1 : p]} has p distinct values in R). We have to show that Ξ is dense in E.
Take ξ 0 / ∈ Ξ, so that ξ 0 , v i = ξ 0 , v j for some i = j in [1 : p]. By continuity of the scalar product, for any ε 0 > 0 sufficiently small, the vector

ξ 1 := ξ 0 -ε 0 (v i -v j ) guarantees ξ 1 , v i1 < ξ 1 , v i2 for all i 1 and i 2 ∈ [1 : p] such that ξ 0 , v i1 < ξ 0 , v i2 (
in other words, ξ 1 maintains strict the inequalities that are strict with ξ 0 ). In addition

ξ 1 , v i -ξ 1 , v j = ξ 0 , v i -v j =0 -ε 0 v i -v j 2 >0 < 0.
Therefore, one gets one more strict inequality with ξ 1 than with ξ 0 . Pursuing like this, one can finally obtain a vector ξ in Ξ. This vector is arbitrarily close to ξ 0 by taking the ε i 's positive and sufficienty small.

⊓ ⊔

Equivalent problems

The problem of determining the B-differential of the piecewise affine function, that is the minimum (1.3) of two affine functions, appears in various contexts, sometimes with non straightforward connections with it (this one is recalled in section 3.1). We review some equivalent formulations in this section (see also [START_REF] Winder | Partitions of N-space by hyperplanes[END_REF][START_REF] Avis | Reverse search for enumeration[END_REF]7] and the references therein) and give a few properties of the Bdifferential in this piecewise affine case. As suggested by proposition 2.2, these problems have an enumeration nature, since a finite list of mathematical objects has to be determined. This list may have a number of elements exponential in p, which makes its content difficult to specify (in this respect, the particular case where the B-differential is complete is a trivial exception). Some formulations, such as the one related to the arrangement of hyperplanes containing the origin (section 3.4), have been extensively explored, others much less. Each formulation sheds a particular light on the problem and is therefore interesting to mention and keep in mind. They also offer the possibility of introducing new algorithmic approaches to describe the B-differential.

B-differential of the minimum of two affine functions

The problem of this section was already presented in the introduction and is sometimes referred to, in this paper, as the original problem.

⊓ ⊔

When E = (x) = ∅, the rows of B -A with indices in E = (x) will play a key role below. We denote its transpose by

V := (B -A) T E = (x), : ∈ R n×p . (3.1)
Note that, due to their indices in E = (x) = [1 : p] and the definition of this index set, the columns of V are nonzero. This matrix may not always have full rank, however.

The following example will accompany us throughout this section.

Example 3.2 (a simple example) Consider the trivial linear complementarity problem 0 x ⊥ (Mx + q) 0 defined by

M =   2 0 0 -α 1+β 0 -α -β 1   and q = 0,
where α :=cos(2π/3) = 1/2 > 0 and β := sin(2π/3) ∈ (α, 2α). Note that, at the unique solution x = 0 to the problem, one has

A(x) = B(x) = E = (x) = ∅ and E (x) = E = (x) = [1 : 3], so that p = 3 and V =   1 -α -α 0 β -β 0 0 0   . ⊓ ⊔

Linear algebra problems

Signed feasibility of strict inequality systems

We call sign vector a vector whose components are +1 or -1. Many proofs below leverage the equivalence between the original problem 3.1 and the following one. The reason is that working on problem 3.3 often allows us to propose shorter proofs. In addition, the algorithms of section 5 all focus on the generation of the sign vectors s forming the set S in (3.2) below. Recall the definition of the Hadamard product:

(u q v) i = u i v i .
Problem 3.3 (signed feasibility of strict inequality systems) Let be given two positive integers n and p ∈ N * and a matrix V in R n×p with nonzero columns. It is requested to determine the set

S := {s ∈ {±1} p : s q (V T d) > 0 holds for some d ∈ R n }. (3.2) ⊓ ⊔
By routine verification, one can see that the sign vectors s in S for example 3.2 are given by the columns of the matrix S below and possible associated directions d such that s q (V T d) > 0 are given by the corresponding columns of the matrix D:

S =   1 1 1 -1 -1 -1 1 -1 -1 -1 1 1 -1 -1 1 1 1 -1   and D =   2 2 2 -2 -2 -2 2 1 -2 -2 -1 2 0 0 0 0 0 0   . (3.3)
The sign vectors ±e := ±(1, 1, 1) are not in S since V e = 0 (there is not d ± such that (±e) q (V T d ± ) > 0, since this would imply that 0 < ±e T V T d ± = 0, a contradiction). Therefore, there are only 6 sign vectors in S instead of the 8 sign vectors in {±1} 3 .

The link between problems 3.1 and 3.3 is established by the following map:

σ : J ∈ ∂ × B H(x) → s ∈ {±1} p , where s i = +1 if i ∈ E = (x), J i, : = A i, : , -1 if i ∈ E = (x), J i, : = B i, : , (3.4a)
where we have used the definition (2.2) of p. The map is well defined since A i, : = B i, : when i ∈ E = (x). Furthermore, σ is bijective since two Jacobians in ∂ × B H(x) only differ by their rows with index in E = (x) and that these rows can take any of the values A i, : or B i, : . Actually, its reverse map is

σ -1 : s ∈ {±1} p → J ∈ ∂ × B H(x), where J i, : = A i, : if i ∈ E = (x), s i = +1, B i, : if i ∈ E = (x), s i = -1. (3.4b)
The question that arises is whether σ is also a bijection between ∂ B H(x) and S.

Proposition 3.4 (bijection ∂ B H(x) ↔ S) Let H : R n → R m be given by (1.3), x be a point in R n such that p = 0 and V be given by (3.1). Then, the map σ is a bijection from ∂ B H(x) onto S. In particular, the following properties hold.

1)

If J ∈ ∂ B H(x), then ∃ d ∈ R n such that σ(J) q (V T d) > 0. 2) If s ∈ {±1} p and ∃ d ∈ R n such that s q (V T d) > 0, then σ -1 (s) ∈ ∂ B H(x).
3)

Let J ∈ ∂ × B H(x). Then, J ∈ ∂ B H(x) ⇐⇒ σ(J) q (V T d) > 0 holds for some d ∈ R n .
Proof The properties 1, 2 and 3 in the statement of the proposition are straightforward consequences of the bijectivity of σ : ∂ B H(x) → S. Now, the discussion before the proposition has shown that σ : ∂ × B H(x) → {±1} p is a bijection. Therefore, σ : ∂ B H(x) → {±1} p is injective and it suffices to prove that

σ(∂ B H(x)) = S. (3.5a) [⊆ or point 1] Let J ∈ ∂ B H(x).
We have to show that s := σ(J) ∈ S, which means that one can find a d ∈ R n such that s q (V T d) > 0. By J ∈ ∂ B H(x), there exists a sequence {x k } ⊆ D H converging to x such that

H ′ (x k ) → J. (3.5b) For i ∈ E = (x), one cannot have (Ax k + a) i = (Bx k + b) i , since A i, : = B i, : would imply that x k / ∈ D H (lemma 2.1
). Therefore, one can find a subsequence K of indices k and a partition

(A 0 , B 0 ) of E = (x)
such that for all k ∈ K:

(Ax k + a) A0 < (Bx k + b) A0 and (Ax k + a) B0 > (Bx k + b) B0 . (3.5c) Now, fix k ∈ K and set d := x k -x. Since (Ax + a) i = (Bx + b) i for i ∈ E = (x), one deduces from (3.5c) that (B -A) A0, : d > 0 and (B -A) B0, : d < 0.
Recalling the definitions of V in (3.1) and S in (3.2), we see that, to conclude the proof of the membership σ(J) ∈ S, it suffices to show that [σ(J)] A0 = +1 and [σ(J)] B0 = -1 or, equivalently, by the definition of σ, (J i, : = A i, : for i ∈ A 0 ) and (J i, : = B i, : for i ∈ B 0 ). This is indeed the case, since by (3.5c), for all k ∈ K, one has (

H ′ i (x k ) = A i, : for i ∈ A 0 ) and (H ′ i (x k ) = B i, : for i ∈ B 0 );
now, use the convergence (3.5b) to conclude.

[⊇ or point 2] Let s ∈ S. We have to find a J ∈ ∂ B H(x) such that σ(J) = s, that is, which satisfies for i ∈ [1 : p]:

(J i, : = A i, : if s i = +1) and (J i, : = B i, : if s i = -1).

(3.5d) Since s ∈ S, there is a d ∈ R n such that s q (V T d) > 0.

(3.5e) Take a real sequence {t k } ↓ 0 and define the sequence {x k } ⊆ R n by

x k := x + t k d.
Then, x k → x. We claim that, for k sufficiently large, x k ∈ D H and H ′ (x k ) is a constant matrix J satisfying (3.5d), which will conclude the proof. Let i ∈ [1 : m]. The equivalence between problems 3.3 and 3.6 is obtained thanks to the following bijection

r If i ∈ A(x), (Ax k + a) i < (Bx k + b) i for k large, so that x k ∈ D H and H ′ i (x k ) = A i, : . r If i ∈ B(x), (Ax k + a) i > (Bx k + b) i for k large, so that x k ∈ D H and H ′ i (x k ) = B i, : . r If i ∈ E = (x), then A i, : = B i, : , so that x k ∈ D H and H ′ i (x k ) = A i, : = B i, : . r If i ∈ E = (x), subtract side by side (Ax k + a) i = (Ax + a) i + t k A i, : d and (Bx k + b) i = (Bx + b) i + t k B i, : d, use (Ax + a) i = (Bx + b) i and next (3.5e) to get (Bx k + b) i -(Ax k + b) i = t k (B i, : -A i, : )d = t k V T i, : d > 0 if s i = +1, < 0 if s i = -1. Hence, x k ∈ D H , (H ′ i (x k ) = A i, : if s i = +1) and (H ′ i (x k ) = B i, : if s i = -1
ı : s ∈ {±1} p → ı(s) := {i ∈ [1 : p] : s i = +1} ∈ P([1 : p]), (3.6)
whose reverse map is ı -1 :

I ∈ P([1 : p]) → s ∈ {±1} p , where s i = +1 if i ∈ I and s i = -1 if i / ∈ I.
As announced above, this equivalence relies on Gordan's theorem of the alternative [40; 1873]: for

a matrix A ∈ R m×n , ∃ x ∈ R n : Ax > 0 ⇐⇒ ∄ α ∈ R m + \ {0} : A T α = 0.
s ∈ S c ⇐⇒ ∄ x ∈ R n : Ax > 0 [definition of S in (3.2)] ⇐⇒ ∃ α ∈ R m + \ {0} : A T α = 0 [Gordan's alternative (3.7)] ⇐⇒ ∃ α ∈ R m + \ {0} : s q α ∈ N (V ) ⇐⇒ N (V ) ∩ O p I = {0} [see below] (3.8) ⇐⇒ I ∈ I [definition of I].
The implication "⇒" in (3.8) is due to the fact that s q α is nonzero and belongs to both N (V ) and O p I . The reverse implication "⇐" in (3.8) is due to the fact that there is a nonzero y ∈ N (V ) ∩ O p I , implying that α := s q y is nonzero and 0 and is such that s q α = y ∈ N (V ). Recall that the nullity of a matrix A, denoted by null(A), is the dimension of its null space. Let us introduce the following collection of index sets (from now on, J usually denotes a set of indices rather than a Jacobian matrix):

C := {J ⊆ [1 : p] : J = ∅, null(V : ,J ) = 1, V : ,J0 is injective if J 0 J}, (3.9)
where " " is used to denote strict inclusion. In the terminology of the vector matroid formed by the columns of V and its subsets made of linearly independent columns [55; proposition 1.1.1], the elements of C are called the circuits of the matroid [55; proposition 1.3.5(iii)]. The particular expression (3.9) of the circuit set is interesting in the present context, since it readily yields the following implication:

J ∈ C =⇒ any nonzero α ∈ N (V : ,J ) has none zero component. (3.10)
From (3.9) and (3.10), one can associate with J ∈ C a pair of sign vectors ±s ∈ {±1} J by s := sgn(α) for some nonzero α ∈ N (V : ,J ); the sign vectors ±s do not depend on the chosen α ∈ N (V : ,J ) \ {0} since null(V : ,J ) = 1. We call such a sign vector a stem vector, because of proposition 3.10 below, which shows that any s ∈ S c can be generated from such a stem vector.

Definition 3.9 (stem vector) A stem vector is a sign vector s = sgn(α), where α ∈ N (V : ,J ) for some J ∈ C. ⊓ ⊔

Note that there are twice as many stem vectors as circuits and that the stem vectors do not have all the same size. The matrix V in example 3.2 has J = [1 : 3] as single circuit. Since V e = 0, the associated stem vectors are ±e = ±(1, 1, 1). The next proposition now confirms that ±(1, 1, 1) are the only elements of S c . Proposition 3.10 (generating S c from the stem vectors) For s ∈ {±1} p , s ∈ S c ⇐⇒ s J = s for some J ⊆ [1 : p] and some stem vector s.

{d ∈ R n : s j v T j d > 0 for all j ∈ J } = ∅, (3.12a) ∀ J 0 J, {d ∈ R n : s j v T j d > 0 for all j ∈ J 0 } = ∅. (3.12b)
To determine such a J, start with J = [1 : p], which verifies (3.12a), since s ∈ S c . Next, remove an index j from [1 : p] if (3.12a) holds for J = [1 : p] \ {j}. Pursuing the elimination of indices j in this way, one arrives to an index set J satisfying (3.12a) and {d ∈ R n : s j v T j d > 0 for all j ∈ J \ {j 0 }} = ∅ for all j 0 ∈ J. Then, (3.12b) clearly holds. We claim that, for a J satisfying (3.12a) and (3.12b), s J is a stem vector, which will conclude the proof of the implication.

To stick to definition 3.9, we start by showing that J is a matroid circuit. By (3.12a), J = ∅. By Gordan's alternative (3.7), (3.12a) and (3.12b) read

∃ α ∈ R J + \ {0} such that j∈J s j v j α j = 0, (3.12c)
∀ J 0 J, ∄ α ′ ∈ R J0 + \ {0} such that j∈J0 s j v j α ′ j = 0. (3.12d)
From these properties, one deduces that α > 0 and that null(V : ,J ) 1. To show that null(V : ,J ) = 1, we proceed by contradiction. Suppose that there is a nonzero α ′′ ∈ R J that is not colinear with α and that verifies j∈J s j v j α ′′ j = 0. One can assume that t := max{α ′′ j /α j : j ∈ J} is > 0 (take -α ′′ otherwise). Set J 0 := {j ∈ J : α ′′ j /α j < t}. By the non-colinearity of α and α ′′ , on the one hand, and the definition of t, on the other hand, one has ∅ J 0 J. Furthermore,

α ′ := α -α ′′ /t 0, α ′ j > 0 for j ∈ J 0 and α ′ j = 0 for j ∈ J \ J 0 . Therefore, j∈J0 s j v j α ′ j = j∈J s j v j α ′ j = 0, yielding a contradiction with (3.12d).
To show that J ∈ C, we still have to prove that V : ,J0 is injective when J 0 J. Equivalently, it suffices to show that any β ∈ N (V : ,J ) with some zero component vanishes. We proceed by contradiction. If there is a β ∈ N (V : ,J ) \ {0} with a zero component, s J q α and β would be two linearly independent vectors in N (V : ,J ) (since s J q α has no zero component), contradicting null(V : ,J ) = 1. Now, since s J = sgn(s J q α), since s J q α ∈ N (V : ,J ) by (3.12c) and since J is a matroid circuit of V , s J is a stem vector.

[⇐] Since s J is a stem vector, it follows that s J := sgn(α) for some α ∈ R J with nonzero components that satisfies V : ,J α = 0. Then, there is no d ∈ R n such that s J q (V T : ,J d) > 0 (otherwise, (s J q α q s J ) q (V T : ,J d) > 0, because s J q α > 0, or α q (V T : ,J d) > 0, implying that 0 = α T (V T : ,J d) > 0, a contradiction). Hence, there exists certainly no d ∈ R n such that s q (V T d) > 0. This implies that s ∈ S c .

⊓ ⊔

To determine the stem vectors, which are based on the matroid circuits of V defined by (3.9), one has to select subsets of columns of V forming a rank one matrix, whose strict subsets form injective matrices. Actually, this last condition can be simplified by the following property. Proposition 3.11 (matroid circuit detection) Suppose that I ⊆ [1 : p] is such that null(V : ,I ) = 1 and that α ∈ N (V : ,I ) \ {0}. Then, J := {i ∈ I : α i = 0} is a matroid circuit of V and the unique one included in I.

Proof 1) Let us show that J is a matroid circuit.

Since α = 0, one has J = ∅. Let us show that null(V : ,J ) = 1. Since J ⊆ I, one has null(V : ,J ) null(V : ,I ) = 1. Furthermore, α J ∈ N (V : ,J ) \ {0} implies that null(V : ,J ) 1.
Now, let J 0 J and suppose that V : ,J0 β = 0. We have to show that β = 0. Since V : ,J (β, 0 J\J0 ) = 0, it follows that (β, 0 J\J0 ) ∈ N (V : ,J ), which is of dimension 1, so that (β, 0 J\J0 ) is colinear to α. Since the components of α are = 0, we get that β = 0.

2) Let us now show that J is the unique matroid circuit of V included in I.

Let J ′ be a matroid circuit of V included in I. Then null(V : ,J ′ ) = 1 and there is a nonzero α ′ ∈ N (V : ,J ′ ). By (3.10), α ′ has nonzero components. Furthermore, (α ′ , 0 I\J ′ ) ∈ N (V : ,I ), which has unit dimension and contains α. Therefore, α and (α ′ , 0 I\J ′ ) are colinear. Since the components of α are = 0, we get that J ′ = J.

⊓ ⊔

Convex analysis problems

The formulation of the original problem 3.1 in the form of the convex analysis problems 3.12 and 3.15 below may be useful to highlight some properties of ∂ B H(x), thanks to the tools of that discipline.

Pointed cones by vector inversions

Recall that a convex cone K of R n is a convex set verifying R ++ K ⊆ K (or, more explicitly, tx ∈ K when t > 0 and x ∈ K). A closed convex cone K is said to be pointed if K ∩ (-K) = {0} [16; p. 54],
which amounts to saying that K does not contain a line (i.e., an affine subspace of dimension one)

or that K has no nonzero direction z such that -z ∈ K. For P ⊆ R n , we also denote by "cone P " the smallest convex cone containing P . Problem 3.12 (pointed cones by vector inversions) Let be given two positive integers n and

p ∈ N * and p vectors v 1 , . . . , vp ∈ R n \ {0}. It is requested to determine all the sign vectors s ∈ {±1} p such that cone{s i v i : i ∈ [1 : p]} is pointed.

⊓ ⊔

The equivalence between the original problem 3.1 and problem 3.12 is obtained thanks to the next proposition, which gives another property ("cone pointedness") that is equivalent to those in (3.7) and that is adapted to the present concern. For a proof, see [41; theorem 2.3.29] or [29]. Proposition 3.13 (pointed polyhedral cone) For a finite collection of nonzero vectors {w i : i ∈ [1 : p]} ⊆ R n , the following properties are equivalent:

(i

) cone{w i : i ∈ [1 : p]} is pointed, (ii) ∄ α ∈ R p + \ {0} : i∈[1 : p] α i w i = 0, (iii) ∃ d ∈ R n , ∀ i ∈ [1 : p] : w T i d > 0.
Equivalence 3.14 (signed linear system feasibility ↔ pointed cone by vector inversion)

The equivalence (i) ⇔ (iii) of the previous proposition shows that the set S defined by (3.2) is also given by

S = {s ∈ {±1} p : cone{s i v i : i ∈ [1 : p]} is pointed}. (3.13)
To put it in words, denoting by v 1 , . . . , vp the columns of the matrix V defined by (3.1), the signed feasibility problem 3.3 is equivalent to problem 3.12. ⊓ ⊔

Linearly separable bipartitions of a finite set

This section extends section 3.3.1 and adopts its concepts and notation. The point of view presented in this section was also shortly considered by Zaslavsky [79;1975, § 6A]. This enumeration problem appears in the study of neural networks [76]. Baldi and Vershynin [7] make the connection with homogeneous linear threshold functions and highlight its impact in deep learning [START_REF] Schmidhuber | Deep learning in neural networks: An overview[END_REF][START_REF] Baldi | Deep learning in biomedical data science[END_REF].

v1 v2 -v3 d v1 -v2 -v3 d v1 v3 -v2 d v3 -v1 -v2 d v2 v3 -v1 d v2 -v1 -v3 d Fig. 3.1
The figure is related to the linear complementarity problem defined by example 3.2: the v i 's are the columns of the matrix V (their third zero components are not represented). Each of the 6 sets of vectors plots the 3 vectors {s i v i : i ∈ [1 : 3]}, for each of the 6 sign vectors s ∈ S (given by the columns of the matrix S in (3.3)), as well as a direction d (given by the columns of D in (3.3), dashed lines) such that

s i v T i d > 0 for all i ∈ [1 : 3]. Each conic hull of these vectors, namely cone{s i v i : i ∈ [1 : 3]}, is pointed. The conic hulls of {v 1 , v 2 , v 3 } and {-v 1 , -v 2 , -v 3 }
are both the space of dimension 2, hence there are not pointed, which confirms the fact that (1, 1, 1) and (-1, -1, -1) are not in S.

Problem 3.15 (linearly separable bipartitioning) Let be given an affine space A and p ∈ N * vectors v1 , . . . , vp ∈ A. Let A 0 := A -A be the vector space parallel to A, endowed with a scalar product •, • . It is requested to find all the ordered bipartitions (i.e., the partitions made of two subsets) (I, J) of [1 : p] for which there exists a vector ξ ∈ A 0 (also called separating covector below) such that

∀ i ∈ I, ∀ j ∈ J : ξ, vi < ξ, vj . ⊓ ⊔
Of course, if (I, J) is an appropriate ordered bipartition to which a separating covector ξ corresponds, then (J, I) is also an appropriate ordered bipartition with separating covector -ξ. Therefore, only half of the appropriate ordered bipartitions (I, J) must be identified, a fact that is related to the symmetry of ∂ B H(x) (proposition 4.1). Figure 3. ing the separating hyperplanes {v ∈ A : ξ T v = t} corresponding to some separating covector ξ and some t ∈ R, for three examples with p = 4. Since it will be shown that |S| is the number of these searched linearly separable bipartitions, this one is denoted that way in the figure. Obviously, |S| not only depends on p and r := dim(vect{v 1 , . . . , vp}) + 1, but it also depends on the arrangement of the vi 's in the affine space A. We also see that |S| cannot take all the even values (proposition 4.1) between its lower bound 2p = 8 and its upper bounds 8 (if r = 2) and 14 (if r = 3) given by propositions 4.7 and 4.10.

The equivalence between the linearly separable bipartitioning problem 3.15 of this section and the vector inversion problem 3.12 (hence, with the original problem 3.1) is grounded on the following construction and proposition. Construction 3.16 1) Let be given two integers n and p ∈ N * and p nonzero vectors v 1 , . . . ,

vp ∈ R n such that K := cone{v k : k ∈ [1 : p]} is a pointed cone. From proposition 3.13, there is a direction d ∈ R n such that d = 1 and ∀ k ∈ [1 : p] : v T k d > 0 . Define A := {v ∈ R n : d T v = 1}, A 0 := A -A = {v ∈ R n : d T v = 0}, ∀ k ∈ [1 : p] : vk := v k /(v T k d) ∈ A.
2) For a given bipartition (I, J) of [1 : p], define

K I := cone{v i : i ∈ I}
and

K J := cone{v j : j ∈ J}, (3.14a)
C I := K I ∩ A
and

C J := K J ∩ A, (3.14b)
with the convention K ∅ = {0} and C ∅ = ∅.

⊓ ⊔ Proposition 3.17 (pointed cone after vector inversions) Adopt the construction 3.16 and take a partition (I, J) of [1 : p]. Then, the following properties are equivalent:

(i) cone((-K I ) ∪ K J ) is pointed, (ii) K I ∩ K J = {0}, (iii) C I ∩ C J = ∅, (iv) there exists a vector ξ ∈ A 0 such that max i∈I ξ T vi < min j∈J ξ T vj . Proof [(i) ⇒ (ii)] We show the contrapositive. If there is v ∈ (K I ∩ K J ) \ {0}, then -v ∈ (-K I ) ⊆ cone((-K I ) ∪ K J ) and v ∈ K J ⊆ cone((-K I ) ∪ K J ). Therefore, cone((-K I ) ∪ K J ) is not pointed. [(ii) ⇒ (iii)] ∅ = A ∩ {0} = A ∩ K I ∩ K J [(ii)] = (A ∩ K I ) ∩ (A ∩ K J ) = C I ∩ C J .
[(iii) ⇒ (iv)] We claim that C I is nonempty, convex and compact.

Indeed, since C I is nonempty (it contains the vectors vi for i ∈ I = ∅), convex (because K I and A are convex) and closed (because K I and A are closed), it suffices to show that C I is bounded or that its asymptotic cone (or recession cone in [66; p. 61]), namely theorem 8.4]. This is indeed the case since v T d > 0 for all v ∈ K I \ {0}. For the same reason, C J is nonempty, convex and compact. Now, since C I ∩ C J = ∅ by (iii), one can strictly separate the convex sets C I and C J in A [66; corollary 11.4.2]: there exists ξ ∈ A 0 such that ξ T v < ξ T w, for all v ∈ C I and all w ∈ C J . This shows that (iv) holds.

C ∞ I = K I ∩ A 0 , is reduced to {0} [66;
[(iv)

⇒ (i)] Since cone((-K I ) ∪ K J ) = cone({-v i : i ∈ I} ∪ {v j : j ∈ J}), by proposition 3.13, it suffices to find d (I,J) ∈ R n such that -v T i d (I,J) > 0, ∀ i ∈ I and v T j d (I,J) > 0, ∀ j ∈ J . (3.15)
By (iv) and the fact that θ ∈ (0,

π) → cot θ ∈ R is surjective, one can determine θ ∈ (0, π) such that max i∈I ξ T v i v T i d < -cot θ < min j∈J ξ T v j v T j d . (3.16) Since sin θ > 0 for θ ∈ (0, π) and since v T k d > 0 for all k ∈ [1 : p], this is equivalent to max i∈I v T i [(cos θ)d + (sin θ)ξ] < 0 < min j∈J v T j [(cos θ)d + (sin θ)ξ].
Therefore, (3.15) is satisfied with d (I,J) := (cos θ)d + (sin θ)ξ.

⊓ ⊔

One can now establish the link between the pointed cone problem of section 3.3.1 (problem 3.12) and the linearly separable bipartitioning problem (problem 3.15). Equivalence 3.18 (pointed cone ↔ linearly separable bipartitioning) Let be given a matrix V ∈ R n×p with nonzero columns denoted by v 1 , . . . , vp and take s ∈ S, which is nonempty. By (3.13)

, cone{s i v i : i ∈ [1 : p]} is pointed. Use the construction 3.16(1) with v i s i v i .
For s ∈ {±1} p , define a partition (I, J) of [1 : p] by

I := {i ∈ [1 : p] : si s i = -1} and J := {i ∈ [1 : p] : si s i = +1}.
Define also K I and K J by (3.14a) with v i s i v i . We claim that

cone{s i v i : i ∈ [1 : p]} is pointed ⇐⇒ ∃ ξ ∈ A 0 : max i∈I ξ T vi < min j∈J ξ T vj . (3.17)
Indeed, one has

cone{s i v i : i ∈ [1 : p]} is pointed ⇐⇒ cone{s i s i (s i v i ) : i ∈ [1 : p]} is pointed ⇐⇒ cone((-K I ) ∪ K J ) is pointed ⇐⇒ ∃ ξ ∈ A 0 : max i∈I ξ T vi < min j∈J ξ T vj ,
where we have used the equivalence (i) ⇔ (iv) of proposition 3.17 (v i s i v i ). The equivalence (3.17) establishes the expected equivalence between the pointed cone problem 3.12 (in which one looks for all the s ∈ {±1} p such that cone{s i v i : i ∈ [1 : p]} is pointed) and the linearly separable bipartitioning problem 3.15 of the vectors vi =

s i v i /(s i v T i d) = v i /(v T i d), i ∈ [1 : p],
where d is associated with the pointed cone cone{s i v i : i ∈ [1 : p]} by the equivalence (i) ⇔ (iii) of proposition 3.13.

⊓ ⊔

Discrete geometry: hyperplane arrangements

The equivalent problem examined in this section has a long history, going back at least to the XIXth century [START_REF] Steiner | Einige Gesetze über die Theilung der Ebene und des Raumes[END_REF][START_REF] Roberts | On the figures formed by the intercepts of a system of straight lines in a plane, and on analogous relations in space of three dimensions[END_REF]. More recently, it appears in Computational Discrete Geometry (the discipline has many other names), under the name of hyperplane arrangements. Contributions to this problem, or a more general version of it, with a discrete mathematics point of view, have been reviewed in [START_REF] Grünbaum | Arrangements and Spreads[END_REF][START_REF] Edelsbrunner | Algorithms in Combinatorial Geometry[END_REF][START_REF] Stanley | An introduction to hyperplane arrangements[END_REF]2,[START_REF] Halperin | Arrangements[END_REF]. It has many applications [START_REF] Edelsbrunner | Constructing arrangements of lines and hyperplanes with applications[END_REF][START_REF] Helena | Hyperplane arrangements -Construction, visualization and application[END_REF]18]. From an algorithmic point of view, the algorithms developed in this domain can immediately be used to compute S defined by (

H i := {d ∈ R n : v T i d = 0}.
H + i := {d ∈ R n : v T i d > 0}
and

H - i := {d ∈ R n : v T i d < 0}.
The problem is to determine the following set of open sectors or cells of R n , indexed by the bipartitions (I + , I -) of [1 : p]:

C := (I + , I -) ∈ B([1 : p]) : (∩ i∈I+ H + i ) ∩ (∩ i∈I-H - i ) = ∅ , (3.19)
where B([

1 : p]) denotes the set of bipartitions of [1 : p]. ⊓ ⊔ v1 v2 v3 H1 H2 H3 (+--) (-++) (-+-) (+-+) (--+) (++-)
Fig. 3.3 Illustration of problem 3.19 (arrangement of hyperplanes containing the origin) for the 3 vectors that are the columns on the matrix V in example 3.2 (since the last components of these v i 's vanish, only the first two ones are represented above). The hyperplanes H i are defined by (3.18). The regions to determine are represented by the sign vectors here denoted (s 1 s 2 s 3 ) with s i = ±: if d ∈ R 2 belongs to the region (s 1 s 2 s 3 ),

then s i = + if v T i d > 0 and s i = -if v T i d < 0.
We see that there are only 6 = 2p regions among the 8 = 2 p possible ones; the regions (+++) and (---) are missing, which reflects the fact that

+ v 1 + v 2 + v 3 = 0 and -v 1 -v 2 -v 3 = 0 (see problem 3.6).
The link between problem 3.19 and the signed feasibility of strict linear inequality systems of section 3.2.1 is obtained from the bijection

η : (I + , I -) ∈ B([1 : p]) → s ∈ {±1} p , where s i = +1 if i ∈ I + , -1 if i ∈ I - (3.20)
and the setting V = v 1 • • • vp , whose columns are nonzero by assumption, here and in section 3.2.1.

Recall the definition (3.2) of the set of sign vectors S.

Proposition 3.20 (bijection C ↔ S) For the matrix V ∈ R n×p , with nonzero columns v i 's, the map η given by (3.20) is a bijection from C onto S.

Proof Let (I + , I -) ∈ B([1 : p]) and s := η((I + , I -)). Then, Besides their theoretical relevance, the properties of the B-differential of H given in this section will also be useful to design the algorithms presented in section 5 and to check the correctness of their implementation.

(I + , I -) ∈ C ⇐⇒ ∃ d ∈ (∩ i∈I+ H + i ) ∩ (∩ i∈I-H - i ) ⇐⇒ ∃ d ∈ R n : (v T i d > 0 for i ∈ I + ) and (v T i d < 0 for i ∈ I -) ⇐⇒ ∃ d ∈ R n : s q (V T d) > 0 ⇐⇒ s ∈ S.

Some properties of the B-differential

Let us start with a basic property of ∂ B H(x), which is its symmetry in the sense of definitions 2.5. This property has been observed by many in other contexts [2; § 1.1.4], so that we leave its short proof, based on the equivalence 3.5, to [29]. It is useful for the algorithms since it implies that only half of the B-differential has to be computed. We now give a necessary and sufficient condition ensuring the completeness of ∂ B H(x) in the sense of definition 2.4. The condition was shown to be sufficient in [78; corollary 2.1(i)] for the nonlinear case (1.6), using a different proof, but we shall see in [START_REF] Dussault | Partial description of the B-differential of the componentwise minimum of two vector functions by linearization[END_REF] that it is an easy consequence of that property in the affine case (1.3). Thanks to the equivalence 3.5, the present proof is short. This property is also useful in the development of algorithms, as a test that these must pass:

|∂ B H(x)| = 2 p if and only if V ∈ R n×p is injective. Proposition 4.2 (completeness of the B-differential) The B-differential ∂ B H(x) of H at x
is complete if and only if the matrix V ∈ R n×p in (3.1) is injective. Hence, this property can hold only if p n.

Proof [⇒]

We show the contrapositive. Assume that V is not injective, so that V α = 0 for some nonzero α ∈ R p . With s ∈ sgn(α), one can write

i∈[1 : p] |α i |s i v i = 0.
By Gordan's alternative (3.7), it follows that there is no d ∈ R n such that s q (V T d) > 0. By (3.2), this implies that s / ∈ S. According to the equivalence 3.5, σ -1 (s) / ∈ ∂ B H(x), showing that the B-differential is not complete.

[⇐] Assume the injectivity of V . Let s ∈ {±1} p . Since V T is surjective, the system V T d = s holds for some d ∈ R n . For this d, s q (V T d) = e, so that s q (V T d) > 0 holds for some d ∈ R n , which implies that the selected s is in S. We have shown that

S = {±1} p or that ∂ B H(x) = σ -1 ({±1} p ) (σ -1 is defined by (3.4b)) is complete. ⊓ ⊔
We focus now on the connectivity of ∂ B H(x), a notion that is more easily presented in terms of S ⊆ {±1} p but that can be transferred straightforwardly to ∂ B H(x) by the bijection σ defined in (3.4). This property was implicitly used, for instance, in the algorithms proposed by Avis, Fukuda and Sleumer [START_REF] Avis | Reverse search for enumeration[END_REF][START_REF] Helena | Output-sensitive cell enumeration in hyperplane arrangements[END_REF] for hyperplane arrangements. Proof [⇒] We prove the contrapositive. Suppose that the columns v i and v j of V are colinear:

v j = αv i , for some α ∈ R * . Assume that α > 0 (resp. α < 0). By (3.2), for any s ∈ S = ∅, one can find d ∈ R n such that s q (V T d) > 0, implying that s i = s j (resp. s i = -s j ). Therefore, one cannot find a path in S joining s ∈ S and -s ∈ S (proposition 4.1), since one would have to change the two components with index in {i, j} and that these components must be changed simultaneously for the sign vectors in S, while the adjacency property along a path prevents from changing more than one sign at a time.

[⇐] We leave to [29] the proof of this implication and of the last claim of the proposition, since the conclusion of the implication is given in [2; section 1.10.4] as a simple observation with a very different point of view, related to graph theory.

⊓ ⊔

For k ∈ [1 : p], we introduce

S k := {s ∈ {±1} k : ∃ d ∈ R n such that s i v T i d > 0 for i ∈ [1 : k]}. (4.3)
We also note S c k := {±1} k \ S k . Hence S = Sp and S c = S c p . Point 1 of the next proposition will be used to motivate an improvement of algorithm 5.5 in section 5.2.4 and its points 2 and 3 will be used to get the equivalence in proposition 4.13, related to a fan arrangement. 3) If v k+1 is not colinear to any of the vectors v 1 , . . . , v k , then, [(s, ±1) and (-s, ±1) ∈ S k+1 for one s ∈ S k ] and [(s ′ , +1) or (s ′ , -1) ∈ S k+1 for any

s ′ ∈ S k ]. In particular, |S k+1 | |S k | + 2. Proof 1) If s ∈ S c k , there is no d ∈ R n such that s i v T i d > 0 for i ∈ [1 : k]. Therefore, there is no d ∈ R n such that (s i v T i d > 0 for i ∈ [1 : k]) and ±v T k+1 d > 0. Therefore, (s, ±1) ∈ S c k+1 . This implies that |S c k+1 | 2|S c k |.
2) Let P be the orthogonal projector on vect{v 1 , . . . , v k } ⊥ for the Euclidean scalar product. By assumption, P v k+1 = 0. Let s ∈ S k , so that there is a direction

d ∈ R n such that s i v T i d > 0 for i ∈ [1 : k]. For any t ∈ R and i ∈ [1 : k], the directions d ± := d ± t P v k+1 verify s i v T i d ± = s i v T i d > 0 (because v T i P v k+1 = 0). In addition, for t > 0 sufficiently large, one has ±v T k+1 d ± = ±v T k+1 d + t P v k+1
2 > 0 (because P 2 = P and P T = P). We have shown that both (s, +1) and (s, -1) are in S k+1 . Therefore,

|S k+1 | 2|S k |. Now, |S k | + |S c k | = 2 k , |S k+1 | + |S c k+1 | = 2 k+1 and |S c k+1 | 2|S c k | by point 1. Therefore, one must have |S k+1 | = 2|S k | and |S c k+1 | = 2|S c k |.
3) We claim that one can find a direction d ∈ R n such that

∀ i ∈ [1 : k] : v T i d = 0 and v T k+1 d = 0. (4.4)
Indeed, let E := {d ∈ R n : v T k+1 d = 0} and P be the orthogonal projector on E for the Euclidean scalar product. By lemma 2.6, one can find a direction

d ∈ E (hence v T k+1 d = 0) such that |{(P v i ) T d : i ∈ [1 : k + 1]}| = |{P v i : i ∈ [1 : k + 1]}|. Since P v k+1 = 0 and P v i = 0 for i ∈ [1 : k] (because the v i 's are not colinear with v k+1 ), one has (P v i ) T d = 0 for i ∈ [1 : k]. Since, 0 = (P v i ) T d = v T i P d = v T i d, (4.4) follows.
Taking

s i := sgn(v T i d) for i ∈ [1 : k],
one deduces from (4.4) that there is a direction d ∈ R n such that

∀ i ∈ [1 : k] : s i v T i d > 0 and v T k+1 d = 0.
It follows that, for ε > 0 sufficiently small, the directions d ± := d ± εv k+1 satisfy

∀ i ∈ [1 : k] : s i v T i d ± > 0 and ± v T k+1 d ± > 0.
This means that (s, ±1) ∈ S k+1 . By symmetry (proposition 4.1), one also has (-s, ±1) ∈ S k+1 , so that we have found 4 vectors in S k+1 . Now, since, for any

s ′ ∈ S k \ {±s} (in number |S k | -2), either (s ′ , +1) ∈ S k+1 or (s ′ , -1) ∈ S k+1 , it follows that |S k+1 | 4 + (|S k | -2) = |S k | + 2. ⊓ ⊔

Cardinality of the B-differential

Information on the cardinality of ∂ B H(x) can be useful to check the correctness of the number of elements computed by the algorithms presented in section 5.2.

Winder's formula

Giving the exact number of elements in [77; 1966] and reads for the matrix V with nonzero columns given by (3.1)

∂ B H(x), that is |∂ B H(x)| = |S| = |C| = 2 p -|S c | = 2 p -
|∂ B H(x)| = I⊆[1 : p] (-1) null(V : ,I ) , (4.5)
where null(V : ,I ) is the nullity of V : ,I and the term in the right-hand side corresponding to I = ∅ is 1 (one takes the convention that null(V : ,∅ ) = 0). Note that, in this formula, the columns of V can be colinear with each other. This amazing expression, with its only algebraic nature, potentially made of positive and negative terms, is explicit but, to our knowledge, has not been at the origin of a method to list the elements of ∂ B H(x). We give in [29] a proof of (4.5) that follows the same line of reasoning as the one of Winder [START_REF] Winder | Partitions of N-space by hyperplanes[END_REF], but that is more analytic in that it uses the sign vectors introduced in section 3.2.1 rather than geometric arguments (i.e., the hyperplane arrangements of section 3.4).

Bounds

When p is large, computing the cardinality |∂ B H(x)| from (4.5) by evaluating the 2 p ranks rank(V : ,I ) for I ⊆ The upper bound was already mentioned in proposition 2.2. ⊓ ⊔ Proposition 4.10 below provides a refinement of the upper bound given by proposition 4.7. The next proposition will be useful for this purpose. Recall that a function ϕ : x ∈ T → ϕ(x) ∈ R, defined on a topological space T, is said to be lower semicontinuous if, for any x ∈ T and any ε > 0, there is a neighborhood V of x such that, for all x ∈ V, one has ϕ(x) ϕ(x) + ε. It is known that the rank of a matrix can only increase in the neighborhood of a given matrix, which implies its lower semicontinuity. The next lemma shows that the same property holds for |S| ∈ N * , viewed as a function of V . Recall that the bijection σ is defined by (3.4).

is viewed as a function of V ∈ R n×p . Then, S(V ) ⊆ S( Ṽ ) for Ṽ near V in R n×p . In particular, V ∈ R n×p → |S(V )| ∈ N * is lower semicontinuous.
Proof By the definition (3.2) of S(V ), for all s ∈ S(V ), there is a ds ∈ R n such that s q (V T ds) > 0.

Clearly, one still has s q ( Ṽ T ds) > 0, for Ṽ near V . Since S(V ) is finite, there is a neighborhood V of V , such that, for Ṽ ∈ V and s ∈ S(V ), there is a d ∈ R n such that s q ( Ṽ T d) > 0 or s ∈ S( Ṽ ).

We have shown that S(V ) ⊆ S( Ṽ ) for Ṽ near V . As a direct consequence of this inclusion, we have that where r := rank(V ).

|S(V )| |S( Ṽ )| for Ṽ near V . The lower semicontinuity of V → |S(V )| follows.

⊓ ⊔

This notion is used by Winder [START_REF] Winder | Partitions of N-space by hyperplanes[END_REF] when r = n. Example of vectors in general position are those in the left-hand side and right-hand side panes in figure 3.2 (the points are the normalized vectors vi 's so that the v i 's are actually in R 3 ); note that in the first case 2 = r < n = 3. Those in the middle pane are not in general position. This is due to the fact that r := rank(V ) = 3 while for the 3 bottom vectors, with indices in I say, one has min(|I|, r)rank(V : ,I ) = 3 -2 = 0.

Equality in the upper estimate (4.8) of the next proposition was shown by Winder [77;1966, corollary] when the columns of V are in general position and r = n, thanks to the identity (4.5).

Long before him, the Swiss mathematician Ludwig Schläfli [67; p. 211] established the identity under the same assumptions, before 1852 [67; p. 174], without reference to (4.5), which was probably not known at that time. Note that equality does not hold in (4.8) for the middle configuration in 32 )] = 14 (we have seen that the vectors in this pane are not in general position). The bound (4.8) is also useful to check the behavior of the algorithms for test-cases in which the columns of V are in general position. This is likely to be so for randomly generated V , and it was verified by all our random test-cases in section 5.2.6(B.1).

Proposition 4.10 (upper bound on |∂ B H(x)|) For V given by (3.1) and r := rank(V ), one has

|∂ B H(x)| 2 i∈[0 : r-1] p -1 i , (4.8)
with equality if and only if (4.7) holds.

Proof 1) The proof of the implication "(4.7) ⇒ (4.8) with equality" is established in [77; corollary], using the identity (4.5). See also [29].

2) Let us now show that (4.8) holds. Below, we systematically identify ∂ B H(x) and S, thanks to the equivalence 3.5. We also note S ≡ S(V ) to stress the dependence of S on V . Let β be the right-hand side of (4.8). We proceed by contradiction, assuming that there is a matrix V ∈ R n×p of rank r such that

|S(V )| > β.
|S( Ṽ )| = β < |S(V )|.
To find Ṽ of rank r arbitrarily close to V verifying (4.9b), we proceed as follows. Since (4.9a) holds, the first part of the proof implies that V does not satify (4.7). Our goal is to construct from V a matrix Ṽ of rank r arbitrarily close to V with columns in general position. Then, Ṽ satisfies (4.9b) by the first part of the proof.

In view of (4.6) and since V does not satisfy (4.7), there is some I ⊆ [1 : p] such that rank(V : ,I ) < min(|I|, r). By linear algebra arguments (see [29] for more details), one can get an arbitrarily small perturbation Ṽ : ,I of V : ,I , such that rank( Ṽ : ,I ) = min(|I|, r) and R( Ṽ : ,I ) ⊆ R(V ). Next, one forms Ṽ ∈ R n×p by setting Ṽ : ,I c = Ṽ : ,I c , so that Ṽ is as close to V as desired and verifies R( Ṽ ) ⊆ R(V ). The perturbation Ṽ : ,I of V : ,I can also perturb V : ,I ′ for other index sets I ′ ⊆ [1 : p]. However, one has rank( Ṽ : ,I ′ ) min(|I ′ |, r) by (4.6). Now, by the property of the rank, which can only increase in a neighborhood of a given matrix, if the perturbation taken above is sufficiently small, one has rank(V : ,I ′ ) rank( Ṽ : ,I ′ ) min(|I ′ |, r) for any I ′ ⊆ [1 : p]. Therefore, rank(V : ,I ′ ) = min(|I ′ |, r) implies that rank( Ṽ : ,I ′ ) = min(|I ′ |, r). As a result, the modification of V into Ṽ described above increases by at least one the number of intervals

I ′ ⊆ [1 : p] such that rank( Ṽ : ,I ′ ) = min(|I ′ |, r).
Since the number of such intervals is finite, proceeding similarly with all the nonempty index sets 3) One still has to show that "(4.8) with equality ⇒ (4.7)". We proceed by contradiction, assuming that (4.8) holds with equality for ∂ B H(x) ≡ S(V ), but that (4.7) does not hold. By (4.6), there exists I ⊆ [1 : p] such that rank(V : ,I ) < min(|I|, r).

I ′′ ⊆ [1 : p] such that rank( Ṽ : ,I ′′ ) < min(|I ′′ |,

(4.9c)

Let β = |S(V )| be the right-hand side of (4.8). It certainly suffices to show that, thanks to (4.9c), one can find a matrix Ṽ ∈ R n×p such that rank( Ṽ ) r and |S( Ṽ )| > β, since this would be in contradiction with what has been shown in part 2 of the proof. This matrix Ṽ is obtained by perturbing V . By proposition 4.8, if the perturbation is sufficiently small, one has S(V ) ⊆ S( Ṽ ), so that it suffices to show that S( Ṽ ) contains a sign vector s that is not in S(V ).

We claim that (4.9c) implies that one can find an index set J ⊆ I such that V : ,J is not injective and |J| r. (4.9d) Indeed, if |I| r, one can take J = I to satisfy (4.9d), since rank(V : ,I ) < |I| by (4.9c), so that V : ,I is not injective. If |I| > r, then rank(V : ,I ) < r by (4.9c), which implies that any J ⊆ I such that |J| = r satisfies (4.9d). Since V : ,J is not injective, one can find

α J ∈ R J \ {0} such that 0 = j∈J α j v j = j∈J sj |α j |v j ,
for some sJ ∈ {±1} J satisfying sj ∈ sgn(α j ) for all j ∈ J. Then, by Gordan's alternative (3.7),

∄ d ∈ R n : sj v j d > 0, for all j ∈ J.
This implies that there is no s ∈ S(V ) such that s J = sJ . To conclude the proof, it suffices now to show that one can construct an arbitrarily small perturbation Ṽ of V , such that R( Ṽ ) ⊆ R(V ) and with an s ∈ S( Ṽ ) satisfying s J = sJ . Let J c := [1 : p] \ J. By (4.9d), |J| r n so that one can find vectors {ṽ j : j ∈ [1 : p]}, such that ṽj = v j for j ∈ J c , the vectors {ṽ j : j ∈ J} are linearly independent, ṽj -v j is arbitrarily small and {ṽ j : j ∈ [1 : p]} ⊆ R(V ). Since the vectors {ṽ j : j ∈ J} are linearly independent, one can find a direction d 0 ∈ R n such that ṽT j d 0 = sj for j ∈ J, hence sj ṽT j d 0 > 0, ∀ j ∈ J.

(4.9e)

Set sj = 1 for j ∈ J c . Let d be a discriminating covector given by lemma 2.6 (there denoted ξ) for the vectors {0} ∪ {s i v i : i ∈ [1 : p]} sufficiently close to d 0 . It results that sj ṽT j d > 0 for j ∈ J (by (4.9e)) and that sj ṽT j d = 0 for j ∈ J c . Finally, we see that the sign vector s ∈ {±1} p defined by s i = sgn(ṽ T i d) for all i ∈ [1 : p] is in S( Ṽ ) and satisfies s J = sJ , as desired.

⊓ ⊔ Corollary 4.11 (stability of the sign vector set) The sign vector set S ⊆ {±1} p defined by

(3.
2) is unchanged by small variations of the matrix V ∈ R n×p preserving its rank, provided the columns v 1 , . . . , vp ∈ R n of V are in general position in the sense of definition 4.9. ⊓ ⊔

Proof If Ṽ is near V , S(V ) ⊆ S( Ṽ )

Particular configurations

We consider in this section some particular matrices V ∈ R n×p given by (3.1), which may be useful to get familiar with the B-differential of H. For these V 's, |∂ B H(x)| can be computed easily. We consider two matrices V with the property that r := rank(V )

rank(V : ,I ) = |I| if |I| 2 2 if |I| > 2.
Therefore (4.7) holds. By proposition 4.10, this implies that equality holds in (4.8), that is, with

r := rank(V ) = 2: |∂ B H(x)| = 2 i∈[0 : 1] ( p-1 i ) = 2p. [⇐] If |∂ B H(x)| =
2p, proposition 4.7 yields 2p max(2p, 2 r ) 2 r + 2(pr) 2p, so that equality holds in these inequalities. By the last one, 2 r = 2r, which only occurs for r ∈ {1, 2}. Since p 2 and the vectors are not colinear, one has r = 2.

⊓ ⊔

A glance at the C-differential

The section presents two links between the B-differential and the C-differential of the function H given by (1.3). The first proposition tells us that, whilst ∂ C H(x) can be obtained from ∂ B H(x) by taking its convex hull (it is its definition (1.2)), the latter can be obtained from the former by taking its extreme points. For a proof, see [29].

Proposition 4.14 (a link with the C-differential) ∂ B H(x) = ext ∂ C H(x).
The second proposition restates theorem 2.2 of Xiang and Chen [78; 2011], which applies to the more general nonlinear function (1.6). The interest of this restatement comes from its proof that is short, thanks to the use of the symmetry of the B-differential (proposition 4.1), and from the fact that proposition 4.15 can be used, straightforwardly, to recover Xiang and Chen's central C-Jacobian of H, given by (1.6); see [START_REF] Dussault | Partial description of the B-differential of the componentwise minimum of two vector functions by linearization[END_REF]. Recall the notation (2.1) of the index sets.

J i, : =    A i, : if i ∈ A(x), 1 2 [A i, : + B i, : ] if i ∈ E (x), B i, : if i ∈ B(x).
(4.10)

Proof Let M ∈ ∂ B H(x)
, which is known to be nonempty. By proposition 2.2, M i, : = A i, : for i ∈ A(x), M i, : = B i, : for i ∈ B(x) and M i, : = A i, : = B i, : for i ∈ E = (x). By the symmetry of

∂ B H(x) (proposition 4.1), M ′ defined by M ′ : ,i = M : ,i if i ∈ A(x) ∪ E = (x) ∪ B(x)
and by

M ′ i, : = B i, : if i ∈ E = (x) and M i, : = A i, : A i, : if i ∈ E = (x) and M i, : = B i, : is also in ∂ B H(x). Therefore, J = (M + M ′ )/2 is in co ∂ B H(x) = ∂ C H(x), by (1.
2). This is the formula of J given in the statement of the proposition.

⊓ ⊔

Instead of taking J 1/2 := 1 2 (M + M ′ ) in the preceeding proof, one could also have taken

J t := (1 -t)M + tM ′ , which is also in co ∂ B H(x) = ∂ C H(x) for any t ∈ [0, 1]
. The inconvenient of this latter choice, when t = 1/2, is that M is usually not known. In particular, it is not necessarily known whether M i, : may be A i, : or B i, : , for a particular i ∈ E = (x), while J t depends on this value when t = 1/2. In contrast, J 1/2 has an explicit formula that does not require the knowledge of the value of M i, : for i ∈ E = (x).

Computation of the B-differential

This section describes techniques for computing a single Jacobian (section 5.1) or all the Jacobians (section 5.2) of the B-differential ∂ B H(x), in exact arithmetic, when H is the piecewise affine function given by (1.3). The algorithms are presented as tools for computing the sign vector set S ≡ S(V ), defined by (3.2) from a matrix V ∈ R n×p , which makes them appropriate, even when p > n.

When V is defined by (3.1), one has p n and the equivalence 3.5 tells us that S is then in bijection with ∂ B H(x), so that the algorithms actually compute Jacobians of the B-differential ∂ B H(x). The piece of software isf has been written to test the algorithms [START_REF] Dussault | ISF and BDIFFMIN -Matlab functions for central hyperplane arrangements and the computation of the B-differential of the componentwise minimum of two affine vector functions[END_REF][START_REF] Dussault | ISF and BDIFFMIN[END_REF].

Computation of a single Jacobian

An interest of the problem equivalence highlighted in proposition 3.4(3) is to provide a method to find rapidly an element of ∂ B H(x), which complements Qi's [60; 1993, final remarks (1)]. It is shown in [START_REF] Dussault | Partial description of the B-differential of the componentwise minimum of two vector functions by linearization[END_REF], that this method extends to the computation of an element of the B-differential in the nonlinear case, i.e., when H is the function H given by (1.6). The method is based on the following algorithm, which associates with p nonzero vectors v 1 , . . . , vp, which may be identical or colinear, a direction d such that v T i d = 0 for all i ∈ [1 : p]; it is a variant of the technique used in the proof of [78; lemma 2.1]. When the v i 's are also distinct, the direction d can also be derived from lemma 2.6, by adding the vector v 0 = 0.

Algorithm 5.1 (computes d ∈ R n such that v T i d = 0 for all i)
Let be given p nonzero vectors v 1 , . . . , vp in R n and take d ∈ R n \ {0}.

Repeat:

1.

If I := {i ∈ [1 : p] : v T i d = 0} = ∅, exit. 2. Let i ∈ I.
3. Take t > 0 sufficiently small such that, for all j / ∈ I, (v

T j d)(v T j [d + tv i ]) > 0. 4. Update d := d + tv i .
Explanation. In step 3, any sufficiently small t > 0 is appropriate (the proof of [78; lemma 2.1] computes bounds explicitely), since (v

T j d)(v T j [d + tv i ]) is positive for t = 0. The new direction d set in step 4 is such that v T i (d + tv i ) = t v i
2 > 0, so that this direction makes at least one more v T j d nonzero than the previous one. This implies that the algorithm finds an appropriate direction in at most p loops.

⊓ ⊔

The next procedure uses a direction d computed by algorithm 5.1 to obtain a single element of ∂ B H(x). Recall that the map σ is defined by (3.4a) and is a bijection from ∂ B H(x) onto S, defined by (3.2) (proposition 3.4).

Algorithm 5.2 (computes a single Jacobian in ∂ B H(x))

Let H be given by (1.3), x ∈ R n and suppose that p = 0.

1. Compute V ∈ R n×p by (3.1) and denote its columns by v 1 , . . . , vp ∈ R n .

By algorithm

5.1, compute d ∈ R n such that v T i d = 0 for all i ∈ [1 : p]. 3. Define s ∈ S by s i := sgn(v T i d), for i ∈ [1 : p]. 4. Then, σ -1 (s) ∈ ∂ B H(x).
1. Execute the recursive procedure stree-rec(V, +1).

Algorithm 5.4 (stree-rec (V, s)) Let be given V ∈ R n×p , with n and p ∈ N * , having nonzero columns, and a sign vector s ∈ S 1

k for some k ∈ [1 : p]. 1. If k = p, print s and return. 2. If (s, +1) ∈ S 1
k+1 , execute stree-rec(V, (s, +1)). 3. If (s, -1) ∈ S 1 k+1 , execute stree-rec(V, (s, -1)).

The method used to determine whether (s, ±1) is in S 1 k+1 depends on the specific algorithm and may or may not use a direction d intervening in (4.3). Note that, as emphasized in proposition 4.6(3), at least one of the sign vectors (s, +1) and (s, -1) belongs to S 1 k+1 (maybe both). It is justified not to explore the S-tree below an (s, ±1) that is not in S 1 k+1 , since then (s, ±1, s ′ ) / ∈ S for any s ′ ∈ {±1} p-k-1 . By construction, the algorithm stree prints all the elements of S 1 p ≡ S 1 := {s ∈ S : s 1 = +1} in step 1 of the stree-rec procedure. B. Rada and Černý's algorithm. The algorithm proposed by Rada and Černý [62; 2018], which is referenced below as the rc algorithm, deals with the determination of the cells associated with a general hyperplane arrangement. We describe it below for an arrangement of hyperplanes containing all zero (see section 3.4), which is the case when V result from (3.1) in the computation of the Bdifferential ∂ B H(x). We also use the linear algebra language of section 3.2.1, viewing the problem as the one of determining the set S defined by (3.1)); in contrast, the language used in [62] is more geometric. The algorithm builds the S-tree of the previous section A and, for each s ∈ S 1 k , it solves a single problem (LOP) to determine whether (s, +1) or (s, -1) is in S 1 k+1 .

+-- ←-S 1 2 ←-S 1 1 + +- ++ ←-S 1 3 = S 1 ++- +-+
The rc algorithm succeeds in solving only one LOP to determine whether (s, +1) and (s, -1)

are in S 1 k+1 , at the node s ∈ S 1 k , thanks to the memorization of a direction d such that s q (V T k d) > 0 (we note V k := V : ,[1 : k] ). Indeed, one has v T k+1 d < 0 =⇒ (s, -1) ∈ S 1 k+1 , v T k+1 d > 0 =⇒ (s, +1) ∈ S 1 k+1 ,
and one of these two cases takes place if we exclude the case where v T k+1 d = 0. In [62; Algorithm 1], the case where v T k+1 d = 0 is not dealt with completely since (s, +1) is declared to belong to S 1 k+1 in that case, while it is clear that (s, -1) is also in S 1 k+1 . Indeed, in our implementation of the rc algorithm, we modify slightly d by adding a small positive or negative multiple of v k+1 to d when v T k+1 d ≃ 0, so that both (s, ±1) are accepted in S 1 k+1 in that case. This choice may be at the origin of the differences that one observes in table 5.1 below between the statistics of the original rc algorithm in [62] and those of our implementation.

Next, when (s, s k+1 ) ∈ {±1} k+1 is observed to belong to S 1 k+1 , the question of whether (s, -s k+1 ) also belongs to S 1 k+1 arises. In the rc algorithm, the answer to this question is obtained by solving a LOP similar to

       min (d,t)∈R n ×R t s i v T i d 1, ∀ i ∈ [1 : k] -s k+1 v T k+1 d -t t -1. (5.2) When s ∈ S 1 k , this problem is feasible (take d satisfying s i v T i d 1, for all i ∈ [1 : k],
and t sufficiently large) and bounded (its optimal value is -1), so that it has a solution [START_REF] Chvátal | Linear Programming[END_REF][START_REF] Bonnans | Optimisation Numérique -Aspects théoriques et pratiques[END_REF][START_REF] Bonnans | Numerical Optimization -Theoretical and Practical Aspects[END_REF][START_REF] Charles | Selected Topics on Continuous Optimization -Version 2[END_REF]. Solving these LOPs is a time consuming part of the algorithms and in the numerical experiments of section 5.2.6, in particular in table 5.2, following [62], we measure the efficiency of the algorithms by the number of LOPs they solve.

One can now formally describe our version of the rc algorithm (the change is in step 2 of the rc-rec algorithm, which is not considered in the original rc algorithm). Algorithm 5.5 (rc (V )) Let be given V ∈ R n×p , with n and p ∈ N * , having nonzero columns.

1. Execute the recursive procedure rc-rec(V, v 1 , +1).

Algorithm 5.6 (rc-rec (V, d, s)) Let be given V ∈ R n×p , with n and p ∈ N * , having nonzero columns, a direction d ∈ R n and a sign vector

s ∈ {±1} k for some k ∈ [1 : p], such that s i v T i d > 0 for all i ∈ [1 : k]. 1. If k = p, print s and return. 2. If v T k+1 d ≃ 0, then 2.1. Execute rc-rec(V, d + , (s, +1)), where d + := d + t + v k+1 with t + > 0 chosen in the nonempty open interval   0, min i∈[1 : k] siv T i v k+1 <0 -v T i d v T i v k+1    . 2.2. Execute rc-rec(V, d -, (s, -1)), where d -:= d + t -v k+1 with t -< 0 chosen in the nonempty open interval    max i∈[1 : k] siv T i v k+1 >0 -v T i d v T i v k+1 , 0    .

3.

Else s k+1 := sgn(v T k+1 d).

3.1.

Execute rc-rec(V, d, (s, s k+1 )). 3.2. Solve the LOP (5.2) and denote by (d, t) a solution.

If t = -1, execute rc-rec(V, d, (s, -s k+1 )).

In steps 2.1 and 2.2, the minimum and maximum are supposed to be infinite if their feasible set is empty. One can check that the directions d ± computed in steps 2.1 and 2.2 are such that s i v T i d ± > 0 for i ∈ [1 : k + 1] and s k+1 = ±1, provided |v T k+1 d| is sufficiently small, which justifies the recursive call to rc-rec with the given arguments. The test v T k+1 d ≃ 0 done at the beginning of step 2 is supposed to take into account floating point arithmetic; admittedly it is not very rigorous, but the algorithm is designed to be as close as possible to the original rc algorithm in [62]; a more careful treatment of this situation is presented in section 5.2.4(B). The most time-consuming part of the rc algorithm comes from the possible numerous LOPs to solve in step 3.2 of rc-rec.

An algorithm using stem vectors

When s ∈ S k , it is conceptually easy to check whether (s, ±1) is in S k+1 , provided a list of all the stem vectors associated with V is known. Indeed, by proposition 3.10, if no subvector of (s, +1) (resp. (s, -1)) is a stem vector, then (s, +1) (resp. (s, -1)) belongs to S k+1 . Note also that, because any s ∈ S k has at least one descendant in the S-tree (proposition 4.6(3)), if it is observed that (s, +1) / ∈ S k+1 , then, necessarily, (s, -1) ∈ S k+1 . This observation prevents the algorithm from checking whether (s, -1) contains a stem vector, which is a time consuming operation when the list of stem vectors is large. For future reference, we formalize this algorithm below. Algorithm 5.7 (stem (V )) Let be given V ∈ R n×p , with n and p ∈ N * , having nonzero columns.

The second bottleneck is linked to the detection of a stem vector is the current sign vectors (s, ±1). This operation requires to examine the long list of stem vectors, which is a time consuming operation.

We shall see in the numerical experiments of section 5.2.6 that algorithm 5.7 is generally the fastest, provided the number of stem vectors is not too large.

Linear optimization problem and stem vector

The property described in this section will be useful for the improvement D 2 of the isf algorithm, described in section 5.2.4(D). It shows that a stem vector can be obtained easily from the dual solution of the linear optimization (LOP) (5.2), when (s, -s k+1 ) / ∈ S k+1 . Consider indeed the LOP (5.2) and denote by (d, t) one of its solutions (these have been shown to exist). Then, either t 0 (equivalently, (s, -s k+1 ) / ∈ S k+1 ) or t = -1 (equivalently, (s, -s k+1 ) ∈ S k+1 ). Let σ i , i ∈ [1 : k + 1], be the multipliers associated with the first k + 1 constraints of (5.2) and τ be the multiplier associated with its last constraint. Then, the Lagrangian dual of (5.2) reads [START_REF] Bonnans | Optimisation Numérique -Aspects théoriques et pratiques[END_REF][START_REF] Bertsekas | Nonlinear Programming[END_REF][START_REF] Bonnans | Numerical Optimization -Theoretical and Practical Aspects[END_REF][START_REF] Charles | Fragments d'Optimisation Différentiable -Théorie et Algorithmes. Lecture Notes[END_REF]

           max (σ,τ )∈R k+1 ×R i∈[1 : k] σ i -τ σ 0 τ 0 σ k+1 + τ = 1 σ k+1 s k+1 v k+1 = i∈[1 : k] σ i s i v i . ≡        max σ∈R k+1 i∈[1 : k+1] σ i -1 σ 0 σ k+1 1 σ k+1 s k+1 v k+1 = i∈[1 : k] σ i s i v i , (5.3)
where the second form of the dual is obtained by eliminating τ from the first form. By strong duality in linear optimization, the dual problems in (5.3) are feasible, have a solution and have the same optimal value as the primal problem. Let (σ, τ ) ∈ R k+1 × R be a dual solution. Then, (s, -s k+1 ) ∈ S k+1 if and only if t = -1 if and only if i∈[1 : k] σ i = 0 and σ k+1 = 0. We have shown that (s, -s k+1 ) ∈ S k+1 ⇐⇒ σ = 0. Therefore, (s, -s k+1 ) / ∈ S k+1 if and only if σ = 0 if and only if σ k+1 = 1 (if σ k+1 = 0, one can make the dual objective value as large as desired by multiplying σ by a factor going to +∞; if σ k+1 ∈ (0, 1), the dual objective would by increased by replacing σ by σ/σ k+1 ; in both cases the optimality of σ would be contradicted) if and only if τ = 0. We have shown that

(s, -s k+1 ) / ∈ S k+1 ⇐⇒ s k+1 v k+1 ∈ cone{s i v i : i ∈ [1 : k]}.
The next proposition shows how a matroid circuit can be detected from the dual solution σ when (s, -s k+1 ) / ∈ S k+1 .

Proposition 5.9 (matroid circuit detection) Suppose that (s, -s k+1 ) / ∈ S k+1 and that (σ, τ )

is a solution to the dual problem in the left-hand side of (5.3) located at an extreme point of its feasible set. Then, {i ∈ [1 :

k + 1] : σ i > 0} is a matroid circuit of V .
Proof We have seen that σ k+1 = 1 and τ = 0 when (s, -s k+1 ) / ∈ S k+1 . The fact that (σ, 0) is an extreme point of the feasible set of the problem in the left-hand side of (5.3) implies that the vectors [START_REF] Chvátal | Linear Programming[END_REF][START_REF] Charles | Fragments d'Optimisation Différentiable -Théorie et Algorithmes. Lecture Notes[END_REF] 0

s i v i i∈[1 : k], σi>0 ,

1

-s k+1 v k+1 are linearly independent.

In particular, the vectors

{s i v i : i ∈ [1 : k], σ i > 0} are linearly independent.
Proposition 5.10 (two descendants without optimization) Suppose that s ∈ {±1} k verifies

s q (V T k d) > 0, that v k+1 = 0 and that max i∈[1 : k] siv T i v k+1 >0 -v T i d v T i v k+1 < -v T k+1 d v k+1 2 < min i∈[1 : k] siv T i v k+1 <0 -v T i d v T i v k+1
.

(5.6)

1) The direction

d + := d + t + v k+1 verifies s q (V T k d + ) > 0 and v T k+1 d + > 0 if and only if t + is in the nonempty open interval    -v T k+1 d v k+1 2 , min i∈[1 : k] siv T i v k+1 <0 -v T i d v T i v k+1    .
(5.7a)

2) The direction

d -:= d + t -v k+1 verifies s q (V T k d -) > 0 and -v T k+1 d -> 0 if and only if t -is in the nonempty open interval    max i∈[1 : k] siv T i v k+1 >0 -v T i d v T i v k+1 , -v T k+1 d v k+1 2    .
(5.7b) C. Changing the order of the vectors v i 's. Each node s of the S-tree described in section 5.2.1(A) has one or two descendants: (s, +1) and/or (s, -1). Since there is at most one LOP solved per node of the S-tree, decreasing the number of nodes should decrease the number of LOPs to solve, which significantly count in the computing time. To reach that goal, one can try to get as much as possible at the top of the tree the nodes having a single descendant. As shown below, this can be achieved by changing the order in which the vectors v i 's, the columns of V , are considered in the depth-first search of the tree; previously, the order was imposed by the modification A, taking into account the rank of V . As we shall see, a new order is not fixed once and for all, but is determined during the construction of the S-tree, is reconsidered at each node and depends on the path going from the root of the S-tree to its leaves.

To implement this strategy, one associates with each node s ∈ S 1 k of the S-tree, k ∈ [1 : p -1], the list of vectors considered so far at that node, denoted by Ts := {i 1 , . . . , i k } ⊆ [1 : p]. Hence, we have to choose the next vector v i k+1 be selecting an index i k+1 in T c s := [1 : p] \ Ts. Now, a natural idea is to restrict the set of possible indices to T b s , the set of indices j of T c s for which one of the intervals (5.7a) or (5.7b), with v k+1 ≡ v j , is empty (implying that the technique used in the modification B will not give two descendants), if there is such an index, or T c s otherwise. To determine the index in T b s , we take

i k+1 = arg max i∈T b s |v T i d| v i , (5.8)
which favors the vectors v i for which |v T i d|/ v i is away from zero. As table 5.2 indicates (section 5.2.6(C.3)), this modification has a significant impact on the decrease of the number of LOPs to solve. D. Using stem vectors. We present in this section various modifications that use the concept of stem vector, introduced in the second part of section 3.2.2. These stem vectors are used to detect infeasible sign vectors, i.e., elements of S c , thanks to proposition 3.10. If s ∈ S 1 k and (s, s k+1 ) ∈ S c for s k+1 ∈ {±1}, s has no descendant in S along (s, s k+1 ), so that this part of the S-tree does not need to be explored. From this point of view, computing all the stem vectors looks attractive, but, to our knowledge, this is a time consuming process, so that this option is not necessarily the most efficient one. The modifications presented below use more and more stem vectors, whose computation requires more and more time. D 1 ) Natural candidates as stem vectors are those obtained from the matroid circuits I made of r + 1 columns of V (r = rank(V )) formed of the r linear independent columns selected by the QR factorization of section 5.2.4(A) and one of the remaining pr columns of V . By proposition 3.11, such I contains exactly one circuit. Therefore, one detects in this way pr circuits and 2(pr) stem vectors. This is not much compared to the total number of stem vectors, which may depend exponentially on p, so that the number of infeasible sign vectors detected by these stem vectors is usually relatively small (see table 5.2).

D 2 ) With this option, when a LOP (5.2) is solved at a certain node s ∈ S 1 k to see whether (s, s k+1 ) belongs to S 1 k+1 , for s k+1 ∈ {±1}, the dual solution is used to determine a matroid circuit, as shown by proposition 5.9. For this purpose, the isf code solves the LOP with the dual-simplex algorithm, so that the computed dual solution is at a vertex of the dual feasible set.

D 3 ) With this option, all the stem vectors are computed, before running the recursive process that builds the S-tree. At each note s ∈ S 1 k , the algorithm still computes a direction d ∈ R n such that s i v T i d > 0 for all i ∈ Ts (the set of vector indices considered so far at s). The advantage of this direction is to allow the algorithm to use the beneficial modifications B and C and to easily determine one or two signs s k+1 ∈ {±1} such that (s, s k+1 ) ∈ S 1 k+1 . If a single sign s k+1 ∈ {±1} is selected, the stem vectors can decide whether (s, -s k+1 ) ∈ S 1 k+1 . If this is the case, this option D 3 has the inconvenient of still requiring to solve a LOP to get a direction associated with (s, -s k+1 ). These LOPs (5.2) have an optimal value -1 and should not be solved exactly. Indeed, as soon as a feasible direction d for (5.2) gives a negative value to the objective of the problem, one could stop solving it, since this d verifies s i v T i d > 0 for all i ∈ T (s,-s k+1 ) . We have not implemented that inexact solve of the LOPs, by lack of flexibility of the solver Linprog in Matlab. D 4 ) Like with the option D 3 , all the stem vectors are computed, before running the recursive process that builds the S-tree. But now, unlike with option D 3 , the algorithm computes no direction d ∈ R n . When option A is also activated, the resulting approach can be viewed as an improvement of the algorithm 5.7 (stem) presented in section 5.2.2.

Note that, knowing all the stem vectors, one could compute the complementary set S c rather easily by completing with ±1 the unspecified components of the stem vectors. Next, S could be obtained from S c by taking its complementary set in {±1} p , but a straightforward implementation of this last operation looks rather expensive, so that we have not experimented it numerically.

Isf algorithm

We have named isf (for Incremental Signed Feasibility) the algorithm that improves the rc algorithm 5.5 or the stem algorithm 5.7 with the enhancements described in section 5.2.4. For the purpose of precision and reference, we formally state it in this section. It would be cumbersome and confusing, hence inappropriate, to mention all the options in its description, in particular because all of them have been specified separately in the previous section. As an example of algorithm, we provide a description with the options ABCD 2 . It starts with a hat procedure isf, similar to that of the rc algorithm but with the additional easy determination of Sr (modification A) and the computation of some stem vectors (modification D 1 ). Then, the hat procedure calls the recursive procedure isf-rec.

Algorithm 5.11 (isf (V ), with options ABCD 2 ) Let be given V ∈ R n×p , with n and p ∈ N * , having nonzero columns.

algorithm and the simulated rc algorithm implemented in the isf code respectively. The latter code will be used next, in the comparison with its improved versions, both regarding the LOP counters (table 5.2) and the CPU times (table 5.3). Table 5.1 Description of the test-problems and comparison of the "original rc algorithm in [62]", written in Python, and the "simulated rc algorithm 5.5", written in Matlab: "(n, p, r, ς)" are the dimensions of the problem (V ∈ R n×p is of rank r and has ς circuits), "|∂ B H(x)|" is the cardinality of the B-differential of H given by (1.3), "Schläfli's bound" is the right-hand side of (4.8), "Original rc" gives the number of linear optimization problems (LOPs) solved by the original piece of software in Python of Rada and Černý [62], "Simulated rc" gives the number of LOPs solved by the implementation in the Matlab code isf of the Rada and Černý algorithm (see algorithm 5.5), "Difference" is the difference between the two previous columns. Note (1): interruption of the run after several days of computation.

1)

The randomly generated problems rand are likely to provide vectors v i 's (the columns of V ) in general position, in the sense of definition 4.9. This can be seen indirectly on the numbers in table 5.1.

r It is known from proposition 4.10 that (4.7) implies equality in (4.8). This equality indeed holds, as we can observe by comparing columns 6 and 7.

r The number of circuits is also predictable for the random problems. Indeed, by the random generation of V , a subset of columns is likely to have nullity 1 (i.e., to form a circuit, in the matroid terminology) if and only if it contains r + 1 columns (r being the rank of V ). Therefore, their number should be ( p r+1 ) (see also [43; footnote 1]), which is indeed the number displayed in the 5th column of table 5.1 for the randomly generated problems.

2) The number of matroid circuits, given in the column labeled by ς, depends on the determination of the nonzero elements of the normalized vector α ∈ N (V : ,I ) \ {0} for the selected index set I (proposition 3.11). This operation is sensitive to a threshold value that is set to 10 5 ε, where ε > 0 is the machine epsilon; smaller values for this threshold have occasionally given larger numbers of matroid circuits. In other words, due to the floating point calculation, there is no certainty that the given number of circuits is the one that would be obtained in exact arithmetic. With a computation in rational numbers, this difficulty is avoided [START_REF] Dussault | Primal and dual approaches for the chamber enumeration of hyperplane arrangements[END_REF].

3) A comparison between the "Original rc code" in Python and its "Simulated rc code" in Matlab

shows that the latter is slightly more effective in terms of the number of LOPs solved. This is probably due to the special treatment in step 2 of the case where v T k+1 d ≃ 0 in algorithm 5.6, which is not considered in the original code.

C. Observations on table 5.2. Table 5.2 shows the effect of the modifications discussed in section 5.2.4 on the number of LOPs solved, which significantly counts in the computing time. This will lead us to select three algorithms, those which bring the best profit on the LOP counter. The columns labeled "Ratio" show the acceleration ratio with respect to the simulated rc code in terms of LOPs, that is the ratio of the LOP counter of the considered algorithm divided by the LOP counter of the simulated rc algorithm. On the last two lines of the table, one finds the mean and median values of these acceleration ratios, which may be viewed as a summary of the effect of the considered modification. These mean/median values must be taken with caution when a solver fails to solve a problem as is the case with isf(ABCD 3 ) and isf(AD 4 ) on problem bek-resonance-6.

1)

The modification A, proposed in section 5.2.4(A), which uses the QR factorization to get r linearly independent columns of V , does not bring a large benefit ("Ratio" is close to 1) and sometimes increases the number of LOPs to solve. The benefit is not important since it "only"

prevents i∈[0 : r-1] 2 i = 2 r -1 nodes from running the LOP solver, which is usually a small fraction of the total number of nodes of the S-tree. One also observes that the number of solved LOPs may increase (acceleration ratio < 1), which is sometimes due to the fact that the number 2 r-1 of nodes at level r with modification A is larger than the one without modification A, which contributes to increase the total number of nodes of the constructed S-tree and, therefore, tends to increase the number of LOPs to solve. Furthermore, the order in which the vectors are considered without/with modification A is not identical, which has also an impact on the number of solved LOPs (see section 5.2.4(C)).

2) The modification B, proposed in section 5.2.4(B), which is able to detect two descendants of an S-tree node, without solving any LOP, has a significant impact on the total number of these problems. We see, indeed, that the (mean, median) acceleration ratio is raised to (1.28, 1.17).

3) Consider now the modification C, described in section 5.2.4(C), which changes the order in which the vectors v i 's are considered. We use the test-problem rand-7-13-5 to show its effect in the next table. The table gives the number of nodes for each level in the S-tree, with the modifications AB and with the modifications ABC. Since rank(V ) = 5 for this problem and since the modification A is used in both cases, the number of nodes per level, only starts to differ from level 6 (before that it is equal to 2 l-1 , where l is the S-tree level). The final level is 13 (since there are p = 13 vectors) and its number of leaves is |S|/2 = 794 (an observation from the table above), necessary identical in both cases. The effect of the modification C can be seen on the smaller number of nodes per level and in all the S-tree (rightmost column). This contributes to the decrease of the number of LOPs to solve: the (mean, median) acceleration ratio is raised to (3.45, 2.51). 5.2 Evaluation of the efficiency of the solvers by the number of LOPs they solve: A (taking the rank of V into account), B (special handling of the case where v T k+1 d ≃ 0), C (changing the order of the vectors v i 's by taking i k+1 by (5.8)), D 1 (pre-computation of 2(p-r) stem vectors after the QR factorization), D 2 (D 1 and 2 additional stem vectors computed after solving a LOP, whose optimal value is nonnegative), D 3 (all the stem vectors are first computed and, for (s, ±1) ∈ S k+1 , a LOP is solved to get a handle d), D 4 (all the stem vectors are first computed and no LOP is solved). The "Ratio" (acceleration ratio) columns give for each considered problem the ratio (LOPs of the considered isf version)/(LOPs of simulated rc). Note (1): interruption of the run after several days of computation. The Mean/Median rows give the mean and medan values of the ratios.

4) The modifications D, described in section 5.2.4(D), deal with the contribution of the computed stem vectors, whose number increases from modification D 1 (2(p-r) stem vectors after the QR factorization of V ), D 2 (more stem vectors from the dual solution of the LOP (5.2) when this one has a nonnegative optimal value), D 3 and D 4 (all the stem vectors).

r We see that the option D 1 yields already some improvement (less LOPs to solve), but not much, raising the (mean, median) acceleration ratio from (3.45, 2.51) to (3.88, 2.95).

r The use of the option D 2 is more beneficial since the (mean, median) acceleration ratio now goes up to (31.54, 11.65). We understand this fact to have its origin in the increase in the number of stem vectors detected from the dual solutions of some solved LOP. Note that this last operation does not require much computation time.

r With option D 3 , only the LOPs (5.2) with the optimal value -1 are solved, while, with option D 4 , no LOP is solved. The efficiency of these modifications largely depends on the total number 2ς of stem vectors. If this one is not too large, the modifications have an important benefit. Otherwise, it can lead to execution failure, as for problem bek-resonance-6, which requires days of computation.

In conclusion of these observations, one could retain the following three solvers for a comparison on their computing time.

r Isf(ABCD 2 ) is the most efficient solver that does not compute all the stem vectors. r The solvers isf(ABCD 3 ) and isf(AD 4 ) cannot be compared with the other solvers on the results of table 5.2 since both use all the stem vectors, so that the time to compute and use these must be taken into account, and isf(AD 4 ) does not solve any LOP, which is the measure of efficiency in table 5.2. D. Observations on table 5.3. Measuring the efficiency of the algorithms by the number of LOPs solved during execution, like in table 5.2, is sometimes misleading. If this is the main cost item for some algorithms, it is no longer the case when a large amount of stem vectors is computed. For two reasons. First, the time spent in the computation of these stem vectors is not negligible, far from it, at least in our implementation, in which each of them requires the computation of the nullity of a matrix and a null space vector. Next, verifying that a sign vector contains a stem vector (proposition 3.10) is also time consuming when there are many stem vectors. Therefore a comparison of the CPU time of the runs is welcome. This is done for a selection of versions of the isf codes in table 5.3, those selected at the end of section 5.2.6(C). Here are some observations on the statistics of this table.

1) A first observation is that the good behavior of the selected versions of the isf codes is confirmed, even though the acceleration ratios are not as large as the one based on the number of LOPs solved. This can be explained by the fact that the time spent in solving LOPs is counterbalanced by the handling of stem vectors for the versions ABCD 3 and AD 4 . Anyway, one observes that the CPU time acceleration ratios have (mean, median) values in the ranges (15..31, 5..20), which is significant.

2) The most effective combination of code options depends actually on the considered problems.

It is difficult to state a rule that would predict which code behaves best because some solvers are better on some phases of the run, but worse on others (the three main phases are the detection of the stem vectors, the execution of LOPs and the search for stem vectors covered by a given sign vector). However, an inductive rule manifests itself: the purely dual method AD 4 is ahead for problems with a reasonable number of stem vectors (or matrix circuits), but can require a too large number of computing time if this number becomes large (this is the case of problems bek-threshold-6, bek-resonance-6 and bek-demicube-7). This conclusion could be invalidated if better techniques are used to enumerate and use the stem vectors. Table 5.3 Evaluation of the efficiency of the solvers by their computing times. The "Ratio" (acceleration ratio) columns give for each considered problem the ratio (Time of the considered isf version)/(Time of simulated rc). Note (1): interruption of the run after several days of computation. The Mean/Median rows give the mean and medan values of the ratios.

Discussion

This paper deals with the description and computation of the B-differential of the componentwise minimum of two affine vector functions. The fact that this problem has many equivalent formulations, some of them being highlighted in section 3, implies that the present contribution has an impact on several domains, including on the description of the arrangement of hyperplanes in the space. To this respect, a singular aspect of this contribution is to propose a dual approach to solve the problem, using some or all the stem vectors, a concept made useful thanks to the convex analysis tool that is Gordan's alternative. Besides this contribution, the paper also brings various improvements of an algorithm of Rada and Černý [62], which was designed to determine the cells of an arrangement of hyperplanes in the space. Even in the spirit of the methods proposed in this article, there is still room for improvement, in relation to three identified bottlenecks: (i) we have mentioned that with the option D 3 , the LOP (5.2) can be solved inexactly, since, in that case, the optimal value is -1, while any negative objective value for a feasible unknown would suffice, but this requires a better tuning of the linear optimization solver, (ii) computing more efficiently all the stem vectors (or matroid circuits) of the matrix V is certainly a source of improvement, (iii) a better algorithm to decide more rapidly that a sign vector contains a stem vector is also welcome. Some of these possible improvements are also linked to a better choice of programming language, probably one using a compilation phase. This contribution has also various possible extensions. A first one would be to develop a dual approach to the problem of the arrangement in the space of hyperplanes having no point in common [START_REF] Dussault | Primal and dual approaches for the chamber enumeration of hyperplane arrangements[END_REF]. Another natural extension would be to see the implications of this work for computing the B-differential of the componentwise minimum of nonlinear vector functions [START_REF] Dussault | Partial description of the B-differential of the componentwise minimum of two vector functions by linearization[END_REF]. Finally, the possibility to take profit of the computation of the full B-differential of the function H in (1.3) in a Newton-like approach to solve (1.4) is a subject that deserves reflection.

( 1 . 4 )

 14 is a reformulation of the balanced [28] Linear Complementarity Problem (LCP) 0 (Ax + a) ⊥ (Bx + b) 0.

Problem 3 . 1 (

 31 B-differential of the minimum of two affine functions) Let be given two positive integers n and m ∈ N * , two matrices A, B ∈ R m×n and two vectors a, b ∈ R m . It is requested to compute the B-differential at some x ∈ R n of the function H : R n → R m defined by (1.3).

Problem 3 . 6 (

 36 orthants encountered by the null space of a matrix) Let be given two positive integers n and p ∈ N * and a matrix V in R n×p with nonzero columns. Associate with I ⊆ [1 : p] the following orthant of R p : O p I := {y ∈ R p : y I 0, y I c 0}, where I c := [1 : p] \ I. It is requested to determine the set I := {I ⊆ [1 : p] : N (V ) ∩ O p I = {0}}. ⊓ ⊔ Note that, if I ∈ I, then I c ∈ I (because y ∈ (N (V ) ∩ O p I ) \ {0} implies that -y ∈ (N (V ) ∩ O p I c ) \ {0}), so that |I| is even (just like |S| and |S c |, see proposition 4.1).

(3. 7 )

 7 Proposition 3.7 (bijection S c ↔ I) The map ı defined by (3.6) is a bijection from S c onto I. Proof Let s ∈ {±1} p and set I := ı(s) = {i ∈ [1 : p] : s i = +1}. Define A := Diag(s)V T to make the link with Gordan's alternative (3.7). One has the equivalences

Since ı : 3 ∅ and O 3 [ 1 : 3 ]

 ı3313 {±1} p → P([1 : p]) is a bijection, the above equivalences show that ı is also a bijection from S c onto I. ⊓ ⊔ Equivalence 3.8 (S c ↔ I) The equivalence between problems 3.3 and 3.6 is a consequence of the bijectivity of ı : S c → I, established in proposition 3.7: to determine S, it suffices to determine S c = ı -1 (I), hence to determine I, and vice versa. ⊓ ⊔ In example 3.2, one has N (V ) = Re, which only encounters the orthants O outside the origin; hence I = {∅, [1 : 3]}. We have seen that S c = {±(1, 1, 1)} for this problem. Clearly, ı maps S c onto I bijectively, as claimed in proposition 3.7.

(3. 11 )

 11 Proof [⇒] The index set J ⊆ [1 : p] in the right-hand side of (3.11) can be determined as one satisfying the following two properties:

8 r = 3 , 3 , |S| = 14 Fig. 3 . 2

 8331432 Fig.3.2 Linearly separable bipartitions of a set of p = 4 points vi in R 2 (the dots in the figure). Possible separating hyperplanes are the drawn lines. We have not represented any separating line associated with the partition (∅, [1 : p]) or ([1 : p], ∅), so that |S| = 2(ns + 1), where ns is the number of represented separating lines. We have set r := dim(vect{v 1 , . . . , vp}) + 1.

Figure 3 .

 3 Figure 3.3 illustrates problem 3.19 for the linear complementarity problem 3.2. It is requested to list the regions of R n that are separated by these hyperplanes, which are the connected components of R n \ ( i∈[1 : p] H i ). Such a region is called a cell or a chamber, depending on the authors [5, 69, 2]. More specifically, let us define the half-spaces

  These equivalences show the bijectivity of η from C onto S. ⊓ ⊔ Equivalence 3.21 (signed linear system feasibility ↔ hyperplane arrangement) The equivalence between problems 3.3 and 3.19 follows from the bijection of the map η : C → S claimed in proposition 3.20. ⊓ ⊔ 4 Description of the B-differential This section gives some elements of description of the B-differential ∂ B H(x), when H is the piecewise affine function given by (1.3) and x ∈ R n . This description is often carried out in terms of the matrix V defined by (3.1), whose p columns are denoted by v 1 , . . . , vp ∈ R n and are nonzero by construction. When the properties are given for S, one may have p n and the referenced matrix V ∈ R n×p is assumed to have nonzero columns, which implies that S = ∅. Some properties of ∂ B H(x) are given in section 4.1, including those that are useful in [31]. Section 4.2 deals with the cardinality |∂ B H(x)| of the B-differential. Section 4.3 analyzes more precisely two particular configurations. Section 4.4 highlights two links between the B-differential and the C-differential of H.

Proposition 4 . 1 (

 41 symmetry of ∂ B H(x)) Suppose that p > 0. Then, the B-differential ∂ B H(x)is symmetric and |∂ B H(x)| is even.

Definition 4 . 3 (

 43 adjacency in {±1} p ) Two sign vectors s 1 and s 2 ∈ {±1} p are said to be adjacent if they differ by a single component (i.e., the vertices s 1 and s 2 of the cube co{±1} p can be joined by a single edge).⊓ ⊔ Definitions 4.4 (connectivity in {±1} p ) A path of length l in a subset S of {±1} p is a finite set of sign vectors s 0 , . . . , s l ∈ S such that s i and s i+1 are adjacent for all i ∈ [0 : l -1]; in which case the path is said to be joining s 0 to s l . One says that a subset S of {±1} p is connected if any pair of points of S can be joined by a path in S. ⊓ ⊔ Proposition 4.5 (connectivity of the B-differential) The set S defined by (3.2) is connected if and only if V has no colinear columns. In this case, any points s and s of S can be joined by a path of length l := i∈[1 : p] |s is i |/2 p in S.

Proposition 4 . 6 (incrementation) 1 )

 461 If s ∈ S c k , then (s, ±1) ∈ S c k+1 . In particular, |S c k+1 | 2|S c k |. 2) If v k+1 / ∈ vect{v 1 , . . . , v k }, then, (s, ±1) ∈ S k+1 for all s ∈ S k . In particular, |S k+1 | = 2|S k | and |S c k+1 | = 2|S c k |.

[ 1 :. 3 .

 13 p] could be excessively expensive. Therefore, having simple-to-compute lower and upper bounds on |∂ B H(x)| may be useful in some circumstances, including theoretical ones. Proposition 4.7 gives elementary lower and upper bounds, while proposition 4.10 reinforces the upper bound, thanks to a lower semicontinuity argument (proposition 4.8). Necessary and sufficient conditions ensuring equality in the left-hand side or right-hand side inequalities in the next proposition are given in section 4Proposition 4.7 (lower and upper bounds on |∂ B H(x)|) For V given by (3.1) and r := rank(V ), one has max(2p, 2 r ) 2 r + 2(pr) |∂ B H(x)| 2 p . Proof The first inequality is clear since p r 1 and 2r 2 r . Consider the second inequality. One can assume that the first r columns of V are linearly independent, so that |Sr| = 2 r (notation (4.3) and proposition 4.6(2)). Next, by proposition 4.6(3), |S r+1 | 2 r + 2. By induction, the given lower bound holds for |Sp| = |S| = |∂ B H(x)|.

Proposition 4 . 8 (

 48 lower semicontinuity of |∂ B H(x)|) Suppose that the set S, defined by (3.2),

  a necessary and sufficient condition to have completeness of ∂ B H(x).Here follows a less restrictive assumption, called general position, which is equivalent to have equality in (4.8) below. In connection with this assumption, it is worth noting that, for a matrix V ∈ R n×p of rank r, one has ∀ I ⊆ [1 : p] : rank(V : ,I ) min(|I|, r).

( 4 . 6 )

 46 Definition 4.9 (general position) The vectors v 1 , . . . , vp ∈ R n are said to be in general position, if the matrix V := v 1 • • • vp verifies ∀ I ⊆ [1 : p] : rank(V : ,I ) = min(|I|, r), (4.7)

figure 3 .

 3 figure 3.2 since |∂ B H(x)| = 12, while the right-hand side of (4.8) reads 2[( 3 0 ) + ( 3 1 ) + ( 3 2 )] = 14 (we

( 4 .

 4 9a) It certainly suffices to show that one can find a matrice Ṽ ⊆ R n×p of rank r arbitrarily close to V that satisfies |S( Ṽ )| = β, (4.9b) since then one would have the expected contradiction with the lower semicontinuity of V → |S(V )| ensured by proposition 4.8:

  r), one finally obtains a matrix Ṽ , arbitrarily close to V , such that (4.7) holds: rank(( Ṽ ) : ,I ) = min(|I|, r) for all I ⊆ [1 : p].

by proposition 4 . 8 .

 48 If the columns of V are in general position, proposition 4.10 tells us that |S(V )| = β, where β is the right-hand side of Schläfli's bound (4.8) with r = rank(V ). Now, by the fact that rank( Ṽ ) = r, proposition 4.10 ensures that |S( Ṽ )| β. Therefore, one must have S( Ṽ ) = S(V ).

Proposition 4 .

 4 [START_REF] Bonnans | Numerical Optimization -Theoretical and Practical Aspects[END_REF] (the central C-Jacobian) One has J ∈ ∂ C H(x) for the Jacobian whose ith row, i ∈ [1 : m], is defined by

Fig. 5 . 1

 51 Fig. 5.1 Half of the S-tree for example 3.2 (the other half is obtained by swapping the +'s and the -'s). Top-down arrows indicate descendence; the sign sets S 1 k are defined by (5.1).

  •, • ). It is a refinement of [78; lemma 2.1]. Suppose that (E, •, • ) is a Euclidean vector space, p ∈ N * and v 1 , . . . , vp are p distinct vectors of E. Then, the set of vectors ξ ∈ E such that |{ ξ, v i

	Lemma 2.6 (discriminating covectors)

  Problem 3.19 (arrangement of hyperplanes containing the origin) Let be given two positive integers n and p ∈ N * and p nonzero vectors v 1 , . . . , vp ∈ R n . Consider the hyperplanes containing the origin:

	by (1.1) and (1.3).	3.2) or ∂ B H(x) defined

  |I|,with the notation (3.2),(3.19) and(3.6), is a tricky task, even in the present affine case, since it subtly depends on the arrangement of the vectors v i 's in the space (see figure3.2). Many contributions have been done on this subject; the earliest we cite dates from 1826[START_REF] Steiner | Einige Gesetze über die Theilung der Ebene und des Raumes[END_REF][START_REF] Roberts | On the figures formed by the intercepts of a system of straight lines in a plane, and on analogous relations in space of three dimensions[END_REF][START_REF] Grünbaum | Arrangements and Spreads[END_REF][START_REF] Zaslavsky | Facing up to arrangements: face-count formulas for partitions of space by hyperplanes[END_REF][START_REF] Las | Matroïdes orientables[END_REF] 3, 4,[START_REF] Edelsbrunner | Algorithms in Combinatorial Geometry[END_REF][START_REF] Kenneth | Applications of random sampling in computational geometry[END_REF][START_REF] Stanley | An introduction to hyperplane arrangements[END_REF] 2]. The formula (4.5) for |∂ B H(x)| is due to Winder

  If p 2 and the vectors v i 's are not two by two colinear, one has rank(V ) = 2 if and only if |∂ B H(x)| = 2p.

	⊓ ⊔
	Proposition 4.13 (fan arrangement) Proof [⇒] A short proof leverages Schläfli's bound (4.8) with equality. Since the v i 's are not two by two colinear, one has for any I ⊆ [1 : p]:

takes the value 2 or p; they yield the lower and upper bounds on |∂ B H(x)| given by proposition 4.7. The lower bound 2p applies to the left-hand side pane of figure 3.2. As shown by the intermediate pane in figure 3.2, however, |∂ B H(x)| does not only depend on r. Proposition 4.12 (injective matrix) The matrix V ∈ R n×p given by (3.1) is injective if and only if |∂ B H(x)| = 2 p . Proof Indeed, by proposition 4.2, the B-differential ∂ B H(x) is complete (meaning that it is equal to ∂ × B H(x), given by (2.3)) if and only if V is injective. Clearly, the completeness of ∂ B H(x) is equivalent to |∂ B H(x)| = 2 p .

Table

			Number of linear optimization problems (LOPs) solved and acceleration ratio (Ratio) for various options
		Simulated	isf (A)		isf (AB)	isf (ABC)	isf (ABCD1)	isf (ABCD2)	isf (ABCD3)	isf (AD4)
	Problem	rc	LOP Ratio	LOP Ratio	LOP Ratio	LOP Ratio	LOP Ratio	LOP Ratio LOP Ratio
	rand-8-15-7	9907	9844	1.01	7641	1.30	5210	1.90	5199	1.91 3.19	0	-
	bek-resonance-4	635	672	0.94	546	1.16	171	3.71	138	4.60	31	20.48	0	-	0	-
	bek-resonance-5	36311	37607	0.97	34056	1.07	6700	5.42	6569	5.53	1141	31.82	810	44.83	0	-
	bek-resonance-6	6164040 6269410	0.98 5956586	1.03 760930	8.10 760457	8.11 155555	39.63	(1)	-	0	-
	bek-crosspolytope-11	86526	110418	0.78	58954	1.47	17569	4.92	15265	5.67	6085	14.22	6049	14.30	0	-
	bek-crosspolytope-12	260601	337910	0.77	182575	1.43	46900	5.56	41780	6.24	18785	13.87 18740	13.91	0	-
	bek-crosspolytope-13	788970 1028066	0.77	560013	1.41 124828	6.32 113564	6.95	57299	13.77 57244	13.78	0	-
	bek-demicube-5	99	90	1.10	33	3.00	24	4.12	12	8.25	3	33.00	0	-	0	-
	bek-demicube-6	4719	4761	0.99	3659	1.29	1882	2.51	1741	2.71	665	7.10	588	8.03	0	-
	bek-demicube-7	674663	704553	0.96	623160	1.08 175870	3.84 175595	3.84	60876	11.08 58333	11.57	0	-
	Mean			0.97		1.28		3.45		3.88		31.54		24.52		-
	Median			1.00		1.17		2.51		2.95		11.65		5.54		-

Acknowledgments

We thank Černý and Rada for providing their code and test problems, those used in [62]; part of these were used in the numerical experiments. We also thank the referees for their remarks and recommendations, which have helped us make the paper more readable.

Explanation. When p = 0, ∂ B H(x) = ∂ × B H(x) contains a single Jacobian that is given by (2.3), which explains why algorithm 5.2 focuses on the case when p > 0. The sign vector s computed in step 3 is such that s i v T i d > 0 for all i ∈ [1 : p], so that it is indeed in S and, by proposition 3.4, σ -1 (s) is a Jacobian in ∂ B H(x).

⊓ ⊔

Computation of all the Jacobians

This section presents two basic algorithms, and some more efficient variants, for computing all the B-differential of H. They use the notion of S-tree presented in section 5.2.1(A). The first algorithm is grounded on the notion of stem vector (section 3.2.2) and is described in section 5.2.2. The second algorithm is the outcome of a series of improvements brought to an algorithm by Rada and Černý [62; 2018] (section 5.2.1(B)) for computing the cells of a hyperplane arrangement, which is known to be an equivalent problem to the one of computing the B-differential of H when the hyperplanes contain zero (see section 3.4). The improvements are detailed in section 5.2.4 and the resulting algorithm is described in section 5.2.5. Finally, numerical experiments are presented in section 5.2.6 to compare the efficiency of the algorithms.

Algorithms for listing the elements of the finite set ∂ B H(x) can be designed by looking at one of the various forms of the problem, those described in section 3 and others [START_REF] Avis | Reverse search for enumeration[END_REF]; this is what we shall do. Most algorithms we have found in the scientific literature take the point of view of hyperplane arrangements of section 3.4 and can be used for more general arrangements than those needed to describe ∂ B H(x) (i.e., in which case the hyperplanes pass through zero). One can quote the contributions by Bieri and Nef [13; 1982], Edelsbrunner, O'Rourke and Seidel [35; 1986], Avis and Fukuda [5; 1996], improved by Sleumer [69; 1998], and, more recently, Rada and Černý [62; 2018], which is described in section 5.2.1(B). See also [START_REF] Dussault | Primal and dual approaches for the chamber enumeration of hyperplane arrangements[END_REF].

Incremental-recursive algorithms

The algorithms described in this section are incremental in the sense that the considered sign vectors have their length increased by one at each step. Furthermore, the algorithms explore the S-tree described in subsection A below by recursive procedures, whose names are recognizable by their suffix "-rec". All the procedures end by returning to their calling program.

A. The S-tree. A common feature of the algorithms considered in this paper is the construction of the S-tree described below, incrementally and recursively. This idea was probably introduced by Rada and Černý [62; 2018].

The level k of the S-tree is formed of a set of sign vectors denoted by

where S k is the subset of {±1} k defined by (4.3). In particular, the level 1 or root of the S-tree contains the unique sign vector +1 ∈ {±1} 1 . The S-tree has p levels, where p is the number of vectors v i , or columns of the given matrix V ∈ R n×p . Note that there is no reason to compute {s ∈ S : s 1 = -1} since this part of S is equal to -{s ∈ S : s 1 = 1} by the symmetry property of S (proposition 4.1). In order to avoid the memorization of the elements of S 1 k , the S-tree is constructed by a depth-first search, which can be schematized as follows.

Algorithm 5.3 (stree (V )) Let be given V ∈ R n×p , with n and p ∈ N * , having nonzero columns.

This algorithm is improved below, as the option AD 4 of the isf algorithm (see paragraphs A and D of section 5.2.4).

Note that, this algorithm need not generate directions d satisfying s q (V T k d) > 0, like the rc algorithm and need not solve any linear optimization problem. Nevertheless, regarding the computation time, the algorithm has two bottlenecks that we now describe.

The first bottleneck comes from the fact that the algorithm must compute all the stem vectors (or the set C of matroid circuits in (3.9)) associated with V . This is usually an expensive operation [START_REF] Khachiyan | On the complexity of some enumeration problems for matroids[END_REF][START_REF] Mary | Efficient enumeration of solutions produced by closure operations[END_REF]64]. For example, if V is randomly generated and of rank r, like in the test-cases data_rand_* in the experiments of section 5.2.6, any selection of r columns of V is likely to form an independent set of vectors, so that C is likely to be the sets of column indices of size r + 1. In this case, the number of circuits is likely to be the combination ( p r+1 ) (and it is actually that number, see section 5.2.6(B.1)), which can be exponential in p (this number is bounded below by 2 p/2 /(p + 1) if p is even and r + 1 = p/2 [24; (11.52)]). In the implemented isf code, numerically tested in section 5.2.6, only the sets of columns whose cardinality is in [3 : r + 1] are examined (since any group of two columns of V is supposed to be linearly independent and a group of r + 2 columns or more is of nullity 2, hence such group cannot form a matroid circuit; see (3.9)).

Since

The conclusion of the proposition follows from proposition 3.11.

⊓ ⊔

Recall that the dual-simplex algorithm finds a dual solution at an extreme point of the dual feasible set. For this reason, we use this approach in the isf algorithm with option D 2 (see section 5.2.4(D)).

Improvements of the rc and stem algorithms

This section presents several modifications of the rc algorithm and one modification of the stem algorithm that significantly improve their performance. The modifications are indicated by the letters A, B, C and D, with reference to the sections where they are introduced. Additional numeric indices specify variants of the D option. The version AD 4 (modifications A and D 4 ) can be considered as an improvement of the new algorithm 5.7.

A. Taking the rank of V into account. Instead of starting with the vector s = +1, one can take into account the rank r := rank(V ) to determine 2 r initial vectors s, hence avoiding to solve linear optimization problems (LOPs) to determine these initial s's. This is especially useful when pr is small. In particular, when p = r, S is straightforwardly determined.

The algorithm selects r := rank(V ) linearly independent vectors v i , among the columns of V ∈ R n×p . These vectors can be obtained by a QR factorization of V P = QR, where P ∈ {0, 1} p×p is a permutation matrix, Q ∈ R n×n is orthogonal (i.e., Q T Q = In) and R ∈ R n×p is upper triangular with R [r+1 : n], : = 0. To simplify the presentation, one can assume, without loss of generality, that P = I, in which case the vectors v 1 , . . . , vr are linearly independent (in practice, the vectors are symbolically reordered by using the permutation matrix P ). By proposition 4.2 and with the notation (4.3):

(5.4) Furthermore, for each s ∈ Sr, we have, using S := Diag(s), Qr := Q : ,[1 : r] and Rr := R [1 : r],[1 : r] , that the vector ds = QrR -T r s

(5.5) is such that s q (V T : ,[1 : r] ds) = e > 0, as desired.

For each s ∈ Sr and the associated ds given by (5.5), the modified algorithm 5.5 runs the recursive function rc-rec(V, ds, s) (see algorithm 5.11 below).

B. Special handling of the case where v

k+1 d ≃ 0 like in step 2 of the rc-rec algorithm 5.6, but also when v T k+1 d is in the interval specified by (5.6) below. Note that the lefthand side in (5.6) is negative and the right-hand side is positive (this can be seen by multiplying numerators and denominators by s i and by using s i v T i d > 0 for all i ∈ [1 : k]), so that these inequalities are verified when v T k+1 d = 0. With the additional flexibility that (5.6) offers, the isf algorithm can sometimes avoid solving a significant number of LOPs of the form (5.2). For a proof of the next proposition, see [29]. identical normalized vectors (columns of V ), to identify nonzero components of circuits, etc. The Julia code described in [START_REF] Dussault | Primal and dual approaches for the chamber enumeration of hyperplane arrangements[END_REF], which deals with more general hyperplane arrangements, offers the user the possibility of requiring a computation in rational numbers, so as to have a computation in exact arithmetic.

We have assessed the codes on randomly generated problems (function rand in Matlab, names prefixed by rand and srand) and problems adapted/taken from [62] (names prefixed by rc) and [START_REF] Brysiewicz | Computing characteristic polynomials of hyperplane arrangements with symmetries[END_REF] (names prefixed by bek). Their relevant features are given in table 5.1 and their specifications are now given.

r The rand-n-p-r problems have their data formed of a randomly generated matrix V ∈ R n×p with prescribed rank r.

r For the problems srand-n-p-q, the first n columns of V ∈ R n×p form the identity matrix and the last pn > 0 columns have q nonzero random integer elements (0 < q pn), randomly positioned.

r The matrix V ∈ R n×p of problem rc-2d-n-p is formed of 4 blocs: V 1:2,1:n-2 = 0, V 3:n,n-1:p = 0, and the remaining blocks have random integer data.

r The problems rc-perm-n refer to the hyperplane arrangements that are called permutahedron in [62]: the matrix V ∈ R n×p is such that V : ,[1 : n] is the identity matrix and V : ,[n : p] is a Coxeter matrix [START_REF] Postnikov | Deformations of Coxeter hyperplane arrangements[END_REF] (each column is of the form e ie j for some i = j in [1 : n], where e k is the kth basis vector of R n ).

r The problems rc-ratio-n-p-r refer to the problems that are controlled by a degeneracy ratio ρ in [62]: the first n columns of the matrix V ∈ R n×p are randomly generated, while the other pn > 0 columns can either (with a probability ρ) be linear combination of the previously generated columns or randomly generated.

r The problems bek-threshold-n refer to the threshold arrangements in [17; § 6.2]: for n 2, each column of V ∈ R n×p is formed of the components of (1, w) where w ∈ R n-1 are all the vectors of {0, 1} n-1 (hence p = 2 n-1 ). This arrangement appears in the study of neural networks [76].

r The problems bek-resonance-n refer to the resonance arrangements in [17; § 6.3]: the columns of V ∈ R n×p are all the nonzero vectors with components in {0, 1} (hence p = 2 n -1). Note that, for this arrangement, the number of chambers (i.e., |S| in our notation) is only known for n 9. Our approach, which does not use the particular structure of this arrangement, can get |S| in a reasonable time on a laptop for n 6, which is to be compared to n 9 in [START_REF] Brysiewicz | Computing characteristic polynomials of hyperplane arrangements with symmetries[END_REF]. See [START_REF] Kühne | The universality of the resonance arrangement and its Betti numbers[END_REF] for applications.

r The problems bek-crosspolytope-n refer to the cross-polytope arrangements in [17; § 6.4]: for n 2, each column of V ∈ R n×p is formed of the components of (1, w) where w ∈ R n-1 are all the ±e i for i ∈ [1 : n -1]; hence p = 2(n -1). For these problems, one numerically observes that |S| = 2 1 3 n-1 -2 n-1 for n 12 (this observation is made for n 21 in [START_REF] Brysiewicz | Computing characteristic polynomials of hyperplane arrangements with symmetries[END_REF]).

r The problems bek-demicube-n refer to the demicube arrangements in [17; § 6.6]: the columns of V ∈ R n×p are the components of (1, w) where w ∈ {w ′ ∈ {0, 1} n-1 : i w ′ i is odd}. We have retained 3 problems per family, the most difficult that isf can solve in a reasonable time for the rc-perm and bek families. These test-problems are available on Software Heritage [START_REF] Dussault | ISF and BDIFFMIN[END_REF].

B. Observations on table 5.1. The dimensions n, p and r of the problems are given in columns 2-4 of table 5.1. Column 5 gives the number ς of matroid circuits of V . In columns 6 and 7, one finds the cardinality |∂ B H(x)| = |S| of the B-differential ∂ B H(x) and the Schläfli upper bound (the right-hand side of (4.8)). The codes will be compared on the number of linear optimization problems (LOPs) they solve, which is a good image of their computation effort, measured independently of the computer used to run the codes and the features of the LOP solver. A first example of comparison is given in columns 8 and 9 of table 5.1, where one finds the number of LOPs solved by the original rc