当前位置:首页 >> >> An output-sensitive algorithm for multi-parametric LCPs with sufficient matrices

An output-sensitive algorithm for multi-parametric LCPs with sufficient matrices


An output-sensitive algorithm for multi-parametric LCPs with su?cient matrices

arXiv:0807.2318v1 [math.OC] 15 Jul 2008

S. Columbano Institute for Operations Research ETH Zurich, Switzerland seba@ifor.math.ethz.ch

K. Fukuda Institute for Operations Research and Institute of Theoretical Computer Science ETH Zurich, Switzerland fukuda@ifor.math.ethz.ch

C.N. Jones Automatic Control Laboratory ETH Zurich, Switzerland cjones@ee.ethz.ch July 15, 2008

Abstract This paper considers the multi-parametric linear complementarity problem (pLCP) with suf?cient matrices. The main result is an algorithm to ?nd a polyhedral decomposition of the set of feasible parameters and to construct a piecewise a?ne function that maps each feasible parameter to a solution of the associated LCP in such a way that the function is a?ne over each cell of the decomposition. The algorithm is output-sensive in the sense that its time complexity is polynomial in the size of the input and linear in the size of the output, when the problem is non-degenerate. We give a lexicographic perturbation technique to resolve degeneracy as well. Unlike for the non-parametric case, the resolution turns out to be nontrivial, and in particular, it involves linear programming (LP) duality and multi-objective LP.

1

Introduction

Given a real square matrix M and a vector q , solving a linear complementarity problem (LCP) consists of ?nding two nonnegative vectors w and z that satisfy the conditions w ? M z = q, w ≥ 0, z ≥ 0, wT z = 0 . (1.1)

This simply stated and well-studied problem has far-reaching applications that have been welldocumented in the literature. Rather than give a survey here, the interested reader is referred to the books [22, 6]. Several authors have studied the properties of various parametric versions of this problem (e.g. [17, 2, 6, 19, 7, 26, 8]), but unless there are restrictions placed on the particular parametric LCP (pLCP) considered, it is in general unrealistic to expect an e?cient computational algorithm. We here study the class of pLCPs where the matrix M is su?cient1 and the right hand side (the vector q in (1.1)) is allowed to vary within a given a?ne subspace S . The goal is then to compute
1

Su?cient matrices are de?ned in Section 3.

1

functions z (·) and w(·) that map from the a?ne subspace S to a solution for pLCP (1.1) whenever one exists. This class of pLCP includes the important cases of linear and convex quadratic programs, where parameters appear linearly in the cost and the right hand side of the constraints [22]. In recent years, there has been a great deal of interest in the control community in parametric programming due to the fact that an important class of control algorithms for constrained linear systems, called model predictive controllers (MPC), can be posed as parametric linear or quadratic programs. The o?ine solution of these parametric problems results in an explicit representation of the optimal control action, which in some cases allows the controller to be implemented on systems with sampling rates of milli- and micro-seconds instead of the traditional seconds and minutes [25, 14, 4]. A similar setup results when computing optimal policies in a dynamic programming framework for partially observable Markov decision processes [18]. While parametric programming is widely used for sensitivity analysis, it is also applied in several other applications. In [15] it was shown that polyhedral projection can be reduced to parametric linear programming in polynomial time and of course such projections have uses ranging from the computation of invariant sets [5] and force closures [23] to program analysis [24] and theorem proving [13]. Polyhedral vertex and facet enumeration can also be posed as projection problems, and hence solved with the proposed pLCP approach [12], which as discussed below results in an output sensitive algorithm in the non-degenerate case (although not the most e?cient one for this purpose). In [2] it was shown that if M is a su?cient matrix and S satis?es certain general position 2 assumptions, then z (·) and w(·) are unique piecewise a?ne functions that are de?ned over a polyhedral partition of a convex set. There is, however, no known e?cient method of testing this general position assumption a priori and in fact, it is often not satis?ed even in the simplest case when the pLCP models a parametric linear program [16]. This paper extends the result of [2] by removing the restrictive and untestable general position assumption, allowing the algorithm to operate on any pLCP which is de?ned by a su?cient matrix and an a?ne subspace. This is achieved through a lexicographic perturbation technique, which has the e?ect of symbolically shifting the a?ne subspace an in?nitesimally small amount and into general position. We ?rst demonstrate that this perturbation always results in a problem that is in fact in general position and hence has the favorable uniqueness and partitioning properties discussed above. The challenge then becomes one of doing calculations in this perturbed space. The main optimization problem that arises as a result of the perturbation is a linear program that is polynomially parameterized by a positive variable ?. The decision problem to be tackled is then the determination of the behavior of this parametric problem as the parameter ? tends to zero. Section 5 discusses how this problem can be converted into a multi-objective linear program, which can then be solved e?ciently. The proposed technique should be applicable to other algorithms that rely on lexicographic perturbation to handle degeneracy. The resulting algorithm has the strong property that its complexity is polynomial in the size of the input (the matrix M ) and linear in the size of the output (the number of pieces in the piecewisea?ne functions w and z ). For this reason, we call the algorithm ‘output sensitive’, although it should be noted that the complexity of the functions w and z can be exponential in the worst case and that this complexity result is for the lexicographically shifted a?ne subspace, which may be more complex than the unshifted case. The reminder of the paper is organized as follows. Section 2 gives some basic notations and a formal de?nition of the parametric LCP. Section 3 provides some useful properties of pLCPs on su?cient matrices. Section 4 then presents the proposed method with a general position assumption,
2

General position is de?ned in Section 3.2.

2

and then this is relaxed in Section 5 where the lexicographic perturbation is introduced. Finally, Section 6 analyzes the complexity of the algorithm.

2

Parametric LCP, critical regions and their adjacency

Let us ?rst ?x some useful notations for matrices. For a matrix A ∈ Rm×n and a column index j ∈ {1, . . . , n}, A·j ∈ Rm denotes the j -th column vector of A. Similarly, for a row index i ∈ {1, . . . , m} Ai ∈ RI ×n denotes the i-th row vector of A. For a subset J ? {1, . . . , n}, A·J ∈ Rm×J denotes the matrix formed by the columns of A indexed by J , and for a vector v ∈ Rn , vJ denotes the vector formed by the components of v indexed by J . For I ? {1, . . . , m}, we denote by AI ∈ RI ×n the matrix formed by the rows of A indexed by I . Given a real square matrix M and a vector q of size n, the linear complementarity problem (LCP) is to ?nd two nonnegative vectors w and z that satisfy w ? M z = q, w ≥ 0, z ≥ 0, wT z = 0 . (2.2)

In this paper, we consider the LCP (2.2) where the right-hand side is allowed to vary within some a?ne subspace. Speci?cally, the goal is to ?nd two functions w(·) and z (·) that solve (2.2) over a given a?ne subspace S . De?nition 2.1. Let Q ∈ Rn×d be a matrix of rank d, q ∈ Rn a vector and M ∈ Rn×n a matrix of order n. The functions w(·) and z (·) are a solution to the pLCP (2.3) if for every θ ∈ Θf , w(θ ) and z (θ ) satisfy the relations w(θ ) ? M z (θ ) = q + Qθ , w(θ ), z (θ ) ≥ 0 , w(θ ) z (θ ) = 0 ,
T

(2.3a) (2.3b) (2.3c)

where Θf ? Rd is the set of feasible parameters θ , that is, those for which a solution to (2.3) exists. For the remainder of the paper we assume that the problem data M , q and Q are given and we de?ne A ∈ Rn×2n to be the matrix I ?M . Consider the following system of linear equality constraints in non-negative variables Ax = q , x≥0 . (2.4)

A basis is a set B ? {1, 2, . . . , 2n} such that |B | = n and rank(A·B ) = n; N := {1, . . . , 2n}\B is its complement and we call xB and xN the basic and non-basic variables respectively. Every basis B de?nes a basic solution to the linear system (2.4)
1 xB = A? ·B q ,

xN = 0 .

(2.5)

A basis B is called complementary if |{i, i + n} ∩ B | = 1 for all i = 1, . . . , n, and feasible if the 1 associated basic solution satis?es the nonnegativity constraint in (2.4), i.e. A? ·B q ≥ 0. Every T complementary feasible basis de?nes a solution of the LCP (2.2), by setting (w , z T ) = xT . In the parametric case, each basis is feasible for a set of parameters, which leads to the notion of a critical region. De?nition 2.2. The critical region RB of a complementary basis B is de?ned as the set of all parameter values for which B is feasible, i.e.,
1 RB := {θ ∈ Rd | A? ·B (q + Qθ ) ≥ 0} .

(2.6)

A complementary basis B is called feasible for the pLCP (2.3) if RB is nonempty. 3

