AN INTERIOR-POINT APPROACH TO SENSITIVITY ANALYSIS IN DEGENERATE LINEAR PROGRAMS ∗ E. ALPER YILDIRIM† AND MICHAEL J. TODD‡ Abstract. We consider the interior-point approach to sensitivity analysis in linear programming developed by the authors. We investigate the quality of the interior-point bounds under degeneracy. In the case of a special degeneracy, we show that these bounds have the same nice asymptotic relationship with the optimal partition bounds as in the nondegenerate case. We prove a weaker relationship for general degenerate linear programs. Key words. sensitivity analysis, degeneracy, interior-point methods, linear programming AMS subject classifications. 90C31,90C51,90C05 1. Introduction. Sensitivity analysis (or post-optimality analysis) is the study of how the optimal solution of an optimization problem changes with respect to the changes in the problem data. The possible presence of errors in the problem data often makes sensitivity analysis as important as solving the original problem itself. In the context of linear programming, sensitivity analysis can be performed using an optimal basis approach (as in the simplex method) or an optimal partition approach, where the optimal partition refers to knowing, for each index, whether the corresponding component of an optimal primal solution or of an optimal dual slack vector can be positive. The latter approach has close connections with interior-point methods since such methods, when properly terminated, provide an optimal solution in the relative interior of the optimal face, from which the optimal partition is readily available. In fact, as will shortly be discussed in more detail, the optimal partition approach has been developed by Adler and Monteiro [1] and Jansen, de Jong, Roos and Terlaky [7] as a promising alternative in order to circumvent the drawbacks of the classical optimal basis approach in the presence of degeneracy. Later, Monteiro and Mehrotra [9] extended this approach by relaxing the requirement that the optimal partition be known. They also provided two methods to estimate the range of perturbations, each of which can be performed at any optimal solution, regardless of where it lies on the optimal face. More recently, Greenberg, Holder, Roos and Terlaky [5] related the dimension of the optimal set to the dimension of the set of objective perturbations for which the optimal partition is invariant. Greenberg [4] considered the simultaneous perturbations of the right-hand side and the cost vectors from an optimal partition perspective. In [13], the authors studied perturbations of the right-hand side and the cost parameters in linear programming, motivated by how interior-point methods from a near-optimal pair of strictly feasible solutions for a problem and its dual would compare with the optimal basis approach obtained from a nondegenerate optimal basic solution for such perturbations. The proposed interior-point perspective stems from the objectives of regaining feasibility and maintaining near-optimality in a single it- ∗ Copyright (C) by SIAM. SIAM Journal on Optimization 12 (2002), 692–714. This research was supported in part by NSF through grant DMS-9805602 and ONR through grant N00014-96-1-0050. †Department of Applied Mathematics and Statistics, State University of New York at Stony Brook, Stony Brook, New York 11794, USA (yildirim@ams.sunysb.edu). This research was performed as part of the author’s Ph.D. study at Cornell University. ‡School of Operations Research and Industrial Engineering, Cornell University, Ithaca, New York 14853, USA (miketodd@cs.cornell.edu). 1 2 E. ALPER YILDIRIM AND MICHAEL J. TODD eration of the interior-point method. This requires the setup of the “right” Newton system among many possible choices in order to achieve both objectives simultaneously. Such a perspective provides a basis for the comparison of the interior-point and the simplex approaches to sensitivity analysis. Under the assumption of a unique, nondegenerate optimal solution, the authors showed that the Newton system proposed in [13] is the “right” one in the sense that it yields asymptotically the same bounds on perturbations as those that keep the current basis optimal (after symmetrization with respect to the origin). Similar results, but changing only one of the primal or dual near-optimal solutions, were obtained by Kim, Park and Park [8]. However, most linear programs (LPs) arising from real-life problems are degenerate. Our goal in this paper is to investigate the quality of the bounds from the interior-point perspective in the absence of the strong assumption of nondegeneracy. This will lead to a complete analysis of the interior-point perspective proposed in [13]. In doing so, we need something to compare our interior-point bounds with. In contrast to the nondegenerate case, the presence of multiple optimal bases makes a simplexbased approach unsuitable, as will be explained shortly. We therefore compare our bounds to those obtained from considering how much the right-hand side or the cost vector can change while maintaining the same optimal partition. Consequently, we use different tools for our analysis in this paper. The next section is devoted to preliminaries including the introduction of the tools relevant for the analysis as well as a restatement of our interior-point approach. Section 3 discusses the equivalence between the primal and dual formulations and shows that it suffices to consider perturbations of the right-hand side only. We analyze the interior-point bounds under a special case of degeneracy in Section 4 and extend the analysis to the general degenerate case in Section 5. We apply our interior-point approach to a well-documented, degenerate transportation example in Section 6 and Section 7 concludes the paper. 2. Preliminaries. We consider the LP in the following standard form: min cT x, subject to Ax = b, x ≥ 0. (P) x The associated dual LP is given by max bT y, subject to AT y + s = c, s ≥ 0. (D) y,s Here, A ∈ Rm×n, b ∈ Rm and c ∈ Rn constitute the data, and (x, y, s) ∈ Rn ×Rm ×Rn are the decision variables. Throughout this paper, the coefficient matrix A will be fixed and we will consider one-dimensional perturbations of the right-hand side vector b and the cost vector c, i.e., b will be replaced by b + t∆b and c by c + t∆c, where ∆b and ∆c will be fixed in Rm and Rn, respectively, and t ∈ R will be the parameter. This is also called parametric analysis in the literature. We will make the following assumptions: Assumption 2.1. The coefficient matrix A has full row rank. Assumption 2.2. Both (P) and (D) have strictly feasible solutions, i.e., there exist x > 0, s > 0 and y such that Ax = b and AT y + s = c. While Assumption 2.1 is without loss of generality, Assumption 2.2 is clearly restrictive. However, we will discuss how our approach can be extended to LPs which do not satisfy Assumption 2.2 (but do have optimal solutions) at the end of Section 2.2. SENSITIVITY IN DEGENERATE LINEAR PROGRAMS 3 The classical approach to sensitivity analysis has been based on the simplex method. Assuming that an optimal solution exists, the simplex method terminates with a basic optimal solution along with a corresponding basis. A natural criterion for the allowable perturbations in the data is then given by the range of such perturbations for which the current basis remains optimal for the resulting family of LPs. Let us consider the parametric right-hand side (RHS) problem, i.e., let b be replaced by b + t∆b. Define v(t) = min{cT x : Ax = b + t∆b, x ≥ 0}. It is well-known that v is a convex, piecewise linear, continuous function of t. The parametric RHS problem includes finding out all the “breakpoints” of v(t). Fixing a value of t, say at 0 for the purposes of this paper, the classical approach to sensitivity analysis then provides the set of values of t for which an optimal basis for t = 0 remains optimal for the resulting LPs parametrized by t. This is called the optimality interval associated with an optimal basis. Note that the optimal basis approach indeed yields the breakpoints of v(t) around 0 under primal and dual nondegeneracy (which holds only if 0 itself is not a breakpoint of v(t)). However, the presence of primal and/or dual degeneracies is a shortcoming for this approach since, for example, multiple optimal bases might yield different optimality intervals. This shortcoming has been observed by several researchers. Adler and Monteiro [1], and Jansen, de Jong, Roos and Terlaky [7] developed an optimal partition approach to sensitivity analysis and showed that the optimality intervals associated with the optimal partitions uniquely and unambiguously identify the breakpoints of v(t) and the intervals between the consecutive breakpoints. By the symmetry between (P) and (D), which will be treated in more detail in Section 3, the same conclusions also hold for the parametric analysis of the cost vector c. The idea of the optimal partition is based on a well-known result of Goldman and Tucker [2]. The optimality conditions for (P) and (D) are given by primal and dual feasibility and complementary slackness, that is, a triple (x, y, s) is optimal for (P) and (D) if and only if it satisfies (2.1) Ax = b, AT y + s = c, xisi = 0, i = 1, . . . , n, x ≥ 0, s ≥ 0, where xi and si denote the ith components of x and s, respectively. Let ΩP and ΩD denote the set of optimal solutions for (P) and (D), respectively. Then, we can define two index sets as (2.2) B = {j ∈ {1, . . . , n} : xj > 0 for some x ∈ ΩP }, N = {j ∈ {1, . . . , n} : sj > 0 for some (y, s) ∈ ΩD}. The optimality conditions (2.1) imply that B ∩ N = ∅. The Goldman-Tucker result indicates that B and N actually partition the index set {1, . . . , n}, i.e., B ∪ N = {1, . . . , n}. Therefore, there exist at least one primal optimal solution x ∈ ΩP and one dual optimal solution (y, s) ∈ ΩD such that x + s > 0. Such a solution will be called strictly complementary and (B, N ) will be called the optimal partition. In contrast to the possibility of multiple optimal bases, the optimal partition is unique for a given LP instance. We will denote by B and N the columns of A corresponding to the indices in B and N , respectively, and we will also partition the cost vector c as cB and cN , and the variables x and s as xB and xN , and sB and sN accordingly. Note that if (x, y, s) is a strictly complementary solution, then we have xB > 0, xN = 0, sB = 0 and sN > 0. 4 E. ALPER YILDIRIM AND MICHAEL J. TODD Let us again restrict our attention to one-dimensional perturbations of the righthand side vector b. The optimal partition approach is based on maintaining the whole dual optimal set invariant rather than an optimal basis as in the classical simplex approach. Note that perturbations of b do not affect the dual feasible region. Consequently, the range of t is given by solving two auxiliary LPs. More precisely, if b is replaced by b + t∆b, and if ΩD denotes the dual optimal set for (D) (i.e., t = 0), then the lower and upper bounds on t are given by the optimal values of (AUX1) minx,λ (maxx,λ) λ subject to Ax = b + λ∆b, x ≥ 0, (s∗)T x = 0, ∀ (y∗, s∗) ∈ ΩD. We will call the resulting bounds the optimal partition bounds. Note that both prob- lems are always feasible since λ = 0 together with any x ∈ ΩP satisfy all the con- straints. Fixing the dual optimal set ΩD is equivalent to fixing the optimal partition (B, N ) by the Goldman-Tucker result. Therefore, the (possibly infinite) last constraint set in (AUX1) can be replaced by the equivalent single constraint xT s∗ = 0, where s∗ is any point in the relative interior of ΩD (hence s∗N > 0). This condition, in turn, is the same as setting xN = 0. Consequently, (AUX1) can be written in the following simplified form: (AUX1) minxB,λ (maxxB,λ) subject to λ BxB = b + λ∆b, xB ≥ 0. The analogous derivation for the one-dimensional perturbations of the cost vector c leads to the following auxiliary problems, whose optimal values give the optimal partition bounds for t when c is replaced by c + t∆c: (AUX2) miny,sN ,λ (maxy,sN ,λ) λ subject to BT y = cB + λ∆cB , N T y + sN = cN + λ∆cN , sN ≥ 0. Here, ∆cB and ∆cN constitute the corresponding partition of ∆c. Before getting into the symmetrized bounds we would like to recall an important result about the dimensions of the optimal solution sets ΩP and ΩD. In what follows, dim(·) denotes the dimension and | · | denotes the cardinality of a set. The reader is referred to Lemma IV.44 in [10] for a proof. Proposition 2.1. dim (ΩP ) = |B| − rank (B) ; dim (ΩD) = m − rank (B). 2.1. Symmetrized Bounds. The auxiliary problems (AUX1) and (AUX2) can be reformulated in the following way. Let us consider (AUX1) and let x∗ ∈ ΩP . Then, the equality constraint can be rewritten as BxB = Bx∗B + λ∆b or B(xB − x∗B ) = λ∆b. SENSITIVITY IN DEGENERATE LINEAR PROGRAMS 5 Therefore, by a change of variable, if we let u = xB − x∗B, then (AUX1) is equivalent to (AUX1) minu,λ (maxu,λ) λ subject to Bu = λ∆b, u ≥ −x∗B. Next, we will tighten the constraints in the above formulation by putting upper bounds on u as well, and our choice for the upper bound will be x∗B, which will give the largest symmetric L∞-box around the origin which is contained in the feasible region: (SA1) minu,λ (maxu,λ) λ subject to Bu = λ∆b, −x∗B ≤ u ≤ x∗B . We will call (SA1) the symmetrized LP and the resulting optimal solutions the symmetrized bounds. The formulation of (SA1) reveals that if (u∗, λ∗) solves the maximization problem, then (−u∗, −λ∗) solves the minimization problem. Therefore, it suffices to solve one LP as opposed to solving two LPs to obtain the optimal parti- tion bounds from (AUX1). A similar treatment of (AUX2) gives rise to the following symmetrized LP: (SA2) minv,w,λ (maxv,w,λ) λ subject to BT v = λ∆cB , N T v + w = λ∆cN , −s∗N ≤ w ≤ s∗N , which is obtained by replacing y − y∗ by v and sN − s∗N by w, where (y∗, s∗) ∈ ΩD. Finally, a similar symmetrization has been applied to w. Next, we would like to discuss the relationship between the auxiliary and the symmetrized LPs. First of all, let us assume that both (P) and (D) have unique and nondegenerate optimal solutions. Then, Proposition 2.1 implies that B is actually a square and nonsingular matrix, hence invertible. In fact, B is the optimal basis. Consequently, (AUX1) and (AUX2) are trivial to solve and their optimal solutions coincide with the optimal basis bounds arising from the simplex method. With this observation, the constraints of (AUX1) reduce to (2.3) λB−1∆b ≥ −x∗B or λ(XB∗ )−1B−1∆b ≥ −e, where XB∗ is the diagonal matrix whose components are given by x∗B and e denotes the vector of ones in the appropriate dimension. Similarly, the constraints of (SA1) can be rewritten as (2.4) −e ≤ λ(XB∗ )−1B−1∆b ≤ e or |λ| (XB∗ )−1B−1∆b ≤ 1, ∞ where · ∞ is the L∞-norm. A similar treatment reveals that the constraints of (AUX2) are equivalent to (2.5) λ(SN∗ )−1(∆cN − N T B−T ∆cB ) ≥ −e, 6 E. ALPER YILDIRIM AND MICHAEL J. TODD where SN∗ is defined similarly, and that those of (SA2) to (2.6) |λ| (SN∗ )−1(∆cN − N T B−T ∆cB) ∞ ≤ 1. The derivations (2.3)–(2.6) imply the following relationship between the auxiliary and the symmetrized LPs: let (λ−, λ+) denote the optimal partition bounds given by the optimal solutions of the auxiliary LPs (including possibly ±∞). Then, the symmetrized bounds for t are (−λs, λs), where (2.7) λs = min (|λ−|, λ+). Therefore, the symmetrized bounds are indeed equal to the “symmetrization” of the optimal partition bounds. Note that we used an open interval for the optimal partition bounds above. The reason is that the optimal partition remains the same in this open interval whereas it does change at the left and right limit points (assuming they are finite). Next, let us assume that (P) has a unique but degenerate optimal solution. Then, by Proposition 2.1, B is nonsquare but it has full column rank. Therefore, (AUX1) is still easy to solve. If ∆b does not lie in the range space of B, then the optimal solutions of (AUX1) and (SA1) are all zero (which implies that t = 0 is a breakpoint of v(t)). Otherwise, there exists a unique vector v such that Bv = ∆b, and hence, the constraints of (AUX1) are equivalent to (2.8) λ(XB∗ )−1v ≥ −e. Similarly, the constraints of (SA1) can be stated as (2.9) |λ| (XB∗ )−1v ∞ ≤ 1. Once again, we conclude that a similar symmetry as in (2.7) continues to hold between (SA1) and (AUX1). In a similar manner, one can show that such a relationship holds between (SA2) and (AUX2) if (D) has a unique but degenerate optimal solution. The preceding discussion shows that the optimal solutions of the auxiliary and the symmetrized LPs have the nice relationship (2.7) as long as there is a unique optimal solution that one can use to symmetrize the constraints of the auxiliary LPs to obtain the symmetrized LPs. An interesting question then is whether the same nice relationship continues to hold between the optimal partition bounds and the symmetrized bounds if there are multiple optimal solutions, that is whether the symmetrized bounds are independent of the choice of the optimal solution used to symmetrize the constraints. Unfortunately, the answer is no as shown by the following example. Let (P) be given by min{x2 − x1 : x1 − x2 = 0, x2 + x3 = 1, x ≥ 0}. Then (P) has multiple optimal solutions given by (x1, x2, x3) = (β, β, 1 − β) where β ∈ [0, 1], with an optimal value of 0. If the right-hand side is perturbed to (0, 1)T + t (2, 1)T , then the reader can easily verify that (AUX1) yields (−1/3, +∞) as the optimal partition bounds, whereas the symmetrized bounds are (−β, +β) if one uses the optimal solutions with β < 1/3 to symmetrize the constraints, and (−1/3, 1/3) if those with β ≥ 1/3 are used. This example illustrates that in case of multiple optimal solutions, the symmetrized bounds are dependent on the optimal solution used in the formulation of the symmetrized LPs. Therefore, the relationship (2.7) no longer holds between the optimal partition bounds and the symmetrized bounds. SENSITIVITY IN DEGENERATE LINEAR PROGRAMS 7 However, we will keep using the symmetrized LPs for two reasons. First of all, at least in the unique solution case, they bear a nice relationship to the auxiliary LPs. For our analysis, we will always choose an optimal solution in the relative interior of the optimal set; therefore the symmetrization will hopefully allow more room for the decision variables of the symmetrized LPs. Secondly, the symmetrized LPs are easier to deal with than the auxiliary LPs and the symmetrized bounds will provide a good comparison basis for our interior-point approach proposed in [13], as will be analyzed in the subsequent sections. 2.2. Interior-Point Approach and Central Path Neighborhoods. We will start with a brief review of the primal-dual path-following interior-point methods. The reader is referred to [11] for an extensive treatment. The central path is a path of strictly feasible points (x(ν), y(ν), s(ν)) parametrized by a positive scalar ν. Each point on the central path satisfies the following system for some ν > 0: (2.10) Ax = b, AT y + s = c, XSe = νe, with x > 0 and s > 0. Under Assumptions 2.1 and 2.2, such a solution exists and is unique for each positive ν. Interior-point methods are iterative algorithms that generate iterates which “follow” the central path in the direction of decreasing ν towards the primal-dual optimal set ΩP × ΩD. The iterates generated typically lie in some neighborhood of the central path. For any given feasible iterate (x, y, s), the duality gap is given by cT x − bT y = xT s ≥ 0 and we define the duality measure µ as µ := µ(x, s) := xT s/n. Let S and S0 denote the set of feasible and strictly feasible primal-dual points respectively, that is, (2.11) (2.12) S = {(x, y, s) : Ax = b, AT y + s = c, (x, s) ≥ 0}, S0 = {(x, y, s) ∈ S : (x, s) > 0}. One of the commonly used neighborhoods in interior-point methods is the so-called wide neighborhood, denoted by N−∞(γ): (2.13) N−∞(γ) = {(x, y, s) ∈ S0 : xisi ≥ γµ, ∀ i = 1, 2, . . . , n}, where γ ∈ (0, 1]. At each iteration, given (x, y, s) ∈ N−∞(γ), the algorithm determines a search direction (∆x, ∆y, ∆s). This direction is usually obtained by seeking an approximation to the point on the central path corresponding to some parameter ν ≤ µ, and then applying Newton’s method to the nonlinear system of equations (2.10). Finally, a (damped) step is taken in this direction in such a way that the resulting iterate still lies in N−∞(γ). However, as in the context of target-following methods [10, Part III], one might seek an approximation to a point other than the one on the central path. It suffices to redefine (2.10) by replacing νe in the right-hand side of the third equality by any target vector v > 0. In this case, the Newton step at (x, y, s) for the target vector v is given by: (2.14) A∆x = b − Ax, AT ∆y + ∆s = c − AT y, S∆x + X∆s = v − XSe. 8 E. ALPER YILDIRIM AND MICHAEL J. TODD Next, we describe the interior-point approach proposed by the authors in [13]. Given a primal-dual pair of LPs (P) and (D), let us assume that b or c is perturbed in some fixed direction. Assuming (x, y, s) is strictly primal-dual feasible for (P) and (D), a full Newton step is taken from (x, y, s) for the target vector v := XSe for the perturbed LP pair, thereby aiming to maintain the current xisi products. If (x, y, s) is near-optimal for (P) and (D), then this particular choice is likely to result in a near-optimal solution for the perturbed LP pair since XSe ≈ 0. We state the results formally, referring the reader to [13] for the proofs. Note, in particular, that the duality gap of the resulting feasible iterate for the perturbed LP pair is no greater than that of the original iterate. Proposition 2.2. Assume that (x, y, s) is a strictly feasible point for (P) and (D) and the right-hand side vector b is replaced by b + t∆b, where t ∈ R and ∆b ∈ Rm. Suppose a Newton step is taken from (x, y, s) for the target vector v := XSe for the perturbed problem. Then a full Newton step will yield a feasible iterate for the new problem if and only if (2.15) | t |≤ 1 S−1AT (AD2AT )−1∆b , ∞ where D = X 1 2 S− 1 2 . Moreover, in this case the new iterate will have duality gap at most xT s. Proposition 2.3. Assume that (x, y, s) is a strictly feasible point for (P) and (D) and the cost vector c is replaced by c + t∆c, where t ∈ R and ∆c ∈ Rn. Suppose a Newton step is taken from (x, y, s) for the target vector v := XSe for the perturbed problem. Then a full Newton step will yield a feasible iterate for the new problem if and only if (2.16) | t |≤ S−1(I 1 − AT (AD2AT )−1AD2)∆c , ∞ where D = X 1 2 S− 1 2 . Moreover, in this case the new iterate will have duality gap at most xT s. Under primal-dual nondegeneracy, the bounds arising from Propositions 2.2 and 2.3 computed at near-optimal solutions for (P) and (D) asymptotically equal the symmetrized bounds arising from (SA1) and (SA2) [13]. The goal of this paper is to investigate the quality of these bounds in the absence of the nondegeneracy assumption. We first present a nice characterization of the distance of the strictly feasible primal-dual points (x, y, s) from strictly complementary optimal solutions in terms of the duality gap µn. Using this characterization, we derive some bounds on the components of such points. In what follows, xB, xN , sB and sN denote the partitions of x and s according to the optimal partition (B, N ) as before. Furthermore, we will use the bounds O(µ), Ω(µ) and Θ(µ) interchangeably for scalars as well as vectors and matrices by which we mean each entry satisfies the stated bounds. O(µ) will indicate that the quantity (in absolute value) is bounded above by some positive multiple of µ, where the multiple depends on the primal-dual instance (P) and (D) but does not depend on the particular strictly feasible point or on µ. Similarly, Ω(µ) will indicate a lower bound by some positive multiple of µ and Θ(µ) will mean a lower and upper bound by some positive multiples of µ. The following proposition will be useful for the analysis that follows. Actually, the proposition continues to hold for any feasible solutions and even for a point where SENSITIVITY IN DEGENERATE LINEAR PROGRAMS 9 feasibility is violated by O(µ). The statement below suffices for the purposes of this paper. Proposition 2.4. Let (x, y, s) be a strictly feasible point for (P) and (D) with du- ality gap µn. Then, there exists a strictly complementary optimal solution (x∗, y∗, s∗) of (P) and (D) such that (2.17) (x, y, s) = (x∗, y∗, s∗) + O(µ). Proof. Optimal solutions of (P) and (D) satisfy the linear system Ax = b, AT y + s = c, cT x − bT y = 0, x ≥ 0, s ≥ 0. Any strictly feasible point (x, y, s) satisfies the same linear system with the third equality replaced by cT x − bT y = µn. Hoffman’s lemma [6] indicates that there exists a solution (xˆ, yˆ, sˆ) of the first sys- tem such that (xˆ, yˆ, sˆ) = (x, y, s) + O(µ). The result follows immediately if (xˆ, yˆ, sˆ) is strictly complementary. If not, there exists an arbitrarily small perturbation of (xˆ, yˆ, sˆ) which leads to a strictly complementary solution and (2.17) follows since µ > 0. The following corollary immediately follows from Proposition 2.4 since x∗N = 0 and s∗B = 0 for any optimal solution of (P) and (D). Corollary 2.5. Let (x, y, s) be a strictly feasible point for (P) and (D) with duality gap µn. Then, (2.18) xN = O(µ), sB = O(µ). Note that both Proposition 2.4 and Corollary 2.5 hold for any primal-dual strictly feasible (x, y, s). Next, we derive some more bounds by restricting the iterates to lie in a wide neighborhood given by (2.13). Proposition 2.6. Let (x, y, s) ∈ N−∞(γ) with duality gap µn for (P) and (D). Then, (2.19)XSe = Θ(µ), sN = Ω(1), xB = Ω(1), XN SN−1e = O(µ), SBXB−1e = O(µ). Proof. Since (x, y, s) ∈ N−∞(γ), we have xisi ≥ γµ. Moreover, xT s = µn. Therefore, xisi ≤ µn since x > 0 and s > 0. This proves XSe = Θ(µ). By Corollary 2.5, xN = O(µ). Then, XSe = Θ(µ) implies sN = Ω(1). A similar argument shows xB = Ω(1). Finally, xN = O(µ) together The proof of SBXB−1e = O(µ) is similar. with sN = Ω(1) imply XN SN−1e = O(µ). We mentioned in Section 2 that Assumption 2.2 is somewhat restrictive. Actually, it is possible to extend our interior-point approach to LP instances where that particular assumption is not satisfied. We simply need to define an appropriate viewpoint. Assume that (x, y, s) is such that Ax = b + ξb, AT y + s = c + ξc, x > 0 and s > 0 with ξb = O(µ) and ξc = O(µ), where µn is the duality gap. (Such a pair of solutions will be generated by several infeasible-interior-point methods when applied to problems where optimal solutions exist.) It follows that Proposition 2.4 and Corollary 2.5 hold for such a point as well as Proposition 2.6 with an appropriate definition of the wide neighborhood. One can then take precisely the same Newton step as before for a perturbed LP pair and obtain an iterate with precisely the same primal and dual infeasibilities with a lower duality gap. Then, the resulting interior-point bound can be interpreted as the range of perturbations for which a single Newton step yields 10 E. ALPER YILDIRIM AND MICHAEL J. TODD a point with a smaller duality gap for a nearby LP instance of the perturbed problem. The analysis of Sections 4 and 5 carry over in this case and will reveal that the interior-point bounds will still be related to the desired symmetrized partition bounds as µ tends to 0. Another possible extension of our interior-point approach in this case is to aim to correct for the entire primal and dual feasibilities in a single Newton step by setting the right-hand side of the first and second equations in (2.14) to t∆b − ξb and −ξc, respectively. The analysis remains unchanged for the case of unique and nondegenerate primal-dual optimal solution (since B will be nonsingular). However, ξb and ξc might not be in the “right” spaces in the degenerate case, which would complicate the analysis. Consequently, we adopt the viewpoint described in the previous paragraph for LP instances failing to satisfy Assumption 2.2. 3. Equivalence. In this section, we show that the interior-point bounds are independent of the problem formulation. It is well-known that although (P) and (D) do not look symmetric, they can easily be reformulated so that (D) is in the form of (P) and vice versa [10, pp. 110–112]. More precisely, (D) is equivalent to min xˆT s, subject to Ks = cˆ, s ≥ 0, (D’) s where xˆ satisfies Axˆ = b, cˆ := Kc, and K ∈ R(n−m)×n is a matrix whose rows form a basis for the null space of A. The dual of (D’) is given by max cˆT u, subject to KT u + x = xˆ, x ≥ 0, (P’) u,x which is equivalent to (P). Let us now focus on perturbations of c, i.e., let c be replaced by c + t∆c. By the above reformulation, this is the same as replacing the right-hand side of (D’) by cˆ + tK∆c. Therefore, Proposition 2.2 can be used to evaluate the interior-point bound at a strictly feasible primal-dual pair (s, x) (note that the roles of x and s are interchanged). We need to compute (3.1) X−1KT (KSX−1KT )−1K∆c. On the other hand, one can also use Proposition 2.3 to compute the interior-point bound directly at (x, s), which requires the evaluation of (3.2) S−1(I − AT (AXS−1AT )−1AXS−1)∆c. A simple manipulation of (3.1) gives rise to another equivalent formula: (3.3) X −1/2S−1/2ΨX 1/2S−1/2∆c, where Ψ is the orthogonal projection matrix onto the range space of X −1/2S1/2KT . Similarly, (3.2) is equivalent to (3.4) X −1/2S−1/2ΞX 1/2S−1/2∆c, where Ξ is the orthogonal projection matrix onto the null space of AX 1/2S−1/2. Therefore, in order to prove that (3.1) and (3.2) are equivalent, it suffices to show that Ψ and Ξ project onto the same subspace, or that the null space of AX 1/2S−1/2 equals the range space of X−1/2S1/2KT . This is easily proven by an inclusion argument: if SENSITIVITY IN DEGENERATE LINEAR PROGRAMS 11 w satisfies AX1/2S−1/2w = 0, then X1/2S−1/2w = KT u for some unique u. Thus, w is in the range space of X −1/2S1/2KT . Conversely, if w = X−1/2S1/2KT u for some u, then AX1/2S−1/2w = AKT u = 0. This proves the equivalence of the interior-point bounds. We next argue that the range of t resulting from the optimal partition bounds is also independent of the formulation. If the two LPs are formulated in the form of (P) and (D), then (AUX2) yields the range of t for perturbations of c. Premultiplying the equality constraints of (AUX2) by K = [KB, KN ] leads to (AUX1’) given by (3.5) min (max) λ w,λ w,λ s.t. KN w = λK ∆c, w ≥ −s∗N , (AUX1’) which exactly yields the range of t for perturbations of the right-hand side of (D’) if one uses the form (D’) and (P’). Similarly, if (w, λ) is feasible for (AUX1’), then λ∆cB λ∆cN − w lies in the null space of K. Then, by our previous observation, there exists v such that BT v = λ∆cB, N T v + w = λ∆cN , which is exactly the constraints of (AUX2), completing the proof of the claim. Using this observation, we will carry out our analysis for perturbations of b only in the subsequent sections, and state the corresponding results for changes in c as corollaries. We begin with a special case of degeneracy first and then consider the most general case. 4. Unique Primal Solution. Throughout this section, we assume that (P) has a unique but degenerate optimal solution x∗. Note that by Proposition 2.1, we have |B| = rank (B), i.e., B has linearly independent columns. In this particular case, Proposition 2.4 provides another useful bound on xB for a strictly feasible primal-dual point (x, y, s). Corollary 4.1. Assume that (P) has a unique optimal solution x∗. Let (x, y, s) be primal-dual strictly feasible for (P) and (D) with duality gap µn. Then, (4.1) xB = x∗B + O(µ). An analogous corollary follows if (D) has a unique optimal solution. Corollary 4.2. Assume that (D) has a unique optimal solution (y∗, s∗). Let (x, y, s) be primal-dual strictly feasible for (P) and (D) with duality gap µn. Then, (4.2) sN = s∗N + O(µ). Next, we will analyze one-dimensional perturbations of b. 4.1. Perturbations of b. In this subsection, we assume that the right-hand side vector b is replaced by b + t∆b, where ∆b ∈ Rm and t ∈ R. We also assume that (x, y, s) ∈ N−∞(γ) is a primal-dual strictly feasible point for (P) and (D) for some γ ∈ (0, 1]. We will compare the interior-point bounds arising from Proposition 2.2 at (x, y, s) with the optimal values of (SA1), i.e., the symmetrized bounds. The interior-point bounds are given by the L∞-norm of (4.3) S−1AT (AD2AT )−1∆b, 12 E. ALPER YILDIRIM AND MICHAEL J. TODD where D2 = XS−1. Let us now consider (SA1). Since B has full column rank, ∆b either does not lie in the range space of B, in which case the optimal values of (SA1) as well as (AUX1) are all 0, or there exists a unique v ∈ R|B| such that Bv = ∆b, in which case the constraints of (SA1) reduce to (2.9), from which the symmetrized bounds can be obtained easily. We will consider both situations in turn. Let us start with the second case. Without loss of generality, we can assume that ∆b has unit L2-norm, which implies a bound on v. Then, we need to compute (4.4) u = (AD2AT )−1∆b = (AD2AT )−1Bv in order to obtain (4.3). However, (4.4) is equivalent to (4.5) AD2AT u = Bv or BXBSB−1BT u + N XN SN−1N T u = Bv, where B and N are the partitions of the coefficient matrix A with respect to B and N as before. Since B has linearly independent columns, there exists a matrix C ∈ Rm×(m−|B|) such that the augmented matrix [B C] is square and nonsingular: let W be its inverse. Therefore, premultiplying the second equality in (4.5) by W , we obtain (4.6) I 0 XBSB−1 [I 0] u˜ + N˜ XN SN−1N˜ T u˜ = I 0 v, where u˜ = W −T u, N˜ = W N and I is a |B| × |B| identity matrix. Therefore, if we partition N˜ and u˜ accordingly as N˜ = N˜1 N˜2 , u˜ = u˜1 u˜2 , (4.6) can then be decomposed in the following way: (4.7) DB2 + N˜1DN2 N˜1T N˜1DN2 N˜2T N˜2DN2 N˜1T N˜2DN2 N˜2T u˜1 u˜2 = v 0 , where DB and DN are the corresponding partitions of D. By (4.3), we need to compute (4.8) S−1AT u = S−1(W A)T u˜ = SB−1u˜1 SN−1 N˜1T u˜1 + N˜2T u˜2 . For notational convenience, let us define F := N˜1DN , G := N˜2DN . Note that G has full row rank since A does. The bottom equality in (4.7) can be rewritten as (4.9) GF T u˜1 + GGT u˜2 = 0, so u˜2 = −(GGT )−1GF T u˜1. Substituting (4.9) in the top equality in (4.7) gives (4.10) DB2 + F F T − F GT (GGT )−1GF T u˜1 = v, DB2 + F (I − PG)F T u˜1 = v, or SENSITIVITY IN DEGENERATE LINEAR PROGRAMS 13 where PG is the orthogonal projection matrix onto the range space of GT . Therefore, I − PG is the orthogonal projection matrix onto the null space of G. We briefly review the Neumann lemma now [3]. Let U be an invertible matrix and let V satisfy U −1V ≤ 1/2. (The particular norm being used does not really matter: we will always use · for the Euclidean norm or the operator norm arising from it.) Then, I + U −1V is invertible with I + U −1V ≤ 2. Moreover U + V is invertible and given by (4.11) (U + V )−1 = U −1 − U −1V (I + U −1V )−1U −1. Now, we apply this result to (4.10) with U := DB2 and V := F (I −PG)F T . Proposition 2.6 implies that both U −1 and V are O(µ) since I − PG is a projection matrix and has unit Euclidean norm. Therefore, assuming the duality gap µn is small, (4.12) u˜1 = DB2 + F (I − PG)F T −1 v, = DB−2v − DB−2F (I − PG)F T I + DB−2F (I − PG)F T −1 DB−2v. It then follows that (4.13) SB−1u˜1 = XB−1v − XB−1F (I − PG)F T I + DB−2F (I − PG)F T −1 DB−2v. However, by Proposition 2.6, F is O(µ1/2), DB−2 is O(µ) and XB−1 is O(1). Consequently, the second term on the right hand side of (4.13) is O(µ2) since I − PG ≤ 1. Finally, Corollary 4.1 implies XB−1 = (XB∗ )−1 + O(µ). Therefore, (4.14) SB−1u˜1 = (XB∗ )−1v + O(µ). We have thus obtained the top part of (4.8). For the lower part, we get (4.15) SN−1(N˜1T u˜1 + N˜2T u˜2) = SN−1DN−1(F T u˜1 + GT u˜2), = (XN SN )−1/2(I − PG)F T u˜1, where we substituted (4.9) for u˜2. Proposition 2.6 implies (XN SN )−1/2 = O(µ−1/2) and F = O(µ1/2). By (4.14), u˜1 = O(µ) since sB is. Combining these bounds with I − PG ≤ 1 leads to (4.16) SN−1(N˜1T u˜1 + N˜2T u˜2) = O(µ). Using (4.8), we conclude that the L∞-norm of the quantity (4.3) we need to evaluate is given by (4.17) S−1AT (AD2AT )−1∆b ∞ = (XB∗ )−1v + O(µ) O(µ) . ∞ The reciprocal of (4.17) gives the desired interior-point bound. Consequently, if the duality gap µn is small, we conclude by comparing (4.17) with (2.9) that the interiorpoint approach yields exactly the same bound as the optimal solution to (SA1) asymptotically in µ. Next, we address the situation where ∆b does not lie in the range space of B. In this case, the optimal values of both (AUX1) and (SA1) are clearly 0. ∆b can be uniquely written as (4.18) ∆b = BvB + CvC , 14 E. ALPER YILDIRIM AND MICHAEL J. TODD where [B C] is nonsingular as before and vC is a nonzero vector. Once again, we need to compute (4.3). We follow a similar treatment as before, and corresponding to (4.7) we obtain: (4.19) DB2 + N˜1DN2 N˜1T N˜1DN2 N˜2T N˜2DN2 N˜1T N˜2DN2 N˜2T The bottom part can be expanded as u˜1 u˜2 = vB vC . (4.20) N˜2XN SN−1N˜1T u˜1 + SN−1N˜2T u˜2 = vC . However, (4.8) implies that the term in the brackets is exactly the bottom part of the quantity (4.3) we seek. Let us denote that term by p and let XN p = q. Then, (4.20) is equivalent to N˜2 q = vC. Since vC is nonzero, the norm of q is bounded below, that is, q ≥ α > 0 where α is the norm of the least squares solution. Therefore, q ∞ ≥ β with β := (α/ n − |B|) (see e.g. [3]). (Note that |B| < n since |B| = n can happen only if m = n, in which case ∆b is always in the range of B.) However, q ∞ ≤ XN ∞ p ∞ since XN p = q. This implies (4.21) p ∞≥ q∞ ≥ XN ∞ β XN = Ω(1/µ), ∞ where the last equality follows from Corollary 2.5. Therefore, as µ tends to 0, p ∞ tends to ∞, which implies that the interior-point bound given by its reciprocal tends to 0 as desired. We remark that if B = ∅, then x∗ = 0 is the only optimal solution of (P), which can happen only if b = 0. In this case, the top part of (4.8) disappears. The interior-point bound is then given by the reciprocal of p ∞, where p is as defined after (4.20). By the preceding argument, the interior-point bound tends to 0 as µ approaches 0. This is still in agreement with the optimal partition bounds since any nonzero perturbation of b leads to a change in the optimal partition and hence, the optimal partition bounds in this case are also equal to 0. Therefore, we have proved the following theorem: Theorem 4.3. Let (x, y, s) ∈ N−∞(γ) be a primal-dual strictly feasible point for (P) and (D). Assume that (P) has a unique but degenerate optimal solution and that b is replaced by b + t∆b where t ∈ R and ∆b ∈ Rm. Then the interior-point bound evaluated at (x, y, s) yields exactly the same value as the optimal solution of (SA1) asymptotically in µ, where µ = xT s/n. The following corollary of Theorem 4.3 is an immediate consequence of the equivalence between (P) and (D) as discussed in Section 3. One uses Corollary 4.2 in place of Corollary 4.1 in the preceding analysis. Corollary 4.4. Let (x, y, s) ∈ N−∞(γ) be a primal-dual strictly feasible point for (P) and (D). Assume that (D) has a unique but degenerate optimal solution and that c is replaced by c + t∆c where t ∈ R and ∆c ∈ Rn. Then the interior-point bound evaluated at (x, y, s) yields exactly the same value as the optimal solution of (SA2) asymptotically in µ, where µ = xT s/n. It does not appear that we can obtain better results for perturbations of c in the case of a unique primal optimal solution (but not dual optimal solution) than those arising from the analysis of the general case in the next section. A similar remark holds for perturbations of b in the case of a unique dual optimal solution (but not primal optimal solution). SENSITIVITY IN DEGENERATE LINEAR PROGRAMS 15 5. General Case. In this section, we turn our attention to the most general case where both (P) and (D) may have multiple optimal solutions. As the small example given at the end of Section 2.1 reveals, some complications arise in the presence of multiple optimal solutions. For instance, unlike the previous case, the symmetrized bounds become dependent on the optimal solution of (P) used in the formulation of (SA1) if (P) has multiple optimal solutions. Furthermore, they do not necessarily coincide with the “symmetrizations” of the optimal partition bounds arising from (AUX1). Similar remarks hold for the relationship between (SA2) and (AUX2) if (D) has multiple optimal solutions. Despite this complication arising from the presence of multiple optimal solutions, we aim to be able to say something about the quality of the interior-point bounds at least in comparison with the symmetrized bounds. In the next subsection, we analyze perturbations of b in this general setting. 5.1. Perturbations of b. Let (P) have multiple optimal solutions and let b be replaced by b + t∆b, where t ∈ R and ∆b ∈ Rm. Suppose that (x, y, s) ∈ N−∞(γ) is primal-dual strictly feasible where γ ∈ (0, 1]. For such a point, Proposition 2.4 guarantees the existence of a strictly complementary solution (x∗, y∗, s∗) whose dis- tance from (x, y, s) is O(µ). We will compare the interior-point bounds evaluated at (x, y, s) with the optimal values of (SA1). Among other optimal solutions of (P), the x∗ above will be the particular choice of the primal optimal solution to be used in the formulation of (SA1). The use of such an optimal solution in the relative interior of the primal optimal set is likely to leave more room for the decision variables of (SA1) since x∗B > 0. Let us first consider (SA1). The constraints of (SA1) are (5.1) Bu = λ∆b, −x∗B ≤ u ≤ x∗B . Let rank (B) = r and |B| = k. Clearly we have r ≤ m and r < k since Proposition 2.1 implies dim (ΩP ) = k − r, which is positive by our assumption. This, in turn, implies that r > 0 since r = 0 if and only if B = ∅ (assuming no columns of A are identically zero). A QR factorization of B yields B = QR, where Q ∈ Rm×m is orthogonal and R ∈ Rm×k is upper triangular with (5.2) R= R1 0 , where R1 has r rows. Note that R1 has full row rank. Premultiplying the equality constraints in (5.1) by QT yields (5.3) R1 0 u=λ ∆b1 ∆b2 , − x∗B ≤ u ≤ x∗B , with ∆b = QT ∆b. Clearly, (5.3) reveals that (SA1) has a nontrivial optimal solution λ∗ if and only if ∆b2 = 0. First, we consider the nontrivial case. (Since ∆b is nonzero, this implies that k > 0.) Let (λ∗, u∗) be an optimal solution to the maximization problem with λ∗ = 0. Note that λ∗ is finite since u∗ is bounded (this follows since B = ∅). Then, we have (5.4) ∆b = QT ∆b = (1/λ∗)Ru∗. 16 E. ALPER YILDIRIM AND MICHAEL J. TODD The interior-point approach, on the other hand, requires the evaluation of (4.3) at (x, y, s). By (5.4), we then need to evaluate the L∞-norm of (5.5) (1/λ∗)S−1AT (AD2AT )−1QRu∗. Let (5.6) w = (AD2AT )−1QRu∗ or AD2AT w = QRu∗. Premultiplying the second equality in (5.6) by QT gives (5.7) R1 N˜1 0 N˜2 DB2 0 0 DN2 R1T 0 N˜1T N˜2T w˜1 w˜2 = R1 0 u∗, where w˜1 and w˜2 are the appropriate partitions of w˜ = QT w and N˜1 and N˜2 are those of N˜ = QT N . Let us define (5.8) F := N˜1DN , G := R1DB, H := N˜2DN . (5.7) can then be decomposed into two equations as (5.9) (GGT + F F T )w˜1 + F HT w˜2 = R1u∗, HF T w˜1 + HHT w˜2 = 0. Note, in particular, that both G and H have full row rank since R1 and A do. From the second equation in (5.9), we obtain (5.10) w˜2 = −(HHT )−1HF T w˜1. Substituting (5.10) in the first equation of (5.9) leads to (5.11) (GGT + F F T − F HT (HHT )−1HF T )w˜1 = R1u∗, or (GGT + F (I − PH )F T )w˜1 = R1u∗, where I −PH is the orthogonal projection matrix onto the null space of H. Proposition 2.6 implies that the second term in parentheses in the second equation above is O(µ) since I − PH ≤ 1. In order to apply Neumann’s lemma, we need to show that (GGT )−1 is bounded. Lemma 5.1. (GGT )−1 = O(µ). Proof. We use the “thin” QR factorization of GT = DBR1T = Y Z, where Y has orthonormal columns and Z is upper triangular and nonsingular. Then, (GGT )−1 = Z−1Z−T . Therefore, it suffices to find an upper bound on Z−1. We have (5.12) DBR1T = Y Z, or R1R1T = R1DB−1Y Z. Therefore, I = Proposition 2.6, (DRB−1R1 1T=)−O1(Rµ11D/2B−),1YwhZi,chorimZp−li1es=th(aRt 1ZR−1T1)−=1 R1DB−1Y O(µ1/2 ) . However, completing by the proof. We can now apply Neumann’s lemma to (5.11). Using the same notation as in (4.11) we have U := GGT and V := F (I − PH )F T . Note that both U −1 and V are O(µ). We obtain (5.13) w˜1 = (GGT )−1GDB−1u∗ − (GGT )−1V (I + (GGT )−1V )−1(GGT )−1GDB−1u∗, SENSITIVITY IN DEGENERATE LINEAR PROGRAMS 17 where we used R1 = GDB−1. By (5.5) and (5.6), we need (5.14) (1/λ∗)S−1AT w = (1/λ∗) Let us define SB−1R1T w˜1 SN−1(N˜1T w˜1 + N˜2T w˜2) . (5.15) W = GT (GGT )−1. For the top part of (5.14) we need to evaluate (5.16) SB−1R1T w˜1 = SB−1DB−1GT w˜1, = (SB XB )−1/2PGDB−1u∗ − (SBXB)−1/2W V (I + (GGT )−1V )−1W T DB−1u∗, where we used (5.13), (5.15) and where PG is the orthogonal projection matrix onto the range space of GT . Consider the second term in the right hand side of the second equality. By Proposition 2.6, (SBXB)−1/2 is O(µ−1/2), V is O(µ) and DB−1 is O(µ1/2). Lemma 5.1 implies that W = GT (GGT )−1 = Y Z−T = O(µ1/2) since Y ≤ 1. Therefore, the whole expression is O(µ2). We conclude that the top part of (5.14) is (5.17) (1/λ∗)(SBXB)−1/2PG(SBXB)1/2XB−1u∗ + (1/λ∗)O(µ2). Let us next consider the lower part of (5.14). We need to compute (5.18) SN−1N˜1T w˜1 + SN−1N˜2T w˜2. By (5.13) the first term in (5.18) is given by (5.19) (SN XN )−1/2F T W T − (GGT )−1V (I + (GGT )−1V )−1W T DB−1u∗. Note that W = O(µ1/2) by the preceding discussion. As for the second term in brackets, we have both (GGT )−1 and V are O(µ), which implies the whole second term is O(µ5/2). Thus, the expression in brackets is O(µ1/2). By Proposition 2.6, (SN XN )−1/2 is O(µ−1/2), whereas both F T and DB−1 are O(µ1/2). We therefore conclude that (5.19) is O(µ). For the second term in (5.18), we use (5.10) together with (5.13): (5.20) −(SN XN )−1/2HT (HHT )−1HF T w˜1 = −(SN XN )−1/2PH F T w˜1. Note that w˜1 = O(µ) by the preceding arguments. The fact that PH ≤ 1 together with (SN XN )−1/2 being O(µ−1/2) and F T being O(µ1/2) implies (5.20) is O(µ). Therefore, we conclude that the lower part of (5.14) is O(µ). Combining this result with (5.17) yields the following: (5.21) r := (1/λ∗) (SB XB )−1/2PG(SB XB )1/2XB−1u∗ + O(µ2) O(µ) . Consequently, we need to evaluate the L∞-norm of (5.21) and take its reciprocal. Observe that XB−1 = (XB∗ )−1 + O(µ) by Proposition 2.4. Using this, we derive an upper bound on the L∞-norm of (5.21). r ∞ ≤ |1/λ∗| (SBXB)−1/2 ∞ PG ∞ (SBXB)1/2 ∞ (XB∗ )−1u∗ ∞ + O(µ) . (5.22) 18 E. ALPER YILDIRIM AND MICHAEL J. TODD Since (x, y, s) ∈ N−∞(γ), (5.23) xisi = µn − xj sj ≤ µn − µ(n − 1)γ = µ(n − (n − 1)γ). j=i Thus, (xisi)1/2 ≤ (µ(n − (n − 1)γ)1/2 and (xisi)−1/2 ≤ 1/(γµ)1/2. Furthermore, since PG ≤ 1, we have PG ∞ ≤ k1/2 (see e.g. [3]), where |B| = k. Finally, since u∗ is optimal for (SA1), (XB∗ )−1u∗ ∞ ≤ 1. Therefore, r ∞≤ 1 |λ∗| 1 (γµ)1/2 k1/2(µ(n − (n − 1)γ))1/2 + O(µ) , (5.24) = 1 |λ∗| (n − (n − γ 1)γ)k + O(µ) . We conclude that the interior-point bound, which is the reciprocal of (5.24), is then bounded below by (5.25) 1 r∞ ≥ √γ |λ∗|. (n − (n − 1)γ)k + O(µ) √ Note, in particular, that the lower bound tends to 1/ k, independent of n, as µ → 0 if (x, y, s) is on the central path. We next consider the case where ∆b is not in the range space of B. Again, in this case, the symmetrized bounds as well as the optimal partition bounds are all 0. The QR factorization of B can be rewritten as B = QR = [Q1 Q2]R = Q1R1, where we use (5.2) and [Q1 Q2] is the appropriate partition of Q. Since Q is orthogonal, ∆b can uniquely be expressed as (5.26) ∆b = Q1v1 + Q2v2, where v2 is nonzero. Arguing similarly to Section 4, we need to evaluate (4.3), which in turn requires the computation of (5.27) w = (AD2AT )−1(Q1v1 + Q2v2) or AD2AT w = Q1v1 + Q2v2. Premultiplying (5.27) by QT leads to (5.28) R1 N˜1 0 N˜2 DB2 0 0 DN2 R1T 0 N˜1T N˜2T w˜1 w˜2 = v1 v2 , which looks like (4.19). Essentially the same arguments as in Section 4 reveal that the interior-point bound tends to 0 as µ approaches 0. Therefore, we have proved the following theorem. Theorem 5.2. Let (x, y, s) ∈ N−∞(γ) be a primal-dual strictly feasible point for (P) and (D) with duality gap µn. Assume that (P) has multiple optimal solutions and that b is replaced by b + t∆b where t ∈ R and ∆b ∈ Rm. Let k = |B|. If the strictly feasible solution given by Proposition 2.4 is used for symmetrization in (SA1), then the ratio of the interior-point bound evaluated at (x, y, s) to the value of the optimal solution of (SA1) is at least (5.29) √γ . (n − (n − 1)γ)k + O(µ) SENSITIVITY IN DEGENERATE LINEAR PROGRAMS 19 Note that the presence of multiple primal optimal solutions implies k > 0, there- fore, the expression (5.29) is well-defined. As in Section 4, Theorem 5.2 leads to the following corollary by the discussion in Section 3. Due to the interchange of the roles of the basic and nonbasic variables in the reformulation given in Section 3, k in the denominator of (5.29) is replaced by (n − k). Under the assumption of multiple dual optimal solutions, Proposition 2.1 indicates that m > r, which implies k < n since A has full row rank. Corollary 5.3. Let (x, y, s) ∈ N−∞(γ) be a primal-dual strictly feasible point for (P) and (D) with duality gap µn. Assume that (D) has multiple optimal solutions and that c is replaced by c + t∆c where t ∈ R and ∆c ∈ Rn. Let k = |B|. If the strictly feasible solution given by Proposition 2.4 is used for symmetrization in (SA2), then the ratio of the interior-point bound evaluated at (x, y, s) to the value of the optimal solution of (SA2) is at least (5.30) √γ . (n − (n − 1)γ)(n − k) + O(µ) 6. An Example. In the previous sections, we have provided a theoretical basis for the behavior of the interior-point bounds evaluated at the near-optimal solutions. We present an example in this section to shed some light on the performance of the interior-point bounds in practice. The example we consider in this section is a degenerate transportation problem discussed by Roos, Terlaky and Vial [10, pp. 398–402]. For this problem, they report strikingly different results on sensitivity analysis on the right-hand side and the cost parameters using different commercial solvers. We aim to test our interior-point approach on this instance. The problem is very simple. There are three distribution centers with capacity 2, 6, and 5 units, respectively, which can supply three warehouses each with a demand of 3 units at a cost of 1 per unit. We aim to minimize the total cost while meeting all the demand. This problem can be formulated as a linear program in standard form as follows: min 3 i=1 3 j=1 xij subject to x11 + x12 + x13 + s1 = 2, x21 + x22 + x23 + s2 = 6, x31 + x32 + x33 + s3 = 5, x11 + x21 + x31 − d1 = 3, x12 + x22 + x32 − d2 = 3, x13 + x23 + x33 − d3 = 3, xij , si, dj ≥ 0, i, j = 1, 2, 3, where xij denotes the amount shipped from distribution center i to warehouse j, si is the excess supply at distribution center i and dj is the shortage of demand at warehouse j, i, j = 1, 2, 3. It is easy to verify that any feasible solution with objective value 9 is optimal. Therefore, there exist optimal solutions with all the variables xij and si, i, j = 1, 2, 3, positive, whereas all dj, j = 1, 2, 3, will always be zero in any optimal solution. By Proposition 2.1, the primal optimal set has dimension 6 whereas the dual optimal solution is unique. 20 E. ALPER YILDIRIM AND MICHAEL J. TODD We solved this LP using CPLEX’s standard barrier solver (http://www.ilog.com /products/cplex/product/barrier.cfm). By setting the tolerance level to different values (1e-3, 1e-4, 1e-5) we obtained several strictly feasible, near-optimal solutions with different duality gaps. We then evaluated the interior-point bounds (2.15) and (2.16) at these iterates for perturbations of the form b + tek and c + tel, where ek and el with k = 1, . . . , 6 and l = 1, . . . , 9 denote the unit vectors in the appropriate dimension. The results are reported in Table 6.1 and Table 6.2 for perturbations of the righthand side parameters and the cost parameters, respectively. Rows 1–3 correspond to the three strictly feasible, near-optimal iterates ordered by descending duality gaps. For the iterates (x, y, s), µ denotes the duality measure given by xT s/n, and γ is the parameter of the narrowest wide central-path neighborhood containing the iterate. The columns bi in Table 6.1 correspond to changes in the right-hand side of the ith constraint and the number in each column is the upper interior-point bound evaluated at the corresponding iterate. Since the changes in b4–b6 are symmetric, we report only the results for b4. Similarly, the columns cij in Table 6.2 represent changes in the cost coefficient of the variable xij and the number in that column is the upper interiorpoint bound evaluated at the corresponding iterate. Again, the changes in cij for fixed i are symmetric. The CPLEX row denotes the range obtained from that solver. (The basic variables returned by CPLEX were x11, x21, x22, x23, x33, and s3.) Finally, the last row in each table gives the correct ranges (optimal partition bounds) for the corresponding right-hand side or the cost parameters. Table 6.1 Interior-point bounds for the transportation example (b) 1 2 3 CPLEX Range µ 3.94 e-4 2.62 e-5 1.86 e-6 γ 0.98 0.98 0.74 b1 1.653 1.639 1.487 [−2, 1] (−2, +∞) b2 2.554 2.572 2.641 [−2, 1] (−4, +∞) b3 2.542 2.519 2.332 [−4, +∞) (−4, +∞) b4 2.496 2.480 2.163 -1 (−3, 4) Table 6.1 reveals that the interior-point bounds provide useful information in comparison with the symmetrized optimal partition bounds as µ tends to 0 even if the LP is highly degenerate. Furthermore, the ratio of the interior-point bounds to the symmetrized optimal partition bounds is much better than the theoretical worst-case ratio (5.29). In fact, we experienced similar behavior in our extensive computational tests with randomly generated problems [12]. The results of Table 6.2 indicate the rapid convergence of the interior-point bounds to 0 as µ tends to 0 as expected. Note that the convergence rate is related to the duality measure µ as illustrated by the previous theoretical results. The particular instance we considered has multiple primal optimal solutions and a unique degenerate dual optimal solution. In order to get a complementary view, we slightly modified the data of the transportation problem so that the primal problem has a unique degenerate optimal solution whereas the dual problem has multiple optimal solutions. More specifically, we increased the cost coefficients of x11, x12, x13, x21, x32 and x33 from 1 to 2 in the objective function and maintained the remaining 1CPLEX provided the following different results: b4 : [−1, 2], b5 : [−1, 2], b6 : [−1, 4]. SENSITIVITY IN DEGENERATE LINEAR PROGRAMS Table 6.2 Interior-point bounds for the transportation example (c) 1 2 3 CPLEX Range µ 3.94 e-4 2.62 e-5 1.86 e-6 γ 0.98 0.98 0.74 c1j 1.126 e-3 7.505 e-5 5.425 e-6 -2 {0} c2j 7.825 e-4 5.238 e-5 3.740 e-6 -2 {0} c3j 6.994 e-4 4.685 e-5 3.437 e-6 -2 {0} 21 data. The resulting primal instance has the unique optimal solution given by x22 = 3, x23 = 3, x31 = 3, s1 = 2 and s3 = 2 with all other variables equal to 0. We tested our interior-point approach on this problem instance. The results are tabulated in Tables 6.3 and 6.4 for the right-hand side and the cost parameters, respectively. (The basic variables returned by CPLEX in this example were x22, x23, x31, s1, s2, and s3.) For perturbations of the right-hand side, Table 6.3 illustrates the convergence behavior predicted by the theoretical results. For perturbations of the cost vector, Table 6.4 reveals that the interior-point bounds provide very useful information in comparison with the symmetrized partition bounds at a moderate cost. 7. Conclusion. In this paper, we have studied the quality of the bounds arising from the interior-point perspective on sensitivity analysis developed by the authors in [13]. By relaxing the strong assumption of nondegeneracy, we have been able to consider all possible degeneracy scenarios and to investigate how our bounds compare with those arising from the optimal partition approach to sensitivity analysis. If the primal problem has a degenerate but unique optimal solution, then our approach yields the same bounds as the “symmetrized” optimal partition bounds for perturbations of b. By the equivalence discussed in Section 3, the same relationship holds for perturbations of c if the dual problem has a degenerate but unique optimal solution. This result directly extends the previous result proved in [13] under the assumption of a unique and nondegenerate solution. We then considered general degenerate LPs. In this case, we were able to show that our interior-point approach would yield bounds that are at least a certain fraction of the symmetrized bounds, where the fraction depends on certain characteristics of the problem instance and of the iterate at which the bounds are evaluated. The behavior of the interior-point bounds on a highly degenerate transportation example indicates that useful information can be gained using the interior-point approach. Our extensive computational tests on random problems suggest that the ratio in practice is much better than the predicted worst-case ratio. Although this result is not as strong as the aforementioned results, our interior-point bounds still yield some useful information on the range of allowable perturbations. The fact that the cost of the evaluation of our bounds is simply the same as that of an interior-point iteration makes it more appealing given the cost of solving two LPs to obtain the range from the optimal partition approach. 2CPLEX ranges: c11 : (−∞, 0], c12 : [0, +∞), c13 : [0, +∞), c21 : {0}, c22 : [−1, 0], c23 : {0}, c31 : [0, +∞), c32 : [0, +∞), c33 : {0}. Table 6.3 Interior-point bounds for the modified transportation example (b) µγ b1 b2 b3 b4 b5 b6 1 2.42 e-4 0.99 1.999 1.513 e-3 2.000 2.000 1.513 e-3 1.513 e-3 2 1.58 e-5 0.95 2.000 1.012 e-4 2.000 2.000 1.012 e-4 1.012 e-4 3 1.71 e-6 0.58 2.000 8.075 e-6 2.000 2.000 8.075 e-6 8.075 e-6 CPLEX Range [−2, +∞) [0, +∞) [−2, +∞) [−3, 2] [−3, 0] (−2, +∞) {0} (−2, +∞) (−3, 2) {0} [−3, 0] {0} E. ALPER YILDIRIM AND MICHAEL J. TODD Table 6.4 Interior-point bounds for the modified transportation example (c) µγ c11 c12 c13 c21 c22 c23 c31 c32 c33 1 2.42 e-4 0.99 1.000 0.777 0.777 1.335 0.900 0.900 1.000 0.777 0.777 2 1.58 e-5 0.95 1.000 0.777 0.777 1.334 0.900 0.900 1.000 0.777 0.777 3 1.71 e-6 0.58 1.000 0.785 0.785 1.335 0.922 0.922 1.000 0.785 0.785 CPLEX Range [−1, +∞) [−1, +∞) [−1, +∞) [−1, +∞) [−1, 1] [−1, 1] [−1, 1] [−1, +∞) [−1, +∞) (−1, +∞) (−1, +∞) (−1, +∞) (−2, +∞) (−2, 1) (−2, 1) (−1, 1) (−1, +∞) (−1, +∞) 22 SENSITIVITY IN DEGENERATE LINEAR PROGRAMS 23 Acknowledgements. We would like to thank two anonymous referees whose helpful comments and suggestions led to an improved exposition. REFERENCES [1] I. Adler and R. D. C. Monteiro, A geometric view of parametric linear programming, Algorithmica, 8 (1992), pp. 161–176. [2] A. J. Goldman and A. W. Tucker, Theory of linear programming, in Linear Equalities and Related Systems, H. W. Kuhn and A. W. Tucker, eds., Princeton University Press, Princeton, N. J., 1956, pp. 53–97. [3] G. H. Golub and C. F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, third ed., 1996. [4] H. J. Greenberg, Simultaneous primal-dual right-hand-side sensitivity analysis from a strictly complementary solution of a linear program, SIAM Journal on Optimization, 10(2) (2000), pp. 427–442. [5] H. J. Greenberg, A. G. Holder, C. Roos, and T. Terlaky, On the dimension of the set of rim perturbations for optimal partition invariance, SIAM Journal on Optimization, 9 (1999), pp. 207–216. [6] A. J. Hoffman, On approximate solutions of systems of linear inequalities, Journal of Research of the National Bureau of Standards, 49 (1952), pp. 263–265. [7] B. Jansen, J. Jong, C. Roos, and T. Terlaky, Sensitivity analysis in linear programming : Just be careful !, European Journal of Operational Research, 101 (1997), pp. 15–28. [8] W.-J. Kim, C.-K. Park, and S. Park, An -sensitivity analysis in the primal–dual interior point method, European Journal of Operational Research, 116(3) (1999), pp. 629–639. [9] R. D. C. Monteiro and S. Mehrotra, A general parametric analysis approach and its implication to sensitivity analysis in interior point methods, Mathematical Programming, 72 (1996), pp. 65–82. [10] C. Roos, T. Terlaky, and J.-P. Vial, Theory and Algorithms for Linear Optimization : An Interior Point Approach, Wiley-Interscience Series in Discrete Mathematics and Optimization, John Wiley and Sons, 1997. [11] S. J. Wright, Primal-Dual Interior-Point Methods, SIAM Publications, Philadelphia, 1997. [12] E. A. Yıldırım, An interior-point perspective on sensitivity analysis in linear programming and semidefinite programming, PhD thesis, Cornell University, 2001. [13] E. A. Yıldırım and M. J. Todd, Sensitivity analysis in linear programming and semidefi- nite programming using interior-point methods, Mathematical Programming, 90 (2001), pp. 229–261.