IMPROVED COMPLEXITY FOR MAXIMUM VOLUME INSCRIBED ELLIPSOIDS KURT M. ANSTREICHER
Abstract. Let P = fx j Ax bg, where A is an m n matrix. We assume that P contains a ball of radius one centered at the origin, and is contained in a ball of radius R centered at the origin. We consider the problem of approximating the maximum volume ellipsoid inscribed in P . Such ellipsoids have a number of interesting applications, including the inscribed ellipsoid method for convex optimization. We describe an optimization algorithm that obtains an ellipsoid whose volume is at least a factor e? of the maximum possible in O(m3:5 ln(mR=)) operations. Our result provides an alternative to a saddlepoint-based approach with the same complexity, due to Nemirovskii. We also show that a further reduction in complexitycan be obtained by rst computingan approximation of the analytic center of P . Key words. maximum volume inscribed ellipsoid, inscribed ellisoid method. AMS subject classi cations. 90C25, 90C22
1. Introduction. Let P = fx j Ax bg, where A is an m n matrix. We assume that P is bounded, with a nonempty interior. It is then known [6] that there is a unique ellipsoid E P of maximum volume. We say that an ellipsoid E P is
-maximal if Vol(E) Vol(E ), where 0 < < 1 and Vol() denotes n-dimensional volume. In this paper we consider the complexity of computing a -maximal inscribed ellipsoid for P . For convenience in stating complexity results we often write = e? (as in [12, x6.5]), where > 0. There are a number of interesting applications of -maximal ellipsoids. For example, the computation of a -maximal ellipsoid, with > 0:92, is required on each iteration of the inscribed ellipsoid algorithm (IEM) [14] for convex programming. The IEM minimizes a convex function over an n-dimensional cube to relative accuracy in O(n ln(n=)) iterations, each requiring evaluation of the function and a subgradient. The order of this complexity, also achieved by the volumetric cutting plane algorithm [1, 15], is optimal [13]. Another application of -maximal ellipsoids is to provide a \rounding" of P . It is known that for the maximum volume inscribed ellipsoid (MVIE) E , E P nE ; where for an ellipsoid E and positive scalar , E denotes the dilation of E about its center by the factor . For a -maximal ellipsoid E it can be shown [14] that
p
E P n 1 + 3 1 ? E: Roundings of this type are required in several contexts, including Lenstra's algorithm for integer programming in xed dimension [10], and randomized algorithms for volume computation [7]. Alternative methodologies for obtaining O(n)-roundings of P include the shallow cut ellipsoid algorithm [3, x4.6] and the volumetric cutting plane algorithm [2]. Department of Management Sciences, University of Iowa, Iowa City, IA 52242 (
[email protected] 1
2
K. M. ANSTREICHER
Assume that P contains a ball of radius one centered at the origin, and is contained in a ball of radius R centered at the origin. Using the ellipsoid algorithm an e? { maximal inscribed ellipsoid can be computed in O(n6(n2 + m) ln(nR=)) operations [14]. For reasonable m this complexity was substantially improved to (1) O m2:5 (n2 + m) ln mR operations by Nesterov and Nemirovskii [12], using an interior-point algorithm with a specialized \rescaling" technique to lower the work required on each iteration. A further reduction to mR n lnR 3 : 5 O m ln ln (2) operations was achieved by Khachiyan and Todd [9], who apply an interior-point algorithm to a sequence of problems, each of which requires less work per iteration than the original problem considered by [12]. Nemirovskii [11] lowers the complexity of obtaining an e? {maximal inscribed ellipsoid to (3) O m3:5 ln mR operations. The approach taken in [11] uses Lagrangian duality to reformulate the MVIE problem as a saddlepoint problem, and shows that the self-concordance theory developed in [12] can be adapted to analyze an algorithm for computing an approximate saddlepoint. The advantage of [11] is that the theory developed for saddlepoint problems is very general. However, the analysis required to develop this theory is quite extensive. In this paper we devise an optimization algorithm whose complexity to obtain an e? {maximal inscribed ellipsoid is also (3). Our work can be viewed as a further improvement of the previous optimization-based results of [12] and [9], and provides an alternative to the saddlepoint-based approach of [11]. We also show that by rst computing an approximation of the analytic center of P the eect of the parameter R can be further reduced, resulting in a total complexity of O (mn2 + m1:5n) ln(R) + m3:5 ln m (4) operations. The dierence between (3) and (4) is certainly of interest, since under standard assumptions bounds on R may be exponential in n [3, Lemma 3.1.25]. Several novel formulations of the MVIE problem are considered in [16]. Primaldual algorithms based on two of these formulations are analyzed in [17], and their numerical performance is compared to original and modi ed versions of the algorithm from [9] on instances up to size n = 500, m = 1200. The performance of one of the primal-dual algorithms is found to be superior to the other three methods. A problem related to that of computing an e? {maximal inscribed ellipsoid for P is that of computing an e {minimal circumscribing ellipsoid for the convex hull of m given points in 0. Do Until tk tmax (Outer iteration) t = tk+1 = (1 + )tk , x = xk , X = Xk . Do Until t(x; x; X) 0:15 (Inner iteration) p = pt(x; x; X), P = Pt(x; x; X) x = x + (=2)p, X = X + P.
End xk+1 = x; Xk+1 = X, k = k + 1 End
The steplength on each inner iteration can be taken to be any value that produces at least the descent in Ft(; ; ) obtained using = 1=(1 + t (x; x; X)); see Lemma 4.3 below. The use of =2 in the step for x may seem strange at rst sight, but is a simple consequence of the use of the direction p = pt (x; ; ) based on Ft (x; ; ). In particular, note that for each i, (bi ? aTi x)(bi ? aTi (x + p)) = (bi ? aTi x)2 ? (aTi p)(bi ? aTi x); 2 T 2 bi ? aTi (x + 2 p) = (bi ? aTi x)2 ? (aTi p)(bi ? aTi x) + 2 (ai4p) : As a result, using x+ = x + (=2)p, X + = X + P in a step from (x; x; X) to (x+ ; x+ ; X + ) has the same rst-order eect on Ft (; ; ) as using x++ = x + p and a step to (x; x++ ; X + ). Our analysis of the main stage algorithm for MVIE is based on the well-known analysis of the barrier algorithm from [12] (see also [4]). The following result facilitates the use of directions based on the family of barrier functions Ft(y; ; ).
6
K. M. ANSTREICHER
Lemma 4.1. Suppose that x and y are interior points of P , and let f(X) = ? ldet X , G(y; x) = fX 0 j aTi Xai (bi ? aTi y)(bi ? aTi x); i = 1; : : :; mg. De ne (y; x) = minX 2G(y;x) f(X), and t(y; x) = minX 2G(y;x) Ft(y; x; X), t 1. Then 1. (y; x) 12 [(y; y) + (x; x)], 2. t(y; x) 21 [t(y; y) + t (x; x)].
Proof. Part 1 is proved in [9, p.144], and we use a similar argument to prove part 2 here. Assume that X 0 is the minimizer of Ft(x; x; ), and Y 0 is the minimizer of Ft(y; y; ). By a change of coordinates we may assume without loss of generality that X and Y are diagonal. For each i = 1; : : :; m the inequalities aTi Xai (bi ? aTi x)2 and aTi Y ai (bi ? aTi y)2 together imply that (10) aTi (XY )1=2 ai kX 1=2ai k kY 1=2ai k = [(aTi Xai )(aTi Y ai)]1=2 (bi ? aTi x)(bi ? aTi y); and ldet(XY )1=2 = (1=2)(ldet X + ldet Y ). To prove that t(y; x) 12 [t(y; y) + t(x; x)] it then suces to show that for each i = 1; : : :; m,
? ln (bi ? aTi x)(bi ? aTi y) ? aTi (XY )1=2ai
? 21 ln (bi ? aTi x)2 ? aTi Xai ? 12 ln (bi ? aTi y)2 ? aTi Y ai ; ?
?
which is equivalent to
2
(bi ? aTi x)(bi ? aTi y) ? aTi (XY )1=2 ai ? ? (bi ? aTi x)2 ? aTi Xai (bi ? aTi y)2 ? aTi Y ai : Using the rst inequality in (10), to prove (11) it suces to show that (11)
2
(bi ? aTi x)(bi ? aTi y) ? [(aTi Xai )(aTi Y ai )]1=2 ? ? (bi ? aTi x)2 ? aTi Xai (bi ? aTi y)2 ? aTi Y ai ;
which reduces to [(bi ? aTi x)(aTi Y ai )1=2 ? (bi ? aTi y)(aTi Xai )1=2]2 0. We now use Lemma 4.1, and standard results on self-concordant functions, to bound the possible reduction in Ft(; ; ) and f() when the Newton decrement is suciently small. Lemma 4.2. Suppose that t > 1, x 2 Int(P ), (x; X) 2 Int(G(x)), and t (x; x; X) 1=3. Let = () = 2 =(1 ? 5:06252). Then 1. Ft(y; y; Y ) Ft(x; x; X) ? 2 for all (y; Y ) 2 G. 2. f(X) fmin + [12m + 2]=(t ? 1). Proof. From part 1 of Lemma 2.2 we have (12) Ft (x; y; Y ) Ft(x; x; X) ? for any (y; Y ) 2 G(x). Part 2 of Lemma 4.1 then implies that for any y 2 Int(P ), t (y; y) 2t (y; x) ? t(x; x) 2(Ft (x; x; X) ? ) ? Ft(x; x; X) = Ft (x; x; X) ? 2;
MAXIMUM VOLUME INSCRIBED ELLIPSOIDS
7
which proves part 1. Next, note that Ft(x; x; X) = (t ? 1)f(X) + F1(x; x; X), where F1(x; ; ) is a (2m + n)-self-concordant barrier for G(x). From (12) and a standard argument (see for example [12, p.75]) we conclude that for all (y; Y ) 2 G(x), + ; f(Y ) f(X) ? 2(2mt +? n) 1 and therefore (x; y) f(X) ? (6m + )=(t ? 1) for all y 2 Int(P ), since m > n. Part 1 of Lemma 4.1 then implies (y; y) 2(y; x) ? (x; x) + ? f(X) 2 f(X) ? 6m t?1 12m = f(X) ? t ?+12 ; which proves part 2. The next lemma considers the eect of increasing t on the Newton decrement for Ft(x; x; X), and the reduction in Ft(; ; ) that can be assured if the Newton decrement is not suciently small. Lemma 4.3. For t 1 let t(y; x; X) be the Newton decrement for Ft(y; ; ) at (x; X). 1. Suppose that p t(x; x; X) :15, and let t+ = (1 + )t, 0. Then t+ (x; x; X) :15(1 + ) + 2m. 2. Suppose that t (x; x; X) = > :15. Let = 1=(1+), x+ = x+(=2)pt(x; x; X), X + = X + Pt(x; x; X). Then (x+ ; X + ) 2 Int(G), and Ft (x+ ; x+ ; X + ) Ft(x; x; X) ? 0:01. Proof. Part 1 is proved in [4, Lemma 2.25]. Since Ft(x; ; ) is strongly selfconcordant for t 1, part 2 of Lemma 2.2 implies that if x++ = pt(x; x; X), then (x++ ; X + ) 2 Int(G(x)), and (13)
Ft (x; x++ ; X + ) Ft(x; x; X) ? ( ? ln(1 + )) Ft (x; x; X) ? 0:01:
It is also easy to see that for any x and y in P , and i = 1; : : :; m, (bi ? aTi x)(bi ? aTi y)
bi ? aTi (x +2 y)
2
;
and therefore if (y; Y ) 2 G(x), then (14)
Ft x +2 y ; x +2 y ; Y Ft(x; y; Y ):
Part 2 is completed by combining (13) and (14), with (y; Y ) = (x++ ; X + ). We can now combine Lemmas 4.2 and 4.3 to obtain the nal complexity result for the main stage. p Theorem 4.4. Let = 0:07= m. Then the main stage algorithm requires O(1) inner iterations per outer iteration. Moreover, for tmax = O(m=) the algorithm terminates with an e? {maximal ellipsoid in O(m:5 ln(m=)) outer iterations, requiring a total of O(m3:5 ln(m=)) operations.
8
K. M. ANSTREICHER
Proof. From part 1 of Lemma 4.3, at the start of each sequence of inner iterations
we have
tk+1 (xk ; xk; Xk ) :15 1 + p:07m + :1 < 0:26: Part 1 of Lemma 4.2 then implies that for any (x; X) 2 G, Ftk+1 (x; x; X) > Ftk+1 (xk ; xk ; Xk ) ? :21; and from Part 2 of Lemma 4.3 there can be at most 20 inner iterations on each outer iteration. From part 2 of Lemma 4.2, to obtain an e? {maximal inscribed ellipsoid it suces to terminate the algorithmpusing tmax = O(m=), which requires O(m:5 ln(m=)) outer iterations using = :07= m. Finally, it can be shown using a small modi cation of the argument used in [9, x6] that each inner iteration can be executed in O(m3 ) operations. To close this section we note two details regarding the de nition of Ft (y; x; X), from (9). First, it would be more \standard" to replace the term ?t ldet X with ?(t + 1) ldet X, and work with t 0 instead of t 1. However, in this case the main stage algorithm must be initialized at a suciently large t0 > 0, which would unnecessarily complicate the analysis of the preliminary stage in the next section. Second, the analysis of this section does not require the term (15)
?
m X i=1
?
ln (bi ? aTi y)(bi ? aTi x) ;
in (9), and in fact this term slightly degrades the theoretical performance of the main stage algorithm. However (15) is very helpful for the analysis of the preliminary stage. 5. Preliminary stage. In this section we consider the preliminary stage of our barrier algorithm for the MVIE problem. The goal of the preliminary stage is to produce the initial point (x0 ; X0) required by the main stage algorithm of the previous section. Our preliminary stage is based on the general preliminary stage described in [12, x3.2.3], except that we work with directions based on F1 (y; x; X) so as to keep the work per iteration O(m3 ). The preliminary stage is initialized at x0 = 0, and a suitable X0 to be described below. Let m b X 1 i T (16a) + ai ; c = ?rx F1(0; 0; X0) = ? i=1 i;0 bi m X (16b) C = ?rX F1 (0; 0; X0) = X0?1 ? 1 ai aTi ; i=1 i;0
where i;0 = b2i ? aTi X0 ai . For t 1 de ne the preliminary stage barrier function Ft0(y; x; X) = t(cT (x + y) + C X) + F1(y; x; X): Let [p0t (y; x; X); Pt0(y; x; X)] denote the Newton direction for Ft0(y; ; ) at (x; X), and 0t (y; x; X) the corresponding Newton decrement. Note that by construction 01 (0; 0; X0) = 0. The preliminary stage algorithm, given below, is very similar to the main stage, except that the preliminary stage uses Ft0(; ; ) and t is decreased rather than increased on each outer iteration.
9
MAXIMUM VOLUME INSCRIBED ELLIPSOIDS
ALGORITHM (Preliminary stage for MVIE): Given k = 0, x0 = 0, X0 , t0 = 1, tmin , > 0. Do Until tk tmin (Outer iteration) t = tk+1 = (1 ? )tk , x = xk , X = Xk . Do Until 0t (x; x; X) 0:15 (Inner iteration) p = p0t (x; x; X), P = Pt0(x; x; X) x = x + (=2)p, X = X + P.
End xk+1 = x; Xk+1 = X, k = k + 1 End
To analyze the preliminary stage we require an extension of part 2 of Lemma 4.1 that applies to Ft0(y; x; X). This turns out to be straightforward under the assumption that C 0. Lemma 5.1. For interior points x and y of P , let 0t (y; x) = minX 2G(y;x) Ft0(y; x; X), 0 < t 1. Assume that C 0. Then 0t (y; x) 21 [0t (y; y) + 0t (x; x)]. Proof. The proof is identical to the proof of part 2 of 4.1. Note that the change of coordinates that simultaneously diagonalizes X and Y preserves the semide niteness of C. After this change of coordinates we have C X =
m X i=1
CiiXii ; C Y =
m X i=1
CiiYii ; C (XY )1=2 =
m X i=1
p
Cii Xii Yii;
where Cii 0, i = 1; : : :; m. It immediately follows that C (XY )1=2 (1=2)(C X + C Y ). For C 0 and a given value of tmin , the analysis of the preliminary stage is very similar to the analysis of the main stage given in the previous section, and is omitted here. The nal complexity result has the following form. p Theorem 5.2. Assume that C 0. Then for = = m, where > 0 is an appropriate positive constant, the preliminary stage requires O(m:5 ln(1=tmin )) outer iterations, O(1) inner iterations per outer iteration, and a total of O(m3:5 ln(1=tmin )) operations.
To complete the analysis of the preliminary stage we must show that X0 can be chosen so that C 0, and characterize the value tmin so that termination of the preliminary stage produces a suitable initial point for the main stage. Lemma 5.3. Assume that B(0; 1) P , and let X0 = m1+1 I . Then C 0. Proof. By the assumption that B(0; 1) P we must have bi kaik, i = 1; : : :; m. Then for each i, (17) i;0 = b2i ? aTi X0 ai = b2i ? m 1+ 1 kaik2 b2i mm+ 1 ; and !
m m 1 X X T C (m + 1)I ? mm+ 1 b12 ai aTi (m + 1) I ? m1 2 aiai 0: k a k i i=1 i i=1 follows from (16b). It remains to obtain a lower bound on the required value of tmin . To this end, note that for any strongly self-concordant function F() de ned on the interior of a
10
K. M. ANSTREICHER
compact, convex set G