By de?nition critical regions are convex polyhedra contained in the set of feasible parameters Θf . Each feasible complementary basis B de?nes a solution of the pLCP for each θ ∈ RB as w(θ ) w(θ ) 1 = 0, which is an a?ne function in RB . As a result, if Θf can = A? ·B (q + Qθ ) and z (θ ) z (θ ) B N be partitioned into a set of critical regions whose interiors are disjoint, then we have immediately a piecewise a?ne solution of pLCP (2.3) de?ned over these critical regions. In this paper we de?ne a set of conditions under which such a partitioning can be achieved and introduce an e?cient algorithm for this class of problems. The algorithm is based on the tracing of a graph whose nodes are the full-dimensional critical regions and whose edges are the pairs of adjacent regions (having a (d ? 1)-dimensional intersection). De?nition 2.3. Two critical regions R1 , R2 are called adjacent if their intersection R1 ∩ R2 is of dimension d ? 1. De?nition 2.4. Let V be the set of complementary bases B of pLCP (2.3) such that RB is fulldimensional and let E be the set of pairs of bases in V whose critical regions are adjacent. The graph G := (V, E ) is called the critical region graph of the pLCP (2.3). The proposed algorithm enumerates all full-dimensional critical regions by tracing the above graph. This tracing requires that we are able to enumerate all neighbors of a given complementary basis. The following section discusses the properties of this graph and investigates restrictions on matrices M under which the neighbor search can be done e?ciently.

3

Well behaving matrix classes for parametric LCPs

The goal of solving a parametric LCP is to compute functions w(·) and z (·) that satisfy (2.3) for all feasible values of the parameter θ ∈ Θf . As discussed in the introduction, linear complementarity problems include a very large set of di?cult optimization problems and so we cannot hope for a solution in the general case. In this section, we identify classes of LCPs that are ‘well-behaving’, or that have properties which guarantee that the algorithm given in Section 4 will ?nd a solution. The two key properties that will be needed are convexity of the feasible set Θf and the existence of a “canonical” single-valued mapping from parameters to critical regions. The latter essentially means that the relative interiors of critical regions do not intersect. In this section we will formalize these notions and discuss a well-known matrix class that has the appropriate properties when the a?ne subspace S ? Rn of all possible right hand sides is the whole space Rn . In Section 3.2 we will then generalize this and give conditions such that these properties still hold when the right hand side is restricted to lie in some lower-dimensional a?ne subspace.3

3.1

Complementary cones

We begin by describing the set of right hand sides q in (2.3) that are feasible for a given set of active constraints. For any index i ∈ {1, . . . , 2n} we denote with ? i the complementary index of i, i.e. ? i = (i + n) ? is de?ned as the set of all complementary indices of mod 2n. For a set I ? {1, . . . , 2n}, the set I elements in I . A set J ? {1, . . . , 2n} is called complementary if i ∈ J implies ? i ∈ J.
Throughout the paper, we use the same notation regarding matrix classes as in [6] and we use the properties of each class proved there. At the end of the paper we append an auxiliary section, where the relevant de?nitions and theorems are mentioned.
3

4

De?nition 3.1. For any complementary set J , the cone C (J ) := cone(A·J ) is called a complementary cone (relative to M ), where cone(T ) denotes the cone of all nonnegative combinations of the columns of a matrix T . If B is a complementary basis, then the complementary cone C (B ) is full-dimensional, and conversely if the complementary cone C (J ) is full-dimensional then the submatrix A·J has full rank, i.e. J is a complementary basis. For a complementary basis B , we have
1 C (B ) = {y ∈ Rn | A? ·B y ≥ 0} .

(3.7)

1 In the remainder of the paper we will denote by β the matrix A? ·B , where B is the considered basis. Therefore we will write C (B ) = {y ∈ Rn | βy ≥ 0}. One can see that for a given basis B , the cone C (B ) is the set of all right hand sides that are feasible for LCP (2.2). We are interested in LCPs that have complementary cones with disjoint interiors and so we introduce the class of su?cient matrices, which has this property.

De?nition 3.2. A matrix M ∈ Rn×n is called column su?cient if it satis?es the implication [zi (M z )i ≤ 0 for all i] =? [zi (M z )i = 0 for all i] . (3.8)

The matrix M is called row su?cient if its transpose is column su?cient. If M is both column and row su?cient, then it is called su?cient. Remark 3.3. We note that both positive semide?nite (abbreviated by PSD) and P-matrices are su?cient. For a given matrix M it is possible to test in ?nite time whether it is su?cient, although no polynomial time test is currently known. The class of LCPs with su?cient matrices has been studied extensively, partly because this class appears to capture all critical structures for LCPs to behave nicely. In particular, this class admits many fruitful results ranging from combinatorial algorithms and duality [11, 10] to the e?cient solvability by interior-point methods [20]. We will see that this class is ideal also for the investigation of parametric LCPs. We start with a key fact. Proposition 3.4 ([6, Theorem 6.6.6]). If M is a su?cient matrix, then the relative interiors of any two distinct complementary cones are disjoint. The union of all complementary cones forms a set known as the complementary range K (M ). The complementary range is equal to the set of all right hand sides of the LCP for which a feasible solution exists [6] K (M ) := {q | the LCP (2.2) with matrix M and right hand side q is feasible} . (3.9)

Proposition 3.5. If M is a su?cient matrix, then the complementary range K (M ) is a convex polyhedral cone K (M ) = cone([I ? M ]). Proof. The statement follows from the fact that su?cient matrices are in Q0 , see Theorem A.10. Remark 3.6. Throughout the paper we will draw upon the properties of two matrix classes extensively. The ?rst class is the Q0 -matrices, whose complementary range K (M ) is a convex cone and f the second is the fully semi-monotone matrices, denoted by E0 which have complementary cones f that are all disjoint in their interiors. The class of su?cient matrices is contained in Q0 ∩ E0 and is perhaps the largest known subclass de?ned by a simple set of conditions, which is why su?ciency is assumed for the majority of the results in this paper. It should be noted, however, that many of the results hold under slightly relaxed assumptions. 5

We will study now the adjacency relationship of complementary cones for the case of su?cient matrices. Speci?cally, since our goal is to compute the critical region graph G , ?nding all neighbors of any given region is a crucial issue. We ?rst look at the neighbors of a complementary cone that determine possible candidates for the neighbors for a critical region. De?nition 3.7. Two complementary bases B1 and B2 are called adjacent if their cones C (B1 ) and C (B2 ) are adjacent, that is, the dimension of C (B1 ) ∩ C (B2 ) is n ? 1. The following lemma is important in narrowing down the candidates of the neighbor search. Lemma 3.8. If M ∈ Rn×n is a su?cient matrix and B1 and B2 are adjacent complementary bases, then |B1 ∩ B2 | ≥ n ? 2. Proof. By the de?nition of adjacency the intersection C (B1 )∩C (B2 ) has dimension n ? 1 and therefore there exists a q ∈ K (M ) that lies in the relative interior of a facet of both complementary cones, which means that both basic solutions (w1 , z1 ), (w2 , z2 ) have exactly n ? 1 strictly positive components. Recall that basic solutions can be stated as: w1 z1
1 = A? ·B1 q , B1

w2 z2

1 = A? ·B2 q . B2

T , z T ) > 0 and |J | = n ? 1. By Theorem A.11, we have Let J be a subset of B1 such that (w1 1 J T T T T ? are (w2 , z2 )J ? = 0. Since (w2 , z2 )B2 has exactly one zero component, at least n ? 2 elements of J not in B2 and therefore their complements are. This shows |B1 ∩ B2 | ≥ n ? 2.

Remark 3.9. It can be shown for P-matrices that two bases are adjacent if and only if they di?er by exactly one element. This implies that the set of all complementary cones for a P-matrix LCP together with their faces forms a polyhedral complex. Unfortunately, this polyhedral complex property is not satis?ed in general for su?cient matrices, nor in fact for the proper subclass of PSD matrices. More precisely, the intersection of two critical regions may not be a common face, see [17]. Lemma 3.10. Let M ∈ Rn×n be su?cient and B be a complementary basis. If B ′ = B \{i} ∪ {? i} ′ is a basis then C (B ) and C (B ) intersect in their common facet C (B \{i}). Moreover no other fulldimensional complementary cones intersect the relative interior of C (B \{i}), i.e. C (B ′ ) is the unique complementary cone adjacent to C (B ) along this facet. Proof. Since B \{i} is a subset of B and B ′ , C (B \{i}) is a common facet of C (B ) and C (B ′ ). The second statement follows directly from the fact that the interior of any other complementary cone can intersect neither C (B ) nor C (B ′ ), since M is su?cient. Remark 3.11. Lemmas 3.8 and 3.10 imply that a complementary basis has at most n + adjacent complementary bases.
n2 ?n 2

Given a complementary basis B one can see that replacing any index i ∈ B with its complement ? i preserves complementarity, i.e. B \{i} ∪ {? i} is still a complementary set. This operation is called a diagonal pivot. If we substitute two di?erent indices i, j ∈ B with their complements, then the operation is called an exchange pivot. Lemma 3.8 ensures that for a given basis B we can reach all adjacent bases by a single diagonal pivot or by a single exchange pivot operation. However, for some i ∈ B the set B ∪ {? i}\{i} may not be a basis, or for some pair (i, j ) ∈ B the basis B ∪ {? i, ? j }\{i, j } may not be adjacent to B . Therefore, in order to determine whether a set given by a diagonal or an exchange pivot is in fact an adjacent feasible basis we need a further condition. Such a condition can be easily derived from the dictionary of the basis B . 6

1 De?nition 3.12. Given a complementary basis B and its complement N the matrix D := ?A? ·B A·N ∈ RB ×N is called the dictionary of B .

We begin by examining the diagonal pivot, for which a well-known adjacency condition can be derived. Fact 3.13. If B is a complementary basis, then for any i ∈ B , the set B \{i} ∪ {? i} is a basis if and only if Di,? = 0 . i We now consider the exchange pivot and derive necessary and su?cient conditions for adjacency, which are again based on examining elements of the dictionary. Proposition 3.14. Let M ∈ Rn×n be a su?cient matrix, B be a complementary basis and D be its dictionary. Consider the complementary basis B ′ = B \{i, j } ∪ {? i, ? j }, where i, j ∈ B are distinct. The following condition holds: dim(C (B \{i}) ∩ C (B ′ )) = n ? 1 ?? Di? i = 0 and Dj? i < 0. Proof. De?ne αk := ?Dk? i , then the following holds: A·? i =
k ∈B \{i}

(3.10)

αk A·k ,

(3.11)

Let j ∈ B \{i}, since αj = 0 (we have assumed B ′ to be a basis) we can rewrite (3.11) as ? ? 1 ? A·? αk A·k ? . A·j = i? αj
k ∈B \{i,j }

(3.12)

Let us consider q (λ) =
k ∈B \{i}

λk A·k ,

(3.13)

which lies in the relative interior of C (B \{i}) if and only if λk > 0 for all k. We can express q (λ) in following way by substituting (3.12) in (3.13): q (λ) =
k ∈B \{i,j }

λk A·k + λj (1/αj A·? i?
k ∈B \{i,j }

αk /αj A·k )

(3.14) (3.15)

=

(λk ? αk /αj )A·k + λj /αj A·? i
k ∈B \{i,j }

Su?ciency: if αj > 0 then there exists a q (λ) that lies in the relative interior of both facets C (B \{i}) and C (B ′ \{? j }). Necessity: since B ′ is a basis the unique way to express q (λ) as a linear combination of the vector indexed by B ′ is (3.15). If αj < 0 any q (λ) ∈ rel int(C (B \{i})) can not lie in C (B ′ ). The case αj = 0 is impossible since we have assumed B ′ to be a basis. Corollary 3.15 follows directly from the proposition above and allows the detection of the boundaries of the complementary range.
1 Corollary 3.15. Let M ∈ Rn×n be su?cient, B be a complementary basis and denote A? ·B as β . n Consider the facet C (B \{i}), for any i ∈ B . The hyperplane a?(C (B \{i})) = {y ∈ R | βi y = 0} de?nes a facet of the complementary range K (M ) if and only if Di? i = 0 and Di ≥ 0.

Remark 3.16. Lemma 3.8 and Proposition 3.14 are valid also for column su?cient matrices (i.e. it is not necessary that M is row su?cient). Lemma 3.10 holds for all fully semimonotone matrices. 7

3.2

Critical Domains

We study now the parametric case where the right hand side of LCP (2.2) is restricted to lie within some a?ne subspace S := {q + Qθ | θ ∈ Rd }. We will see, under some assumptions on S , that the properties of the complementary cones discussed in the previous section still hold in this case. De?nition 3.17. If B is a complementary basis, then the critical domain SB is the intersection of the a?ne subspace S with the complementary cone C (B )
1 d SB := C (B ) ∩ S = {y | A? ·B y ≥ 0, y = Qθ + q, θ ∈ R } .

(3.16)

Since we have assumed Q to be full column rank, the parametrisation q + Qθ is an invertible function and it is not hard to see that a critical domain is therefore the image of a critical region, i.e. SB = QRB + q . Since the parametrisation is a bijection, SB and RB have the same combinatorial structure for any complementary basis B . In particular, we have: Remark 3.18. The inequality βi y ≥ 0 is redundant in SB if and only if βi (Qθ + q ) ≥ 0 is redundant 1 in RB , where β = A? ·B . We now de?ne a key assumption, which will allow the extension of the properties of complementary cones to critical domains. De?nition 3.19. The a?ne subspace S is said to lie in general position if for every complementary basis B the following condition holds S intersects C (B ) ? S intersects int(C (B )) . (3.17)

If a critical domain has dimension d = dim(S ), we simply say that it is full-dimensional. By the de?nition above, we have: Remark 3.20. If S lies in a general position, then every critical domain is either full-dimensional or empty. Proposition 3.21. If M is su?cient and S lies in general position, then the relative interiors of critical domains SB1 and SB2 are disjoint for any two distinct complementary bases B1 and B2 . Proof. The statement is a direct consequence of Proposition 3.4, i.e. int(C (B1 )) and int(C (B2 )) are disjoint, and from rel int(SBi ) ? int(C (Bi )) for i = 1, 2. We denote the set of the feasible points of S by Sf . By (3.9) it follows that Sf = S ∩ K (M ). Corollary 3.22. If M ∈ Rn×n is su?cient then Sf is a convex polyhedron. Proof. The statement is a direct consequence of Proposition 3.5. If M is a su?cient matrix and S is in general position, then Proposition 3.21 and Corollary 3.22 ensure that the set K of nonempty critical domains de?nes a polyhedral decomposition of Sf in the sense that ? each member P of K is a convex polyhedron, ? ∪P ∈K P = Sf , ? dim P = d for all P ∈ K , and 8

? dim(P ∩ P ′ ) ≤ d ? 1 for any two distinct members P and P ′ of K . It is important to note that the set K may not induce a polyhedral complex, i.e. the intersection of two critical domains may not be a common face. Nevertheless, because Sf is convex, we can de?ne a graph structure of the decomposition which is connected. De?nition 3.23. Let V be the set of complementary bases B of pLCP (2.3) such that SB is fulldimensional and E consist of edges connecting each pair of bases in V whose critical domains are adjacent. The graph G := (V, E ) is called the critical domain graph of the pLCP (2.3). As stated above, each critical domain SB is the image of a critical region RB under the a?ne map θ → q + Qθ , and a similar statement can be make for the feasible sets Sf = q + QΘf . Since for each complementary basis B the critical domain SB and the critical region RB have the same combinatorial structure, the critical domain graph G = (V, E ) also de?nes the graph of critical regions and vice versa. In the discussion of the algorithm we will mostly consider only critical domains. Corollary 3.24. If M is a su?cient matrix and S lies in general position, then the graph of critical domains G is connected. Proof. The statement follows directly from the convexity of Sf = K (M ) ∩ S , which implies that between every pair of critical domains there exists a path in the graph of critical domains. In the previous section we have seen that for the case of su?cient matrices we can reach all adjacent cones from any complementary cone with a single diagonal or a single exchange pivot operation. Assuming general position of S , this useful property also holds for critical domains. Proposition 3.25. If M is su?cient and S lies in general position, then for any two complementary bases B1 and B2 that have nonempty critical domains the following holds: If SB1 and SB2 are adjacent in S then C (B1 ) and C (B2 ) are adjacent cones. Proof. If SB1 and SB2 are adjacent critical domains, then their intersection is contained in C (B1 ) ∩ C (B2 ). If C (B1 ) and C (B2 ) are not adjacent cones then there exists a complementary cone C adjacent to C (B1 ) that contains C (B1 ) ∩ C (B2 ). In this case C ∩ S would be adjacent to SB1 and would overlap with the relative interior of SB2 , which is a contradiction of Proposition 3.21. Given a complementary feasible basis B , Proposition 3.25 ensures that all the critical domains adjacent to SB can be reached by exploring complementary bases adjacent to B .

4

Description of the generic algorithm

Now we are able to present an algorithm that enumerates all complementary bases whose critical domains de?ne a polyhedral partition of Sf with the following two sets of assumptions: Assumption 4.1 (Regularity). The matrix M is su?cient and the matrix of the parametrisation Q ∈ Rn×d has full column rank. Assumption 4.2 (General Position). The a?ne subspace S := {q + Qθ | θ ∈ Rd } lies in general position with respect to the complementary cones relative to M . Assumption 4.1 is essential for our algorithm to work, whereas Assumption 4.2 will be relaxed in the next section where an extension of the algorithm simulating general position for any given a?ne subspace via a symbolic perturbation is presented. 9

The proposed algorithm given in Algorithm 1 is based on a standard graph search procedure. It assumes a given function neighbors(B ), which returns all bases whose critical domains are adjacent to that of a given basis B . The validity follows immediately from the connectivity of the critical domain graph, Corollary 3.24. As input it takes a matrix M and an a?ne subspace S that satisfy the above assumptions, as well as an initial feasible complementary basis B0 such that SB0 is fulldimensional. The basis B0 is ?agged as “unexplored” and added to the set of discovered bases B . In each iteration of the algorithm an unexplored basis B is selected from B , marked as “explored” and all bases that have adjacent critical domains are enumerated and added to B , marking the new bases as “unexplored”. Once all bases in B have been explored, then we have found all bases with full-dimensional critical domains. Algorithm 1 Enumerate all critical domains by graph search Input: A feasible basis B0 with dim(SB0 ) = d, a su?cient matrix M ∈ Rn×n and an a?ne subspace S that lies in general position. Output: The critical domain graph G = (B , E ).
1: 2: 3: 4: 5: 6: 7: 8: 9: 10:

Initialise the set of nodes B := {B0 } and edges E := ?. Flag B0 as “unexplored”. while there exists an unexplored basis B in B do Flag B as “explored” Bnew := neighbors(B ) //Bnew := neighbors? (B ) if S does not lie in general position Flag each B ′ ∈ Bnew \B as “unexplored” B := B ∪ Bnew E := E ∪ (B, B ′ ) for each B ′ ∈ Bnew end while return G = (B , E )

The remainder of this section describes how the results of the previous sections can be exploited to e?ciently enumerate all adjacent critical domains of a given basis, i.e., how the function neighbors(·) can be properly implemented. The following section will then detail how the method can be extended so that the general position assumption can be relaxed.

4.1

Neighborhood computation of a critical domain

This section details a computational method that enumerates all bases that de?ne adjacent critical domains of a given basis, i.e. how the basis is “explored”, under both Assumption 4.1 and Assumption 4.2. The function neighbors is given as Algorithm 2. Let B be a basis whose critical domain is fulldimensional. By Proposition 3.25, each adjacent critical domain must have a (d ? 1) dimensional intersection with a facet of SB . We begin therefore by ?rst computing all facets of SB and then by determining the critical domains that intersect each one. Given a complementary feasible basis B we determine which facets of C (B ) de?ne the facets of 1 SB by removing the redundant inequalities of SB = {y ∈ Rn | βy ≥ 0, y = Qθ + q }, where β := A? ·B . The hyperplane hi = {y | βi y = 0}, i ∈ B intersected with SB is a facet of SB if there exists a y ? ∈ SB such that βj y ? > 0 for all j ∈ B \{i} and βi y ? = 0. This fact relies on the general position

10

assumption, Assumption 4.2. Therefore hi ∩ SB is a facet of SB if and only if the following LP: t? = max s.t. t ?βj Qθ + t ≤ βj q, ?j ∈ B \{i} ?βi Qθ = βi q (4.18)

has an optimal value t? > 0 strictly positive. Algorithm 2 Function neighbors(B ): returns all bases whose critical domains are adjacent to SB . Input: A complementary basis B , the matrix M and the a?ne subspace S . M is assumed to be su?cient and S to lie in general position. Output: The set B of complementary bases, whose critical domains are adjacent to SB .
1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18:

B := ? 1 β := A? ·B 1 B ×N D := ?A? ·B A·N ∈ R for each i ∈ B do if βi y ≥ 0 is non-redundant in SB then if Di? i > 0 then add B \{i} ∪ {? i} to B else for each j ∈ B with Di? j < 0 do ′′ ? ? B := B \{i, j } ∪ {i, j } if dim(SB ∩ SB ′′ ) = d ? 1 then add B ′′ to B end if end for end if end if end for return B

//where N = {1, . . . , 2n}\B //by solving LP (4.18)

//by solving LP (4.19)

By solving LP (4.18) for each i ∈ B we can determine if C (B \{i}) de?nes a facet of SB or not; see Line 5 of Algorithm 2. If it does, then the goal is to determine which bases, if any, have critical domains that intersect this facet. From the previous section, we saw that there are three possible cases:
′ ′ ? Diagonal pivot. If Di? i > 0 then the cone C (B ) de?ned by B := B \{i} ∪ {i} is the unique complementary cone adjacent to C (B ) along the facet C (B ) ∩ hi , due to Lemma 3.10, and therefore SB ′ is the unique critical domain adjacent to SB along the facet SB ∩ hi . Since SB ′ is nonempty (as SB ∩ hi is included in SB ′ ) it is full-dimensional by Assumption 4.2.

Boundary of K (M ). C (B \{i}) is a facet of the complementary range and therefore no other complementary cones intersect it. From Corollary 3.15, this is the case when Di? i = 0 and Di ≥ 0. Exchange pivot. By looking at the dictionary D of B all complementary cones adjacent to C (B ) that contain the index ? i can be determined (see Theorem 3.14). However, not all such cones intersect SB with dimension d ? 1 and for this reason we need to test for each cone C (B ′′ ) adjacent to C (B ) 11

?M 1

I2 S

C1

C2

C3

S

I1

?M 2
(a) Two-dimensional slice of a three-dimensional example. The cones C2 and C3 are both adjacent to C1 along the same facet. However the a?ne space S does not intersect C2 and therefore S ∩ C2 is not adjacent to S ∩ C1 . This situation arises in the case of an exchange pivot and requires that adjacency much be checked by solving an LP. (b) Example of diagonal pivot: no adjacency check is needed.

Figure 1: Adjacency of critical domains. See Example 1 for details. whether SB ′′ is adjacent to SB , i.e. whether the condition dim(SB ∩ SB ′′ ) = d ? 1 holds. Now assume that the basis B ′′ := B \{i, j } ∪ {? i, ? j }, where j ∈ B \{i}, de?nes such an adjacent complementary cone according to Theorem 3.14. In order to determine whether SB ′′ is adjacent to SB we can solve following LP: max s.t. ?βk Qθ + t ′′ Qθ + t ?βk ?βi Qθ ′′ Qθ ?βj t ≤ βk q ?k ∈ B \{i} ≤ βk q ?k ∈ B ′′ \{j } = βi q ′′ q , = βj

(4.19)

?1 1 ′′ where β = A? ·B and β = A·B ′′ . As in the simpler case (4.18) above, the intersection SB ∩ SB ′′ has dimension d ? 1 if and only if the optimal value of (4.19) is strictly positive. In this case SB ′′ = C (B ′′ ) ∩ S is nonempty and by the general position assumption, it is also full-dimensional.

Remark 4.3. Since C (B ) and C (B ′′ ) are adjacent cones and the hyperplanes {y |βi y = 0} and ′′ y = 0} de?nes their shared facet, the two hyperplanes must be equivalent. Therefore in (4.19) {y |βj ′′ Qθ = β ′′ q or ?β Qθ = β q can be removed. one of the equality constraints ?βj i i j Example 1. In Figure 1(a) a two-dimensional slice of three, three-dimensional cones is shown. The cones C2 and C3 are both adjacent to C1 along the same facet. However the a?ne subspace S does not intersect C2 and therefore S ∩ C2 is not adjacent to S ∩ C1 . In Figure 1(b) we consider the complementary basis B1 = {1, 2} and the cone C (B1 ) = cone(I ) = {y ≥ 0}. The goal is to ?nd the critical domain that is adjacent to SB1 = C (B1 ) ∩ S . The inequality y1 ≥ 0 is not redundant in SB1 and therefore the hyperplane h1 = {y1 = 0} de?nes a facet of SB1 . Since cone(?M1 , I2 ) ∩ h1 is equal to C (B1 ) ∩ h1 , we have that cone(?M1 , I2 ) ∩ S is adjacent to SB1 . The other inequality y2 ≥ 0 is redundant in SB1 and therefore the critical domain cone(I1 , ?M2 ) ∩ S is not adjacent to SB1 .

12

5

Extension of the algorithm for S not in general position

The previous section presented an algorithm that enumerates all feasible bases and returns the graph of critical domains. The algorithm works only under the assumption that the image of the parametrisation S lies in general position; Assumption 4.2. However, this assumption is not realistic and it is highly desirable to remove it. In the case of degeneracy (i.e. S is not in general position), Propositions 3.21 and 3.25 are no longer valid, as can be seen in Example 2. Therefore, during neighborhood computation it is not su?cient to explore only the adjacent complementary cones. In order to extend the algorithm to the degenerate case, we apply a symbolic perturbation technique (the lexicographic perturbation) which has the e?ect of shifting S into general position. The next subsection will demonstrate how to handle the perturbation for neighborhood computation, in particular lines 5 and 11 of Algorithm 2. By using this technique we obtain a graph of critical domains G ? relative to the perturbed a?ne subspace S ? , which can di?er from the graph of critical domains G relative to S . In particular, some full-dimensional critical domains in S ? may be non full-dimensional in S . We will see that there exists a subgraph G of G ? that is a graph of critical domains relative to S and which can be obtained by postprocessing G ? . Example 2. This example demonstrates the e?ect when a parametric LCP is not in general position. Consider the parametric LCP de?ned by the matrices M= 1 ?1 1 , Q= 1 1 ?1 and q = 0 0 .

A ?gure depicting the complementary cones relative to M and of the a?ne subspace S is shown in Figure 5. Let B1 = {1, 2}, B2 = {2, 3}, B3 = {3, 4} and B4 = {1, 4}. For notational simplicity, we denote by Ci = C (Bi ) the complementary cones and by Si = S ∩ Ci for i = 1, . . . , 4 the critical domains. Clearly S does not lie in general position because it intersects C1 , C3 and C4 on their boundary but not in their interiors. Proposition 3.21 is violated because S1 is neither empty nor full-dimensional, and furthermore S3 and S4 are equal and hence rel int(S3 ) = rel int(S4 ). Theorem 3.25 is violated because S2 and S4 are adjacent, but C2 and C4 are not.

5.1

Lexicographic perturbation

This section presents a well-known method that permits the perturbation of the image of the parametrisation S into general position and which can be treated symbolically: the lexicographic perturbation. We introduce the following notation that will be used for the reminder of the paper. De?nition 5.1. The vector ? := (?, ?2 , ?3 , . . . , ?n )T is called the lexicographic perturbation vector and is a function of a positive real number ?. We denote with S ? the a?ne subspace S perturbed by ? S ? := S + ? = {y ∈ Rn | y = Qθ + q + ?, θ ∈ Rd } . (5.20)

Theorem 5.2. Let M be a su?cient matrix and S be the a?ne subspace {Qθ + q | θ ∈ Rd } for a given matrix Q and a vector q . There exists a δ > 0 such that S ? := S + ? lies in general position for each ? ∈ (0, δ).

13

S

I2

C C
2

1

C C
3

I1
4

?M1

?M2

Figure 2: An a?ne subspace S that is not in general position. (See Example 2) Remark 5.3. In the remainder of the paper we will use the standard expression “property A holds for all su?ciently small ? > 0” rather than the more cumbersome “there exists δ > 0 such that property A holds for each ? ∈ (0, δ)”. Therefore the claim of Theorem 5.2 can be written as: S ? lies in general position for all su?ciently small ? > 0. To prove Theorem 5.2, we need the following lemma. Lemma 5.4. Let Q ∈ Rn×d , q ∈ Rn and S ? = {Qθ + q + ?| θ ∈ Rd }. For any complementary cone C , there exist ?nitely many ? such that S ? intersects C but not int(C ). Proof. Let B be a complementary basis. We denote with hi the hyperplane hi := {y |βi y = 0} for all i ∈ B . Therefore C (B ) ∩ hi is a facet of C (B ) = {y ∈ Rn |βy ≥ 0}. We will prove that for any subset J ? B there are ?nitely many ? such that the following condition holds
? ? = SB = C (B ) ∩ S ? ? C (B ) ∩ ( i∈J ? hi ) and SB ? hi for all i ∈ J .

(5.21)

The statement of the lemma will then follow directly, since there are ?nitely many subsets of B . Let J be any nonempty subset of B . For any ? for which (5.21) holds, J contains the indices of all inequalities of S ? that are implicit equalities and it holds that βJ (Qθ + q + ?) = 0 for all θ ∈ R? B , (5.22)

? where R? B = {θ | Qθ + q + ? ∈ SB }. We can distinguish two cases. The case 1: there exists i ∈ J such that βi Q is a zero row vector. ? is nonempty, β (q + ?) = 0 for the condition (5.22) to be valid. This non-trivial polynomial Since SB i equation holds for at most n values of ?. Now consider the case 2: βJ Q has no zero rows. We will prove that the matrix βJ Q does not ? is nonempty and has dimension d ? rank(β Q) < d, the Chebyshev have full row rank. Since SB J

14

center problem has an optimal value of zero, i.e. max{t | βi Qθ ? t ≥ ?βi (q + ?) ?i} = 0. We consider its dual problem min (β (q + ?))T y s.t. ?(β Q)T y = 0 (5.23) yi = 1 y ≥ 0 . From strong duality, there exists a non-zero optimal solution y ? 0 such that (βQ)T y ? = 0 and (β (q + ?))T y ? = 0. We now claim that all indices of the strictly positive components of y ? are ? > 0. Then, for any θ ∈ R? , contained in J . Assume that there is an index i ∈ J with yi B ? βi (Qθ + q + ?) = 0, i.e. SB ? hi , which contradicts the maximality condition in (5.21). Therefore, ? ? = 0, i.e. there exists a non-trivial combination of rows of β Q. yJ 0 and (βQ)T yJ J Let {v1 , . . . , vs } be a set of columns of Q such that βJ v1 , . . ., βJ vs form a basis of the column space of βJ Q. Since βJ Q does not have full row rank, s < |J |. If βJ (q + ?), βJ v1 , . . . , βJ vs are linearly independent then for any θ ∈ Rd the equation (5.22) cannot hold. We claim these vectors are linearly dependent for ?nitely many ?. First we consider the case s = |J | ? 1. Then, these vectors are linearly dependent if and only if det(βJ (q + ?), βJ v1 , . . . , βJ vs ) = 0 . (5.24)

This condition is a polynomial equation in ? and holds for ?nitely many ?. Finally, if s < |J | ? 1, we use the same argument by adding a proper number of vectors v ?1 , . . . , v ?|J |?1?s such that βJ v1 , . . . , βJ vs , v ?1 , . . . , v ?|J |?1?s are |J | ? 1 linearly independent vectors. Proof of Theorem 5.2. We can assume without loss of generality that S does not lie in general position. For any complementary cone C exactly one of the following cases holds: 1. S does not intersect C , 2. S intersects the interior of C , 3. S intersects the boundary of C and S ∩ int(C ) = ? , which can be di?erentiated into two subcases: (a) ?δ > 0 such that C ∩ S ? = ? for all ? ∈ (0, δ) (b) For all δ > 0 there exists an ? ∈ (0, δ) such that C ∩ S ? = ? . For each of these cases we need to prove that there exists δ > 0 such that either S ? intersects int(C ) for all ? ∈ (0, δ) or C ∩ S ? = ? for all ? ∈ (0, δ). Clearly, this condition holds for cases 1, 2 and 3a. Therefore, it su?ces to prove that for case 3b there exists a δ > 0 such that S ? intersects int(C ) for all ? ∈ (0, δ). Let C be any complementary cone that satis?es condition 3b. From Lemma 5.4 and by assumption there exists a δ > 0 such that S δ intersects the interior of C and for any ? ∈ (0, δ) either int(C ) ∩ S ? = ? or C ∩ S ? = ?. More precisely, one can select any δ > 0 smaller than the smallest ? > 0 for which S ? intersects C but not int(C ). Since S ? shifts continuously with ?, there exists no ? ∈ (0, δ) such that C ∩ S ? = ?. Example 3. Figure 5.1 shows two examples in which S does not lie in general position. In the ?rst example (Figure 3(a)) the cones C ({2, 3}) = cone(I2 , ?M1 ) and C ({1, 4}) = cone(I1 , ?M2 ) contain adjacent critical domains, although they are not adjacent cones. In the second example (Figure 3(b)) two di?erent critical domains coincide. In higher dimensions the critical domains can overlap in 15

?M 1

I2

?M 1

I2

(ε 1 , ε2 )=ε I1 ?M 2 S ?M 2
(a)

I1 Sε I2

?M 1

I2

?M 1

S

(ε 1 , ε2 )=ε I1

Sε I1

?M 2

?M 2
(b)

Figure 3: Lexicographic perturbation of the a?ne subspace S several ways and therefore it is not evident how to choose an appropriate decomposition when this situation occurs. In both cases the a?ne subspace S can be arti?cially and symbolically shifted into general position through the use of lexicographic perturbation.

5.2

Neighborhood computation in S ?

Given a complementary basis B that is feasible in S ? , the goal is to determine the adjacent ? . Since S ? lies in general position, it su?ces to explore the adjacent bases of critical domains to SB ? at Line 5 of the basis B . Similarly to the non degenerate case, we ?rst determine the facets of SB Algorithm 3 and then compute the adjacent critical domains that intersect with each facet. ? = {y ∈ Rn | βy ≥ 0, y = Qθ + q + ?}, where Let B be a complementary basis and consider SB ?1 ? forms a facet of β := A·B . The hyperplane hi := {y | βi y = 0}, for some i ∈ B intersected with SB ? if there exists a y ? ∈ S ? such that β y ? > 0 for all j ∈ B \{i} and β y ? = 0. Therefore h ∩ S ? is SB j i i B B ? if and only if a facet of SB t? (?) = max s.t. t ?βj Qθ + t ≤ βj (q + ?), ?j ∈ B \{i} ?βi Qθ = βi (q + ?) (5.25)

has a positive optimal value for all ? > 0 su?ciently small. 16

?. Algorithm 3 function neighbors? (B ): returns all bases whose critical domains are adjacent to SB Input: A complementary basis B , a su?cient matrix M and an a?ne subspace S ? . ? . Output: The set of complementary bases B , whose critical domains are adjacent to SB

1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18:

B := ? 1 β := A? ·B 1 B ×N D := ?A? ·B A·N ∈ R for each i ∈ B do if isLexPositive(redundancy(B, i)) is true then if Di? i > 0 then add B ′ := B \{i} ∪ {? i} to B else for each j ∈ B with Di? j < 0 do B ′′ := B \{i, j } ∪ {? i, ? j} if isLexPositive(adjacency(B, B ′′ )) then add B ′′ to B end if end for end if end if end for return B

//where N = {1, . . . , 2n}\B //if βi y ≥ 0 is non-redundant

? ∩ S? ) = d ? 1 //if dim(SB B ′′

This decision problem is no longer an LP, because the right hand side of the constraints depends on a polynomial in ? and we want to know the behavior of t? (·) in the neighborhood of zero. In the next subsection we will propose an e?cient method for determining if t? (·) is positive for su?ciently small ?. ? , then we can distinguish the same three cases as If the hyperplane hi de?nes a facet of SB discussed in Section 4.1.
? is full-dimensional Diagonal pivot. If there is exactly one adjacent complementary basis B ′ , then SB ′ ? along S ? ∩ h . (by the general position of S ? ) and is the unique adjacent critical domain to SB i B

Boundary of K (M ). If there are no adjacent complementary cones to C (B ) along hi (see Corol? along the facet S ? ∩ h . lary 3.15), then there is no adjacent critical domain to SB i B Exchange pivot. If there are adjacent bases B ′′ with |B ′′ ∩ B | = 2 and dim(C (B ) ∩C (B ′′ )) = d ? 1, ? ∩ S ? ) = d ? 1 (Line 11 of Algorithm 3). then we must check for each such basis B ′′ whether dim(SB ′′ B ′′ Assume B := B \{i, j } ∪ {? i, ? j } where j ∈ B \{i}. As in the non-degenerate case, we formulate a decision problem similar to LP (4.19) to test the dimension of the intersection: t? (?) = max s.t. βk Qθ ? t βi Qθ ′′ Qθ ? t βk β? j Qθ t ≥ = ≥ = ?βk (q + ?), ?k ∈ B \{i} ?βi (q + ?) ′′ (q + ?), ?k ∈ B ′′ \{? j} ?βk ?β? ( q + ? ) , j

(5.26)

17

?1 1 ′′ ? ? where β := A? ·B and β := A·B ′′ . The two critical domains are adjacent, i.e. dim(SB ∩ SB ′′ ) = d ? 1, in S ? for all ? > 0 su?ciently small if and only if t? (?) > 0 for all su?ciently small ? > 0.

5.2.1

Symbolic computation of the parametric Chebyshev center problem

As seen in the previous subsection, the goal is to decide whether the optimal value t? (?) of (5.25) (and of (5.26)) is positive for all su?ciently small ? > 0. We call this decision problem a parametric Chebyshev center problem. Here we introduce a method that can compute exactly the behavior of t? (·) for su?ciently small positive ?; the proposed approach is summarized as Algorithm 4. The procedure is explained only for (5.25), since the same method can be easily applied for (5.26). The goal is to transform (5.25) into a multi-objective LP that can then be solved with any LP-solver. Recall the parametric LP (5.25): t? (?) := max s.t. t ?βj Qθ + t ≤ βj (q + ?) , for all j ∈ B \{i} ?βi Qθ = βi (q + ?) . (5.27)

For a ?xed value of ?, this is a linear program and its dual is: t? (?) = min (β (q + ?))T y s.t. ?(β Q)T y = 0 j =i y j = 1 yj ≥ 0 for all j ∈ B \{i} yi free . Let us denote its feasible region by F (i), that is, F (i) = {y ∈ Rn | ? (β Q)T y = 0,
j =i

(5.28)

yj = 1, yj ≥ 0 for all j ∈ B \{i}}.

(5.29)

In order to solve (5.28) symbolically we introduce the following standard notion. De?nition 5.5 (Lexico-positive). A vector a ∈ Rs is lexico-positive (denoted by a ? 0), if a = 0 and the ?rst non-zero component of a is strictly positive. Given two vectors x and y ∈ Rs , we write x ? y if and only if x ? y ? 0. A matrix is called lexico-positive if all its rows are lexico-positive. If S = {si }i∈I is a set of vectors, then sj is the lexico minimum of S if and only if si sj for each i ∈ I. The following theorem demonstrates that minimizing the polynomial cost function of (5.28) is equivalent to computing the lexicographic minimum of a vector. Theorem 5.6. If y is a feasible vector of the dual problem (5.28), then the two statements below are equivalent: 1. ?δ > 0 such that (β (q + ?))T y > 0 for all ? ∈ (0, δ) , 2. (β [q I ])T y ? 0 . Proof. The statement follows from the equality (β (q + ?))T y = (1, ?, ?2 , . . . , ?n )(β [q I ])T y , which holds for all y and for all ?. For every polynomial p(x) = (1, x, x2 , . . . , xn )(p0 , p1 , p2 , . . . , pn )T the following holds: there exists a δ > 0 such that p(x) > 0 for all x ∈ (0, δ) if and only if the ?rst non-zero coe?cient of (p0 , p1 , p2 , . . . , pn ) is positive, i.e. (p0 , p1 , p2 , . . . , pn ) is lexico positive. 18

We can now consider an equivalent problem that we call a lexicographic linear program (lexLP): redundancy(B, i) : T ? := lexmin (β [q, I ])T y s.t. y ∈ F (i) (5.30)

The cost of this optimization problem is vector valued and the operator lexmin means to compute the lexicographic minimum vector (β [q, I ])T y over all feasible decision variables y . We will denote ? , as the function this particular lexLP, which tests the redundancy of the i-th inequality in SB redundancy(B, i). Theorem 5.7. If t? (·) is the optimal value of (5.27) as a function of ? and T ? is the optimal value of (5.30), then the following holds: ?δ > 0 such that t? (?) > 0 for all ? ∈ (0, δ) ?? T ? ? 0 . Proof. The statement is a direct consequence of Theorem 5.6. Note that as is the case for linear programs, the restrictions and the objective function of (5.30) are linear, although the objective returns a vector instead of a scalar. We say that T ? is the optimal value of (5.30). If the vector βi Q is non-zero then the feasibility region is bounded and therefore the optimal value is always attained if the problem is feasible. For our purposes, it is not necessary to compute the entire vector T ? , but only a su?cient number of its elements in order to determine if it is lexico-positive or not. To this end, the goal is to ?nd the ?rst non-zero component of T ? and therefore the lex min problem (5.30) can be treated as a multi-objective LP in the following way. First (say at step 0) we solve the LP: T0 := min cT 0 y s.t. y ∈ F (i), (5.32) (5.31)

where c0 := β q . If T0 = 0 then we can conclude that the optimal value T ? of the problem (5.30) is lexico-positive or lexico-negative from the sign of T0 . Otherwise, if T0 does equal zero, then we must consider the next objective function c1 := β1 and minimise it while maintaining T0 = 0, and so on. If T0 = T1 = · · · = Tr?1 = 0, then at the step r we solve: Tr := min cT r y s.t. y ∈ F (i) T k = 0, . . . , r ? 1, ck y = 0, (5.33)

where c0 = β q and ck = βk for k = 1, . . . , n. If Tr = 0 is the ?rst non-zero value of T ? = (T0 , . . . , Tn ) then T ? of (5.30) is lexico-positive if Tr > 0 and lexico-negative otherwise. The resulting procedure is Algorithm 4, where the feasible region of the LP (5.33) is denoted by Fr .

19

Algorithm 4 Function isLexPositive(lexLP ) Input: A lex linear program lexLP as redundancy(·, ·) (5.30) or adjacency(·, ·) (5.34). Output: Answer about lexico-positiveness of the optimal value T ? of (5.30) or (5.34) respectively.
1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11:

Let c0 , c1 , . . . , cn be the objective functions of lexLP and F0 be the feasibility region of lexLP for r = 0 to n do T //by solving the LP (5.33) t? r := min{cr y | y ∈ Fr } ? if tr > 0 then return T ? is lexico-positive. else if t? r < 0 then return T ? is lexico-negative. else Fr+1 := Fr ∩ {y | cT k y = 0} end if end for

Remark 5.8. Note that the lexLP adjacency(B, i) always has a non-zero optimal value T ? because zero is not feasible in (5.33) and the optimal solution y ? of the last LP (5.33), with r = n, must be optimal also for all previous LPs. Since β has full rank we must have that βy ? is non-zero and therefore there must be at least one component of T ? that is non-zero. The parametric LP (5.26) that determines whether two critical domains are adjacent can also be solved using the same procedure. LP (5.26) can be rewritten as a lex min LP as follows: 0 1 adjacency(B, B ′′ ) : ? 0, k ∈ B \{i} ? ? ? ? 0 , k ∈ B ′′ \{? j} ? ? free. (5.34) We will call this lexLP adjacency(B, B ′′ ) because it tests the adjacency of SB and SB ′′ . Recall that one of the two equalities in (5.26) can be removed because one is redundant and therefore one of yi and xj is also redundant. Hence, (5.34) has the same structure as (5.30) and can be solved as a multi-objective LP as explained above (see Algorithm 4). ? ? T := lexmin ? ? ? ? s.t. ? ? ? (β [q, I ])T y + (β ′′ [q, I ])T x ?(β Q)T y ? (β ′′ Q)T x k ∈B ′′ \{? j } xk k ∈B \{i} yk + yk xk yi , x? j = = ≥ ≥

5.3

Post-processing of the graph of critical domains G ?

The previous section introduced a computational method for computing the graph of critical domains G ? relative to the lex-perturbed space S ? for all ? > 0 su?ciently small. The goal in this section is to recover the graph G = (V, E ) of critical domains relative to the original space S according to De?nition 3.23. The following theorems will show that one can construct from G ? the graph G of critical domains relative to the unperturbed space S . Theorem 5.9. Let M be a su?cient matrix, S an a?ne subspace and consider the graph of critical domains G? = (V ? , E ? ) relative to the lexicographically perturbed space S ? for su?ciently small ? > 0. Let V be the set of all bases B in V ? with full-dimensional critical domains SB . Then the following statements hold. 20

1. For each B ∈ V ? , the complementary cone C (B ) intersects S and therefore SB is nonempty. 2. For any two distinct bases B1 and B2 in V , SB1 and SB2 have disjoint relative interiors. 3. The set of all full-dimensional critical domains SB for B ∈ V covers Sf , i.e., SB = Sf = K (M ) ∩ S ,
B ∈V

and forms a polyhedral decomposition of Sf . Proof. First, note that complementary cones are closed and so the ?rst statement follows directly. To prove the second, note that for all ? > 0 su?ciently small (say ? < δ) S ? lies in general position. From Proposition 3.25 the critical domains de?ned by S ? are disjoint in their interiors for any positive ? < δ. The second statement follows from the fact that for any basis B , S ? ∩ C (B ) changes continuously in ?, for ? < δ. The third statement is proven in two steps. First we prove that the critical domains whose bases are in V ? de?ne a covering of S ∩ K (M ). Let q ∈ K (M ) ∩ S , since ? ∈ K (M ) and K (M ) is a convex cone, there exists a basis B such that q + ? ∈ C (B ) for all ? > 0 su?ciently small. Therefore q is in C (B ) and hence in SB because complementary cones are closed and thus B ∈V ? SB = Sf = K (M ) ∩ S . Since critical domains are closed and Sf is a convex polyhedron, the full-dimensional critical domains de?ne a covering of Sf . The above theorem demonstrates that the bases in V ? have nonempty critical domains in the unperturbed space S . Moreover, there exists a subset V whose critical domains form a polyhedral decomposition of Sf according to De?nition 3.23. The next theorem discusses how adjacency in G ? relates to adjacency in G . Theorem 5.10. Let M be a su?cient matrix and B1 , B2 be two complementary bases in V ? V ? , i.e. SB1 and SB2 both have dimension d. If SB1 and SB2 are adjacent, then there exists a path ? 1, B ? 2, . . . , B ? r in G ? = (V ? , E ? ) from B1 to B2 with the following property: S ? i intersects SB ∩ SB B 1 2 B with dimension d ? 1, i.e. dim(SB ? i ∩ SB1 ∩ SB2 ) = d ? 1 for all i = 1, . . . , r .
? and S ? are adjacent, (B , B ) is clearly the desidered path. We assume they are Proof. If SB 1 2 B2 1 not adjacent. We choose a q ? ∈ rel int(SB1 ∩ SB2 ) which is not contained in any critical domain ? ∈ Rd be the parameter with q ?. or in any face of dimension d ? 2 or less, and let θ ? = q + Qθ We look now (for a moment) at the parameter space and at the critical regions. The hyperplane f := {θ | aT θ = b} contains the intersecion of the two critical regions, i.e. f ? (RB1 ∩ RB2 ) and ? + ta. The image of θ (t) is q (t) = q consider its perpendicular (normal?) line θ (t) = θ ? + tQa in the ? original space S , respectively q (t) = q ? + ? + tQa in the perturbed space S ? . We know that for each ? > 0 su?ciently small every critical domain becomes either fulldimensional or empty. The full-dimensional ones vary continuously in function with ?. Consider ? and S ? for all a segment [q ? (t1 ), q ? (t2 )] of the line {q ? (t) | t ∈ R} such that it intersects either SB B2 1 ? > 0 su?ciently small. Because of the continuity no critical domain, which has dimension smaller than d ? 1 in the original space S intersects this segment for all ? > 0 su?ciently small. Similarly, ? of dimension smaller than d ? 1 intersects [q ? (t ), q ? (t )] for all ? > 0 for any B ∈ V ? no face of SB 1 2 su?ciently small. The desidered path is given by the critical domains which decompose the line segment between SB1 and SB2 .

Note that the last condition in the above theorem, along with Proposition 3.21, implies that dim(SB ? i ) = d ? 1 for i = 2, . . . , r ? 1 and therefore Theorems 5.9 and 5.10 imply the following corollary. 21

Corollary 5.11. Let M be a su?cient matrix, S an a?ne subspace and let G ? = (V ? , E ? ) be the graph of critical domains relative to S ? . Then, the graph of critical domains G = (V, E ) relative to S is related to G ? as follows: 1. V ? V ? 2. For every basis B ∈ V , SB has dimension d. 3. For each pair of bases B1 and B2 , the critical domains SB1 and SB2 are adjacent if and only ?1 , . . . , B ?r ) in G ? with B ? 1 = B1 , B ?r = B2 and dim(Sk ) = d ? 1 for if there exists a path (B ? k = 2, . . . , r ? 1 (or (B1 , B2 ) ∈ E ). The above corollary provides a simple procedure for computing a critical region graph G = (V, E ) relative to the unperturbed a?ne set S from the perturbed one G ? . We begin from the perturbed critical region graph G = (V, E ) := G ? and remove each node B from G that has a critical domain SB which is not full-dimensional and add all new edges (B1 , B2 ) to E satisfying the statement 3 of Corollary 5.11. From Theorem 5.9, the critical domains of the nodes of the resulting graph will form the desired polyhedral covering of the Sf . Theorem 5.10 states that the resulting graph contains edges for all adjacent bases, but may be overconnected since some of the critical domains SB of removed bases B may have had a dimension less than d ? 1. It remains, therefore, to test each edge in order to determine if the connected bases are in fact adjacent in the unperturbed space. As discussed previously, both operations for testing full-dimensionality and adjacency can be posed as linear programs. Remark 5.12. Note that much of the computation required to test for full-dimensionality of the critical domains for the unperturbed a?ne set has already been done while building the perturbed graph. Speci?cally, one can determine if a region is full-dimensional by examining the ?rst component T0 of the optimizer of LP (5.30).

6

Complexity of the algorithm

In this section we will discuss the complexity of the proposed algorithm, which enumerates all fulldimensional critical domains relative to the lexicographically perturbed a?ne subspace S ? . The well-known example by Murty (see [21] or see Chapter 6 in [22]), which was used to prove the nonpolynomiality of the Lemke and the principal pivoting methods, can be easily seen to demonstrate that the number of critical domains of a pLCP with an a?ne subspace S of dimension 1 and Pmatrix M is exponential in n. Since the complexity of the graph search (Algorithm 1) is a polynomial function of the number of critical domains, no algorithm for pLCP is polynomial in n. It is, however, possible to bound the number of operations required to explore the neighborhood of each critical domain, i.e. the complexity of Algorithms 2 and 3. Since each critical domain will be explored exactly once, we can say that the algorithm is output sensitive in that its complexity is a polynomial function of the number of full-dimensional critical domains and the size of input, provided that a polynomial-time algorithm for linear programming is used. We ?rst consider the general position case and study the complexity of Algorithm 2. Assume that M is of order n, the a?ne subspace S is of dimension d and let B be a complementary basis with a nonempty critical domain SB . The main computations of the function neighbors(B ) (Algorithm 2) are: ? Redundancy checking at Line 5 (Solve LP (5.28) with n variables and d + 1 constraints.) 22

? Checking adjacency for the case of an exchange pivot at Line 11 (Solve LP (4.19) with 2n variables and d + 1 constraints.) We denote the time necessary to solve an LP in standard form by TLP (var, eq), where var denotes the number of (nonnegative) variables and eq is the number of equality constraints. The time necessary to explore a critical domain can then be bounded as follows. Theorem 6.1. Let M ∈ Rn×n be a su?cient matrix and assume that the a?ne subspace S lies in general position. For each complementary basis B with a nonempty critical domain SB the time necessary to explore the neighborhood of SB is bounded by: n TLP (n, d + 1) + n2 ? n TLP (2n, d + 1). 2 (6.35)

Proof. Redundancy checking requires the solution of LP (5.28) once for each of the n inequalities of SB , which takes n TLP(n, d + 1) time. Adjacency checking by solving the LP (4.19) is necessary only in the case that the considered adjacent basis B ′′ di?ers by two elements from the basis B and 2 ?n 2 ?n since there are at most n 2 such bases, the second term n 2 TLP (2n, d + 1) follows. If S does not lie in general position, then the lexLPs (5.30) and (5.34) are solved instead of LPs (4.18) and (4.19). Each lexLP can be solved as a sequence of at most n + 1 LPs with the same variables and constraints (see Algorithm 4), which leads to the following complexity bound. Theorem 6.2. If M ∈ Rn×n is a su?cient matrix, then for each complementary basis B with ? the time necessary to explore the neighborhood of S ? can be bounded nonempty critical domain SB B by (n2 + n) TLP (n, d + 1) + n3 ? n TLP (2n, d + 1). 2 (6.36)

The above theorems bound the complexity of “exploring” a basis of the output in Algorithm 1 (Line 5). The condition at Line 6, which can be veri?ed in time bounded by the logarithm of the size of the output, ensures that each output basis is explored exactly once. As a result, the complexity of the algorithm grows linearly with the size of the output and so is output sensitive.

7

Example

In this section we present a simple illustrative example that arises from control theory. Consider the following discrete time constrained linear time-invariant system: x+ = 1 1 1 u , x+ 0.5 0 1 ||u||∞ ≤ 1 ,

||x||∞ ≤ 5,

where x ∈ R2 is the system state, x+ is the successor state and u ∈ R is the system input. A common method of control for this class of systems is Model Predictive Control, in which we solve

23

4

1
3 2

0.5 U* (x) 0 ?0.5 ?1 ?5 0
?4 ?3 ?2 ?1 0 x1 1 2 3 4 5

1 0 ?1 ?2 ?3 ?4 ?5 x2

1

5 x1

?2 x2

0

2

Figure 4: Polyhedral partition (left) and optimal input (right) for the control problem (7.37) at each point in time the following ?nite horizon optimal control problem:
N ?1

J ? (θ ) = min
k =0

Qxk

2 2

+ Ruk

2 2

+ Q f xN

2 2

subject to xk+1 = 1 1 1 u xk + 0.5 k 0 1 ||uk?1 ||∞ ≤ 1 , ?k ∈ {1, . . . , N }

||xk ||∞ ≤ 5 , x0 = θ ,

(7.37)

where θ is the current state of the system, the prediction horizon N is 5 and the weighting matrices Q, R and Qf are the identity. For high-speed systems, such as electric power converters, the goal is to solve the above quadratic program as rapidly as possible, in some cases at rates exceeding hundreds of kilohertz (e.g. [3]). By computing the optimizer o?ine as an explicit piecewise-a?ne function of the state θ , these speeds can be achieved [25, 14, 4]. The above parametric quadratic program is easily converted to a pLCP with a positive semi-de?nite matrix M [22], which was then solved using the proposed algorithm. The resulting polyhedral partition and mapping from the parameter θ to the optimizer u0 is shown in Figure 4.

8

Conclusion

In this paper an algorithm to enumerate all feasible bases of the parametric LCP de?ned by a su?cient matrix M and a lexicographically perturbed a?ne subspace S was proposed. It has been shown that the perturbed parametric LCP can be solved in a time linearly bounded by the size of the output and moreover, this output can be e?ciently post-processed in order to generate a polyhedral decomposition for the unperturbed original a?ne subspace S . One feature of the algorithm which is not ideal is the space requirement. Namely, the proposed algorithm must store all discovered feasible bases in the memory because it relies on the standard graph search technique. A great improvement can be made if we could apply the reverse search technique [1] which is essentially memory free. For this, it is necessary for the underlying graph to be oriented properly with exactly one sink. Somewhat similar to the present work, the paper 24

[9] proposed an algorithm to compute a polyhedral complex known as the Gr¨ obner fan which was shown to have such a “reverse search property.” Finding such an orientation for the graph of critical domains is an excellent subject of the future research.

References
[1] D. Avis and K. Fukuda. Reverse search for enumeration. Discrete Applied Mathematics, 65:21– 46, 1996. [2] G. Bank, J. Guddat, D. Klatte, B. Kummer, and D. Tammer. Non-linear Parametric Optimization. Akademie-Verlag Berlin, 1983. [3] A.G. Beccuti, G. Papafotiou, R. Frasca, and M. Morari. Explicit Hybrid Model Predictive Control of the dc-dc Boost Converter. In IEEE PESC, Orlando, Florida, USA, June 2007. [4] A. Bemporad, M. Morari, V. Dua, and E. N. Pistikopoulos. The explicit linear quadratic regulator for constrained systems. Automatica, 38(1):3–20, January 2002. [5] F. Blanchini. Set invariance in control - a survey. Automatica, 35(11):1747–1768, November 1999. [6] R.W. Cottle, J.-S. Pang, and R.E. Stone. The Linear Complementarity Problem. Academic Press, 1992. [7] R.A. Danao. On the parametric linear complementarity problem. Journal of Optimization Theory and Applications, 95(2):445–454, November 1997. [8] B. De Schutter and B. De Moor. The extended linear complementarity problem. Mathematical Programming, 71(3):289–325, December 1995. [9] K. Fukuda, A. Jensen, and R. Thomas. Computing Gr¨ obner fans. Mathematics of Computation, 76:2189–2212, 2007. [10] K. Fukuda, M. Namiki, and A. Tamura. EP theorems and linear complementarity problems. Discrete Applied Mathematics, 84:107–119, 1998. [11] K. Fukuda and T. Terlaky. Linear complementarity and oriented matroids. Journal of the Operations Research Society of Japan, 35:45–61, 1992. [12] J.E. Goodman and J. O’Rourke, editors. Handbook of discrete and computational geometry. CRC Press, Inc., Boca Raton, FL, USA, 1997. [13] J. Hooker. Logical inference and polyhedral projection. pages 184–200. 1992. [14] T. A. Johansen, I. Petersen, and O. Slupphaug. On explicit suboptimal LQR with state and input constraints. In Decision and Control, 2000. Proceedings of the 39th IEEE Conference on, volume 1, pages 662–667 vol.1, 2000. [15] C. N. Jones, E. C. Kerrigan, and J. M. Maciejowski. On polyhedral projection and parametric programming. Journal of Optimization Theory and Applications, 137(3), June 2008. [16] C.N. Jones, J.M. Maciejowski, and E.C. Kerrigan. Lexicographic perturbation for multiparametric linear programming with applications to control. Automatica, 43(10):1808–1816, October 2007. 25

[17] C.N. Jones and M. Morari. Multiparametric Linear Complementarity Problems. In IEEE Conference on Decision and Control, December 2006. [18] L.P. Kaelbling, M.L. Littman, and A.R. Cassandra. Planning and acting in partially observable stochastic domains. Arti?cial Intelligence, 101:99–134, 1998. [19] D. Klatte. On the Lipschitz behavior of optimal solutions in parametric problems of quadratic optimization and linear complementarity. Optimization, 16:819–831, 1985. [20] M. Kojima, N. Meggiddo, T. Noma, and A. Yoshise. A uni?ed approach to interior point algorithms for linear complementarity problems, volume 538 of Lecture Notes in Computer Science. Springer-Verlag, 1991. [21] K. G. Murty. Computational complexity of complementary pivot methods. Mathematical Programming Study 7, pages 61–73, 1978. [22] K. G. Murty and F. T. Yu. Linear Complementarity, Linear and Nonlinear Programming. Helderman-Verlag, 1988. [23] J. Ponce, S. Sullivan, A. Sudsang, J. Boissonnat, and J. Merlet. On computing four-?nger equilibrium and force-closure grasps of polyhedral objects. International Journal of Robotics Research, February 1995. [24] S. Sankaranarayanan, F. Ivanˇ ci? c, and A. Gupta. Program analysis using symbolic ranges. pages 366–383. 2007. [25] M. M. Seron, G. C. Goodwin, and J. A. De Don? a. Geometry of model predictive control for constrained linear systems. Technical Report EE0031, The University of Newcastle, Australia, 2000. [26] K. Tammer. Parametric linear complementarity problems. 29T15:56:22Z]. [Online: Stand 2008-04-

A

Useful properties of matrix classes

This section gives an overview of some matrix classes with important properties for linear complementarity problems. The reader is referred to [6] for a thorough survey.

P-Matrices
De?nition A.1. The matrix M ∈ Rn×n is a P-matrix if and only if all principal minors of M are strictly positive. This class characterizes the matrices M for which the corresponding LCP always has a unique solution. Theorem A.2. The following statements are equivalent: 1. M ∈ P, 2. The LCP de?ned by the matrix M has a unique solution for all right hand side vectors q ∈ Rn ,

26

3. M does not reverse the sign of any nonzero vectors, i.e. [zi (M z )i ≤ 0 for all i] ? [z = 0] . Recall that M ∈ Rn×n is a positive de?nite matrix, denoted with P D if for all x ∈ Rn it holds that xT M x > 0. It is then easy to see from the above theorem that positive de?nite matrices belong to the class P.

P0 -matrices
De?nition A.3. The matrix M ∈ Rn×n is a P0 -matrix if and only if all principal minors of M are non-negative. Analogously to the positive-de?nite case above, positive-semide?nite matrices (PSD) are clearly in P0 . The following theorem gives properties of PSD matrices relevant to the solution of pLCPs. Theorem A.4. Let M ∈ Rn×n be a matrix and I be the identity matrix of the same order. The following statements are equivalent: 1. M ∈ P0 , 2. For each vector x = 0 there exists an index k such that zk = 0 and zk (M x)k ≥ 0, 3. (M + ?I ) is a P matrix for all ? > 0.

Semimonotone matrices
De?nition A.5. A matrix M ∈ Rn×n is called semimonotone if the following holds: for all x ≥ 0, x = 0 ? [xk > 0 and (M x)k ≥ 0 for some k] . (1.38)

The class of such matrices is denoted by E0 and by Theorem A.4, every P0 -matrix is semimonotone. De?nition A.6. Let M ∈ Rn×n be a semimonotone matrix. If for all index subsets α ? {1, . . . , n} with det(Mαα ) = 0 the principal pivot transform of M with respect to α M ′ :=
?1 ?1 M Mαα ?Mαα αα ? ?1 ?1 ?1 Mαα ? Mαα Mα ? Mαα Mαα ? ?α ? ? Mαα

is semimonotone, then M is called fully semimonotone. The class of such matrices is denoted with f f E0 and M is said to be an E0 -matrix.

Q0 -Matrices
De?nition A.7. An LCP de?ned by the matrix M and right hand side vector q is called weakly feasible if there exist positive vectors z and w such that w ? M z = q and feasible if z ′ w = 0 also holds. The class of matrices M for which the LCP is feasible whenever it is weakly feasible, is denoted by Q0 . Since P-matrices are feasible for each vector q , they are also Q0 -matrices. Theorem A.8. Let M ∈ Rn×n and I be the identity matrix of same order. The following statements are equivalent: 27

1. M ∈ Q0 , 2. The complementary range K (M ) is convex, 3. K (M ) = cone([I ?M ]) The implications of convexity of the complementary range are discussed in the next section.

Su?cient Matrices
De?nition A.9. A square matrix M is called column su?cient if it satis?es the implication: [zi (M z )i ≤ 0 for all i] =? [zi (M z )i = 0 for all i] . (1.39)

The matrix M is called row su?cient if its transpose is column su?cient. If M is both column and row su?cient, then it is said to be su?cient. Theorem A.10. If M is row su?cient matrix, then 1. M ∈ P0 , 2. M ∈ Q0 . From the theorem above we have that every column su?cient matrix also belongs to P0 . Below we state a characterisation of column su?cient matrices, which has an important implication regarding the structure of the resulting complementary cones. Theorem A.11. Given a matrix M ∈ Rn×n , the following statements are equivalent: 1. M is column su?cient, 2. For each vector q ∈ Rn the following holds: if (z 1 , w1 ), (z 2 , w2 ) are two solutions of the LCP de?ned by the matrix M and the vector q , then (z 1 )T w2 = (z 2 )T w1 = 0.

28


赞助商链接
更多相关文档:
更多相关标签:
网站地图

文档资料共享网 nexoncn.com copyright ©right 2010-2020。
文档资料共享网内容来自网络,如有侵犯请联系客服。email:zhit325@126.com