Coursehero >> New York >> Cornell >> VIVO 22892

Course Hero has millions of student submitted documents similar to the one below including study guides, homework solutions, papers, exam answer keys and textbook solutions.

for Algorithms Ellipsoids Stephen B. Pope Sibley School of Mechanical & Aerospace Engineering Cornell University Ithaca, New York 14853 Report: FDA-08-01 February 2008 Abstract We describe a number of algorithms to perform basic geometric operations on ellipsoids in n spatial dimensions, for n 1. These algorithms are implemented in ELL LIB, a library of Fortran subroutines. With E, E1 and E2 being given ellipsoids, and p a given point, the tasks considered include: determine the point in E which is closest to p or furthest from p; grow or shrink E so that its boundary intersects p; project E onto a given a ne space; determine a separating hyperplane between E1 and E2 ; determine an ellipsoid (of small volume) which covers E1 and E2 . Contents 1 Introduction 2 Representation of ellipsoids 3 Summary of routines 4 Useful preliminary results 4.1 Linear transformation . . 4.2 Quadratic minimization 4.3 Householder matrix . . . 4.4 Rank-one modi cation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3 5 5 6 8 8 9 10 12 12 13 14 E and .... .... .... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p ........... ........... ........... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 14 17 17 18 19 19 19 20 21 22 22 24 27 5 Smallest and largest principal semi-axes of E 6 Is the point x covered by E? 7 Relative distance to the boundary of E 8 Nearest point in E to a given point 9 Furthest point in E to a given point 10 Minimum-volume ellipsoid covering 10.1 Householder matrix algorithm . . . 10.2 Rank-one modi cation algorithm . 10.3 Behavior . . . . . . . . . . . . . . . 11 Shrink E based on a given point 11.1 Maximum-volume algorithm . . . . 11.1.1 Algorithm . . . . . . . . . . 11.1.2 Behavior . . . . . . . . . . . 11.2 Near-content algorithm . . . . . . . 11.2.1 Algorithm . . . . . . . . . . 11.2.2 Behavior . . . . . . . . . . . 11.3 Conservative algorithm . . . . . . . 11.3.1 Algorithm: reduction to 2D 11.3.2 Algorithm: solution in 2D . 1 12 Orthogonal projection of E onto a given line 13 Orthogonal projection of E onto an a ne space 29 30 14 Generate an ellipsoid which does not cover any speci ed points 32 15 Separating hyperplane of two ellipsoids 16 Pair covering query 34 36 17 Shrink ellipsoid so that it is covered by a concentric ellipsoid 36 18 Ellipsoid that covers two given ellipsoids 18.1 Spheroid algorithm . . . . . . . . . . . . 18.2 Covariance algorithm . . . . . . . . . . . 18.3 Iterative algorithm . . . . . . . . . . . . 18.3.1 Stage 1 . . . . . . . . . . . . . . . 18.3.2 Stage 2 . . . . . . . . . . . . . . . 18.3.3 Stage 3 . . . . . . . . . . . . . . . 18.3.4 Stage 4 . . . . . . . . . . . . . . . 18.3.5 Stage 5 . . . . . . . . . . . . . . . 18.3.6 Stage 6 . . . . . . . . . . . . . . . 18.3.7 Mutual covering . . . . . . . . . . 18.3.8 Discussion . . . . . . . . . . . . . 19 Conclusions 20 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 38 39 41 41 41 41 42 44 45 45 46 46 46 2 1 Introduction In this paper we describe a number of algorithms to perform basic geometric operations on ellipsoids. These algorithms have been implemented in Fortran routines which are contained in the library ELL LIB, which is available at http://eccentric.mae.cornell.edu/ tcg/ELL LIB. Ellipsoids arise in numerous computational problems. The algorithms and routines described here have been developed speci cally for use in the in situ adaptive tabulation (ISAT) algorithm (Pope 1997). In Sec. 2 we consider di erent mathematical representations of an ellipsoid E in n , for n 1. For n = 1, E is a line segment; for n = 2, E is an ellipse; for n = 3, E is an ellipsoid; and for n > 3, E is a hyper-ellipsoid. For simplicity we generally refer to E (for all n 1) as an ellipsoid. A summary of the routines in ELL LIB is provided in Sec. 3. Following some preliminary results in Sec. 4, the algorithms used are described in Secs. 5 18. 2 Representation of ellipsoids There are many ways to represent ellipsoids, with the di erent ways arising naturally in di erent circumstances. In this section we show the relations between the di erent representations. Let the ellipsoid E be centered at c; let the columns of the n n orthogonal matrix U be unit vectors in the directions of E s principal axes; and let be the diagonal matrix (with diagonal elements ii = i ) such that 1/ i is the length of the ith principal semi-axis. We assume that the principle axes are nite and strictly positive, i.e., 0 < i < . Then E is given by E {x | (x c)T U 2 UT (x c) 1}. This may alternatively be expressed as E {x | UT (x c) 1}. or, from the de nition w UT (x c), E {x | x = c + U 1 w, w 1}. (3) (2) (1) The above three de nitions can be re-expressed in terms of di erent related matrices. Let A be the matrix 3 A U 2 UT , (4) appearing in Eq.(1). Evidently A is symmetric positive de nite; its eigen2 vectors are the columns of U and its eigenvalues are i = i . We denote by = 2 the diagonal matrix of eigenvalues. Thus Eqs. (1) (3) can be 1 trivially re-written by substituting 2 for ; or less trivially in terms of A as E {x | (x c)T A(x c) 1}, (5) E {x | A 2 (x c) 1}, E {x | x = c + A 2 w, 1 1 (6) (7) w 1}. Let B be a non-singular square matrix, which we use to form A as A = BBT , and let the SVD of B be B = U VT . Note that we have from Eq.(9), A = BBT = U 2 UT , (10) (9) (8) consistent with Eq.(4), and showing that there is a family of matrices B yielding the same matrix A, namely B = U VT for given U and but arbitrary orthogonal V. In terms of B, Eqs. (1) (3) can be reexpressed as E {x | (x c)T BBT (x c) 1}, E {x | BT (x c) 1}, E {x | x = c + B T w, with a di erent de nition of w. The matrix B can also be factored as B = LQ, 4 (14) w 1}, (11) (12) (13) where L is lower triangular with positive diagonal elements and Q is orthogonal. Thus we obtain A = BBT = LLT , (15) showing that L is the Cholesky factorization of A. The de nitions of E in terms of B apply equally in terms of L, i.e., E {x | (x c)T LLT (x c) 1}, E {x | LT (x c) 1}, E {x | x = c + L T w, w 1}. (16) (17) (18) Computationally it is most e cient to use the Cholesky representation of the ellipsoid in terms of c and L, and to store L in packed format. In ELL LIB, routines with names starting ell use this representation, while those with names starting ellu represent L in unpacked format. In the algorithms described below, it often happens that an ellipsoid E1 given by c1 and L1 is modi ed to yield an ellipsoid E2 which is known in terms of c2 and B2 (i.e., a full matrix). The representation in terms of L2 is e ciently and accurately computed by the LQ algorithm. It is stressed that it is much more accurate to obtain L2 from the LQ algorithm than from the Cholesky decomposition of A = B2 BT . All of the algorithms described 2 below use the LQ algorithms and avoid forming A. The following table shows various routines in ELL LIB that can be used to transform between one representation of an ellipsoid and another. In this table, A denotes the lower triangle of the symmetric matrix A, while all other symbols have the meanings given above. 3 Summary of routines Table 2 summarizes the tasks performed by the principal routines in ELL LIB. Here E, E1 and E2 denote given ellipsoids, and p denotes a given point. 4 Useful preliminary results In this Section we give some general results that are used in the subsequent sections. 5 Table 1: Routines to transform between Routine From ell bbt2chol B ellu bbt2chol B ell bbt2eig B ell chol2eig L ellu chol2eig L ell eig2chol (U, ) ell full2eig A ell full2low A ell low2chol A di erent representations To L L (U, ) (U, ) (U, ) L (U, ) A L 4.1 Linear transformation y = MT (x b), (19) We consider the linear transformation where b is a speci ed vector and M is a speci ed non-singular n n matrix. The inverse transformation is x = b + M T y. If the ellipsoid E is de ned by E {x | BT (x c) 1}, then it is readily shown, from Eq.(20), that it is also given by T E = {y | B (y c) 1}, with and c MT (c b). (24) B M 1 B, (22) (23) (21) (20) 6 Table 2: Tasks performed by principal routines in ELL LIB Routine Section Task performed ell rad lower 5 Determine the smallest principal semi-axis of E ell rad upper 5 Determine the largest ellu rad upper principal semi-axis of E ell radii upper 5 Determine the smallest and largest ellu radii upper principal semi-axes of E ell pt in 6 Determine whether p is covered by E ell pt dist 7 Determine the relative distance from p to the boundary of E ell pt near far 8 Determine the point in E which is ellu pt near far 9 nearest to or furthest from p ell pt modify 10.2 Determine EV , the modi cation of E of least change of content whose boundary intersects with p ell pt shrink 11 Shrink E based on p ell pt hyper Determine a separating hyperplane between E and p ell pts uncover 14 Determine an ellipsoid (of large content) which does not cover any given point ell line proj 12 Project E onto a given line ell a pr 13 Project E onto a given a ne space ell pair shrink 17 For concentric E1 and E2 , shrink E2 (if necessary) so that it is covered by E1 ell pair separate 15 Determine if E1 and E2 intersect; and, if they do not, determine a separating hyperplane ell pair cover query 16 Determine wether E1 covers E2 ell pair cover 18 Determine an ellipsoid (of small content) which covers E1 and E2 7 4.2 Quadratic minimization Several of the algorithms described below depend on the solution to the following problem: determine a vector x which minimizes the quadratic function 1 g(x) 2 xT Ax + bT x, (25) (26) subject to xT x 2 , where , b and A are a given scalar, vector, and symmetric matrix, respectively. Such problems are e ciently solved by the routine dgqt from MINPACK-2 (see Averick et al. 1993). 4.3 Householder matrix Given a vector p with p p > 0, the corresponding Householder matrix H(p) is de ned by H = I 2vvT , (27) where v is the unit vector v1 = [(1 + |p1 |)/(2p)] 2 , sign (p1 )pi vi = , for i 2. 2v1 p 1 (28) It is readily shown that H has the following properties (not all independent): 1. H is symmetric: HT = H. 2. H is orthogonal: HT H = I. 3. The rst column of H is parallel (or anti-parallel) to p. 4. All columns of H except the rst are orthogonal to p. 5. H maps p to the rst axis: Hp = pe1 . Given p, the routine ell house returns v. 8 4.4 Rank-one modi cation Let E be an ellipsoid centered at the origin (i.e., c = 0), given in terms of the positive symmetric de nite matrix A, which has SVD A = U 2 UT , and Cholesky decomposition A = LLT . For a given vector w and scalar , let F be the rank-one modi cation of A: F = A + wwT . (29) We require F to be positive symmetric de nite, which in turn requires to exceed a critical value 0 ( 0 < 0) above which all eigenvalues of F are positive (see item 7 below). Then, the modi ed ellipsoid E is de ned by E {x | xT Fx 1}. The following results are readily obtained: 1. For 0, xT Fx is greater than or equal to xT Ax, and hence E covers E. 2. For 0, xT Fx is less than or equal to xT Ax, and hence E covers E. 3. The eigenvalues of F are interlaced with those of A (see Golub and Van Loan 1996, Sec. 8.5.3). Hence the lengths of the principal axes of E and E are also interlaced. 4. With z UT x, we have xT Ax = zT z, so that E is the unit ball in z-space. Correspondingly, E is given by E = {z | zT (I + wwT )z 1}. where w = 1 UT w. Thus in z-space, one principal semi-axis of E is (1 + |w|2 ) and the others are orthogonal unit vectors. 1 2 (30) (31) (32) w/|w|, 5. Similarly, with y LT x, we have xT Ax = yT y, so that E is the unit ball in y-space. Correspondingly, E is given by E = {y | yT (I + wwT )y 1}. where w = L 1 w. 9 (34) (33) 6. The rank-one modi cations to the identity appearing in Eq.(31) has the symmetric square root I + wwT = G2 = (I + wwT )(I + wwT ) where and are related by = 2 + 2 |w|2 , = ( {1 + |w|2 } 2 1)/|w|2 . 1 (35) (36) (37) The matrix in Eq.(33) has a similar square root, with the same value of , since |w| = |w|. 7. The critical value of is 0 = 1/|w|2 , corresponding to the smallest value of for which is real. (38) 5 Smallest and largest principal semi-axes of E As sketched in Fig. 1, given an ellipsoid E, the task is to determine its smallest and largest principal semi-axes, or equivalently, the radii of the inscribed and circumscribed hyper-spheres, rin and rout , respectively. This task is accomplished by the routines ell rad lower and ell rad upper. Considering the ellipsoid E, Eq.(16), let the SVD of the Cholesky factor L be L = U VT , (39) with 1 2 . . . n being the components of the diagonal matrix of singular values . Then the matrix A, Eq.(15), is A = LLT = U 2 UT = U UT . (40) Evidently the columns of U are eigenvectors of A, and the diagonal matrix consists of the eigenvalues of A with components 2 i = i . (41) 10 E rin rout Figure 1: Sketch of an ellipsoid E showing the radii, rin and rout of the inscribed and circumscribed hyper-spheres. The principal semi-axes are given by ri = i principle semi-axis is 1 rin = [max( i )] 2 , and the largest is rout = [min( i )] 2 . 1 1 2 1 = i . Thus the smallest (42) (43) The most stable way to compute rin is as rin = 1/ 1 , where 1 is obtained from the SVD of L, and similarly for rout . The routine ell radii determines both rin and rout using the SVD. The routines ell rad lower and ell rad upper determine rin and rout at signi cantly lower computational cost, but with less accuracy. Given L, the routine ell rad lower determines rin via Eq.(42), using the LAPACK routine dsyevx to compute the largest eigenvalue of A. The routine ell rad upper determines rout via quadratic minimization. Since rout is the furthest distance from the center to any point on the ellipsoid we 2 have that rout is the maximum of rT r subject to rT Ar 1. (44) By de ning y = LT r, this can be put in the standard form of a quadratic 2 minimization problem: rout is the minimum of yT L 1 L T y, 11 (45) subject to yT y 1. (46) 6 Is the point x covered by E? Given a point x and an ellipsoid E (in terms of c and L), the task is to determine whether E covers x. This is readily determined (by the routine ell pt in) through the de nition of E given by Eq.(17). That is, E covers x if the quantity s LT (x c) , (47) is less than or equal to unity. 7 Relative distance to the boundary of E Given an ellipsoid E (in terms of c and L) and a point p (p = c), let b be the intersection of the ray c-p with the boundary of E. The relative distance s to the boundary is de ned by: s |b c|/|p c|. Now the boundary point satis es LT (b c) = 1, and we have b c = s (p c). Hence, s is determined as s = LT (p c) 1 (48) (49) (50) . (51) The routine ell pt dist determines s. It may be noted that the three cases s < 1, s = 1 and s > 1 correspond, respectively, to: p not being covered by E; p being on the boundary of E; and, p being covered by E. 12 f E n p Figure 2: The points n and f in the ellipsoid E which are nearest and furthest, respectively, to the given point p. 8 Nearest point in E to a given point We consider the ellipsoid E Eq.(16) centered at c, and a given point p. The task (performed by the routine ell pt near far)is to nd the point x in E closest to p (see Fig.2). Thus x satis es (x c)T LLT (x c) 1, and minimizes s2 (x p)T (x p). (53) Evidently, if p is contained in E then x = p and s = 0. Otherwise x is on the boundary of E. The above equations are readily transformed into the standard quadratic minimization problem. Let y be de ned by y LT (x c), so that Eq.(52) becomes yT y 1. Equation (54) can be inverted to yield x = c + L T y, so that Eq.(53) becomes s2 = (c + L T y p)T (c + L T y p) = yT L 1 L T y + 2(c p)T y + (c p)T (c p). 13 (57) (56) (55) (54) (52) Thus y is obtained via dgqt by minimizing 1 yT L 1 L T y + (c p)T y, 2 subject to y 1; and then x is obtained from Eq.(56). (58) 9 Furthest point in E to a given point This task is essentially the same as that considered in the previous section, except that the furthest point in E, x, maximizes s2 x p 2 , Eq.(53) (rather than minimizing s2 ). Thus the same algorithm is used, but the quantity minimized is , Eq.(58) (rather than ). Again, this is performed by the routine ell pt near far. For all given points p, the corresponding furthest point x lies on the boundary of the ellipsoid E (see Fig.2). 10 Minimum-volume ellipsoid covering E and p Given an ellipsoid E and a point p, are seek the ellipsoid E (concentric with E) of minimum volume which covers both E and p. In the rst subsection we describe an algorithm based on the geometry of the problem. In the subsequent subsection a simpler algorithm is given, which is implemented in the routine ell pt modify. Then the behavior of the algorithm is examined. 10.1 Householder matrix algorithm As sketched in Fig. 3(a), consider a given point p lying outside an ellipsoid E. A modi ed concentric ellipsoid E is sought which has minimum volume subject to: 1) E intersects p 2) E covers E. This is achieved by 1) performing the linear transformation which transforms E to the unit hypersphere, and p to a point p on the rst axis, see Fig. 3 (b) 14 c E c p E Figure 3: Sketch of: (a) the given ellipsoid E and point p; (b) the transformation of E to the unit ball, and p to a point on the rst axis; (c) the transformed ellipsoid E ; and (d) the grown ellipsoid E . 2) extending the hypersphere in the rst direction to form an ellipsoid with p on its boundary, see Fig. 3 (c) 3) inverting the initial transformation to obtain E , see Fig. 3 (d). The ellipsoid E (Eq.(17)) is de ned by LT r = LT (x c) 1. De ning the transformed variable y by y = LT (x c), evidently the ellipsoid in y-space is the unit hypersphere y 1. The same transformation applied to p yield p = LT (p c). 15 (62) (61) (60) (59) (c) (a) p (b) p (d) p There is an orthogonal matrix Q which transforms p to a point on the 1-axis; that is, p = QT p = QT LT (p c), (63) where pi = p i1 , (64) and p = p . It is readily shown that Q is simply the Householder matrix Q = H( ). The same transformation applied to E yields the unit ball p zT z 1, where z = QT LT (x c). (66) p This ball and the point p ( i p i1 ) are shown in Fig. 3 (b). The next step is to de ne the modi ed ellipse E in z-space, as shown in Fig. 3 (c). This ellipsoid is zT 2 z 1, (67) where 2 = ij + i1 j1 ( 2 1). p ij (68) Transforming back to the original space, we obtain the equation for E (see Fig. 3 (d)): (x c)T LQ 2 QT LT (x c) 1, (69) or (x c)T L L T (x c) 1. (70) Thus the modi ed ellipsoid E has center c and Cholesky matrix L given by L L T = LQ 2 QT LT = (LQ )(LQ )T . (65) (71) In summary, given the ellipsoid E (in terms of c and L) and the point p, then p is determined by Eq.(62), v by Eq.(28), Q by Eq.(27), 2 by Eq.(68), and nally the Cholesky matrix of the modi ed ellipsoid E is determined by Eq.(71). 16 10.2 Rank-one modi cation algorithm Following from Eqs.(61) and (62), in y-space, the modi ed ellipsoid E is the unit ball extended to intersect the point p. Thus, based on the results of Sec. 4.4, we have E = {y | yT G2 y 1}, (72) where G = I + ppT , and (determined by the condition pT G2 p = 1) is = 1 1 1 . | | p | |2 p (74) (73) Thus, in place of Eq.(71), the Cholesky matrix L can be obtained as L L T = (LG)(LG)T . This algorithm is implemented in the routine ell pt modify. (75) 10.3 Behavior Test are reported for the minimum-volume growing algorithm. The tests are in two dimensions and are based on an initial ellipsoid E (centered at the origin and aligned with the coordinate axes). The length of the major semiaxis is rmajor = 1, and the minor semi-axis is rminor . The grow point p is not covered by E, and the vector p is at an angle to the x1 axis. An example of the grown ellipse E given by the minimum-volume algorithm is shown in Fig. 4. The lengths of the principal semi-axes of the grown ellipse E are denoted by Rmajor and Rminor . A gure of demerit F of the growth operation is de ned by: F Rmajor / max(rmajor , |p|) 1. (76) Tests suggest the following behavior: 1. For given E and , the maximum of F occurs for |p| = rmajor , i.e., for p being on the bounding circle. 2. As illustrated in Fig. 5, for given E, F has a unique maximum (Fmax ) at = max between = 0 and = /2. 17 Minimum volume growing for grow point at 70 degrees 1 0.5 0 0.5 1 1.5 1 0.5 0 0.5 1 1.5 Figure 4: Illustration of the minimum-volume algorithm. The initial ellipse E (solid blue line); its bounding circle (dashed blue line); the grow point p (magenta); and the grown ellipse (red line) E: rminor = 0.1, |p| = rmajor = 1, = 70o . 3. rminor decreases, max tends towards /2, and Fmax tends towards As 2. With some justi cation, it is speculated that even in higher dimensions F is bounded above by 2. 11 Shrink E based on a given point Given an ellipsoid E and a point p, the problem considered is to generate a concentric ellipsoid E of smaller volume which covers p. There is no unique solution to this problem. In the following three subsections we describe three algorithms which yield di erent solutions for E , denoted by EV , EN and EC , respectively. The third (which is based on the rst two) is preferred in application to ISAT. The three algorithms are implemented in the routine ell pt shrink with the parameter k ell set to 1, 2 and 3, respectively. 18 Amplification of the major axis against the angle of the grow point 1.4 1.35 1.3 1.25 Rmajor/rmajor 1.2 1.15 1.1 1.05 1 0.95 0 10 20 30 40 50 60 70 80 90 Figure 5: The demerit measure F as a function of the grow-point angle : rminor = 0.1, |p| = rmajor = 1. 11.1 11.1.1 Maximum-volume algorithm Algorithm The algorithm described in section 10 to grow on ellipsoid E based on a point p lying outside of E can also be applied when the point p lies inside E. The resulting modi ed ellipsoid EV has maximal volume subject to: 1) p is on the boundary of EV 2) EV is covered by E. 11.1.2 Behavior We consider the same 2D case as for the minimum-volume growing algorithm, except that here the shrink point p is covered by the ellipse E. The maximumvolume shrinking algorithm determines the ellipse EV of maximum volume which is covered by E and has p on its boundary. This is illustrated in Fig. 6. We de ne a gure of merit F Rminor / min(rminor , |p|) 1, (77) 19 Figure 6: The maximum-volume shrinking algorithm. The initial ellipse E (solid blue line); its inscribed circle (green line); the grow point p (magenta); and the maximum-volume shrunk ellipse EV (red line): rminor = 0.4, rmajor = 1, = 70o , |p| = 0.48. Figure 7: The maximum-volume shrinking algorithm. The initial ellipse E (solid blue line); its inscribed circle (green line); the grow point p (magenta); and the maximum-volume shrunk ellipse EV (red line): rminor = 0.1, |p| = rmajor = 1, = min = 84.3o . where rminor and Rminor are the lengths of the minor semi-axes before and after shrinking. Based on empirical testing we draw the following conclusions. 1. For given E and , the minimum (Fmin ) of F occurs for p = rminor , i.e., for p lying on the inscribed circle. 2. As rminor /rmajor decreases, the values of at which F is minimum (i.e., F ( min ) = Fmin ) tends to /2 rminor /rmajor , and Fmin tends to zero as 2rminor /rmajor . Figure 7 shows the case rminor /rmajor = 0.1, = min = 84.3o for which Fmin = 0.1962. 11.2 Near-content algorithm As shown above, the maximum-volume shrinking algorithm can lead to very small gures of merit F , Eq.(77). As may be seen from Fig. 7, the shrunk 20 ellipsoid can exclude a substantial portion of the original ellipsoid which is closer to the center than the shrink point. Here we describe the alternative near-content shrinking algorithm which yields values of F of one, or close to one. However, unlike the maximum-volume algorithm, the result depends on the metric of the space. It is implemented in the routine ell pt shrink (for the parameter value k ell= 2). 11.2.1 Algorithm We rst describe the algorithm and then give its partial justi cation. The initial ellipsoid E centered at c is given in terms of the matrix A with SVD A = U 2 UT . We transform to principal axes (y-space) by de ning y = UT x, and p = UT (p c). (79) In y-space, the modi ed ellipsoid EN is de ned as the rank-one modi cation to E: EN = {y | yT Fy 1}, (80) F = 2 + wwT , w = D , p where the diagonal matrix D is de ned by Dii = max(0, 1/| |2 2 ), p ii and the positive scalar is determined by the intersection condition p pT F = 1. We make the following observations about the algorithm: 1. The quantity 1/| |2 2 (appearing in Eq.(83)) is positive if, and only p ii if, the length of the i-th principal semi-axis 1 is greater than | |. p ii (84) (83) (81) (82) (78) 21 p 2. Since E covers p, pT D is strictly positive, and hence a positive value of exists which satis es Eq.(84), namely pp = (1 pT 2 p)/( T D ). (85) 3. If p lies inside the ball of radius 1 , then Dii = 1/| |2 2 , p is an p 11 ii eigenvector of F, and hence p is the smallest principal semi-axis of EN . This is a necessary condition for EN to be a maximum nearest content ellipsoid. 4. If for some i (1 i < n), p lies between the balls of radius 1 ii 1 and (i+1)(i+1) , then wj = 0 for j i, and as a consequence the rst i principal axes of EN are the same as those of E. This again is a necessary condition for EN to be a maximum nearest content ellipsoid. 11.2.2 Behavior The behavior of the near-content algorithm is illustrated for a 2D case in Fig. 8; and it is compared to the maximum-volume algorithm in Fig. 9. There is a striking di erence between the two methods, especially for relatively small |p|. 11.3 Conservative algorithm As may be seen from the last pane in Fig.9, the shrunk ellipsoids (denoted by EV and EN ) generated by the maximum-volume and near-content algorithms can be quite di erent. In particular each excludes substantial regions included by the other. Here we describe a third algorithm which is conservative in the sense that the shrunk ellipsoid generated, EC , includes all points in EV and EN . Speci cally, EC is de ned as the ellipsoid of minimum volume which covers both EV and EN . This algorithm is implemented in the routine ell pt shrink (for the parameter value k ell= 3). The behavior of the conservative algorithm is illustrated in Fig. 10. It has the following properties: 1. EC covers p. Generally, p is in the interior of EC . 2. EC covers both EV and EN . 3. In general EC is not covered by E. 22 Max. close content: |p|/rmin=1.5 Max. close content: |p|/rmin=1.01 0.5 0.5 0 0 0.5 0.5 1 0.5 0 0.5 1 1 0.5 0 0.5 1 Max. close content: |p|/r min =0.99 Max. close content: |p|/r min =0.5 0.5 0.5 0 0 0.5 0.5 1 0.5 0 0.5 1 1 0.5 0 0.5 1 Figure 8: The near-content shrinking algorithm. The initial ellipse E (solid blue line); its inscribed circle (light blue line); the grow point p (magenta); and the maximum near content shrunk ellipse EN (red line): rmajor = 1, rminor = 0.4, = 70o . 4. The largest principal axis of EC can exceed that of E. Tests suggest that this ratio of principal axes is seldom greater than 1.3. 5. The volume of EC is no greater than that of E, and in general is less. (This follows from the fact that both E and EC cover EV and EN , but EC is, by de nition, of minimum volume.) 6. If the algorithm is applied a second time to EC based on the original point p, the results is (in general) a di erent ellipsoid. (This is in contrast to the maximum-volume and near-content algorithms in which p intersects the boundaries of EV and EN and hence a re-application of the algorithm no has e ect.) 7. Because the algorithm involves EN , the result EC depends on the metric of the space. 23 Max. close content: |p|/rmin=1.6 Max. close content: |p|/rmin=1.2 0.5 0.5 0 0 0.5 0.5 1 0.5 0 0.5 1 1 0.5 0 0.5 1 Max. close content: |p|/r min =0.6 Max. close content: |p|/r min =0.2 0.5 0.5 0 0 0.5 0.5 1 0.5 0 0.5 1 1 0.5 0 0.5 1 Figure 9: Comparison of the near-content and maximum-volume shrinking algorithms. The initial ellipse E (solid blue line); its inscribed circle (light blue line); the grow point p (magenta); the near content shrunk ellipse EN (red line); and the maximum-volume shrunk ellipse EV (green line): rmajor = 1, rminor = 0.4, = 70o . 8. Tests reveal that the gure of merit F Eq.(77) is never less than unity. The algorithm to generate EC consists of reducing the problem to two dimensions, and then solving a 2D problem. These two parts of the algorithm are described in the following two subsections. 11.3.1 Algorithm: reduction to 2D Following the development in Sec. 11.2, we transform to the principal axes of E (y-space), de ning y, p and w by Eq.(78), Eq.(79) and Eq.(82). Then EN is de ned by EN = {y | yT Fy 1}, (86) with F = 2 + wwT . 24 (87) Max. close content: |p|/rmin=1.6 Max. close content: |p|/rmin=1.2 0.5 0.5 0 0 0.5 0.5 1 0.5 0 0.5 1 1 0.5 0 0.5 1 Max. close content: |p|/r min =0.6 Max. close content: |p|/r min =0.2 0.5 0.5 0 0 0.5 0.5 1 0.5 0 0.5 1 1 0.5 0 0.5 1 Figure 10: Conservative shrinking algorithm. The initial ellipse E (solid blue line); the grow point p (magenta); the maximum-volume shrunk ellipse, EV , (green line); the near-content shrunk ellipse, EN , (red line); and the conservative shrunk ellipse, EC , (magenta): rmajor = 1, rminor = 0.4, = 70o . We now transform to z-space in which E is the unit ball: z = y = UT x. In z-space the de nition of EN , Eq.(86), becomes EN = {z | zT (I + ww)z 1}, with w = 1 w. And the maximum volume shrunk ellipsoid is EV = {z | zT (I + ppT )z 1}, where p is the transform of p 25 (91) (90) (89) (88) z2 1 EC EN 1 EV p 1 z1 E 1 Figure 11: Sketch of E, EN and EV in the z1 - 2 plane (in which E is the unit z disc) and of EC which covers EN and EV . The boundaries of both EN and EV intersect the grow point p. p = = UT (p c), p (92) (see Eq.(72)). The scalars and are determined by the condition that the boundaries of EN and EV intersect p (in z-space). Note that EN and EV are completely determined by two vectors, w and p. If these are co-linear, then EN = EV = EC . Otherwise we perform a rotation z = QT z, w = QT w, p = QT p, (93) such that w is in the z1 -direction, and p is in the z1 - 2 plane. The appropriate z orthogonal matrix Q is obtained from the QR decomposition QR = [w p]. 26 (94) Figure 11 is a sketch of the intersection of the ellipsoids with the z1 - 2 z plane. Note that in the other directions the principal semi-axes are unit vectors aligned with the coordinate axes. The minimum volume ellipsoid EC covering EN and EV is also shown in the gure. Its intersection with the z1 - 2 plane is given by the ellipse z z EC = {( 1 , z2 ) | LT z1 z2 1}, (95) where L is a 2 2 Cholesky matrix which is determined in the next subsection. Thus the ellipsoid EC is given by T EC = { | L z 1}, z (96) where the n n Cholesky matrix L is L= L0 0I . (97) Transforming Eq.(96) back to the original x-space, we obtain EC = {x | LT (x c)|| 1}, C where the Cholesky matrix LC is obtained from T AC = LC LC = BC BT , C (98) (99) (100) BC = U QL. 11.3.2 Algorithm: solution in 2D The problem to be solved is the determination of the 2 2 Cholesky matrix L de ning the covering ellipse EC in the z1 - 2 plane (see Fig. 11). z The scalars and in Eq.(89) and Eq.(91) are determined by the inter section condition (of EN and EV with p) to be = (1 | |2 )/ 2 , p p1 (101) 27 = (1 p|2 )/| |4 , p and the intersection of EN with the z1 axis is at z1 = = | 1 | p 1 p2 2 . (102) (103) A transformation to -space is performed 1 2 =C z1 z2 = 1/ 0 01 z1 z2 , (104) consisting of a stretching in the z1 direction so that EN transforms to the unit disc. With this transformation EV becomes EV = { | T A 1}, with A= The SVD of A is with U= ac cb 2 (1 + p2 ) 1 p2 1 p 1 p2 p 1 + p2 2 2 (105) = . (106) A = U UT , cos sin sin cos , (107) (108) (109) (110) (111) = 1 tan 1 ( 2c/(a b)), 2 2 = a cos2 2c sin cos + b sin2 , 11 2 = a sin2 + 2c sin cos + b cos2 . 22 Now EC is the ellipse which covers EN (which is the unit disc) and EV which has principal semi-axes whose directions are given by the columns of U and whose lengths are 1 and 1 . (In general, one of 11 and 22 is 22 11 less than unity and one greater than unity.) Thus the covering ellipse EC has the same principal directions 28 EC = { | T U 2 UT 1}, and eigenvalues 11 = min(1, 11 ), 22 = min(1, 22 ). (112) (113) (114) Transforming back to z1 - 2 we obtain the required result that EC is given z by Eq.(95), where the 2 2 Cholesky matrix L is obtained as 2T T LL = CU U C. (115) 12 Orthogonal projection of E onto a given line L {x | x = x0 + sv}, (116) We consider a given line L, parameterized by s, de ned by where x0 is a given point and v is a given non-zero vector, see Fig. 12. Given any point x in the space, its orthogonal projection onto L corresponds to s= vT (x x0 ) . vT v (117) Now the given ellipsoid E is given by E = {x | x = c + L T y, ||y|| 1}. Thus the projection of points in E correspond to values of s s= vT (L T y + c x0 ) vT v = s0 + wT y, for y 1, vT (c x0 ) , vT v 29 (118) (119) where s0 (120) s s+ c s0 L E s v x0 Figure 12: The orthogonal projection of the ellipsoid E onto the line L is the interval [s , s+ ]. and L 1 v . (121) vT v Given the condition y 1, it is evident from Eq.(119) that the orthogonal projection of E onto L corresponds to the interval [s , s+ ] in s, with w s s0 |w|. (122) This is sketched in Fig. 12. This method is implemented in the routine ell line proj. 13 Orthogonal projection of E onto an a ne space m We consider an ellipsoid E (in n ), an a ne space A (in and the orthogonal projection P (E) of E onto A. The ellipsoid E is given by E {x | x = c + L T u, The a ne space A is given by 30 u 1}. ), for 1 m n, (123) A {x | x = d + Tt}, (124) where d is a given m-vector, T is a given n m orthogonal matrix, and t is a vector of m parameters. The orthogonal projection of a general point x onto A is P (x) = d + TTT (x d). (125) We thus obtain P (E) = {x | x = d + TTT (c d + L T u), u 1}. Now let the SVD of TT L T be TT L T = U[ 0]VT . Then TT L T u = U w, where w denotes the rst m elements of w VT u. Note that u 1 implies w 1. Thus we obtain P (E) = {x | x = d + T( + U w), c where c TT (c d). (131) The above equation for P (E) is for an m-dimensional ellipsoid in A. It can be put in standard form by de ning B by B T = U , and then L as so that we can write T P (E) = {x | x = d + T( + L w), c (126) (127) (128) (129) w 1}, (130) (132) (133) T LL = BBT , w 1}. (134) The Cholesky matrix L can be computed from the LQ decomposition U 1 = LQ. This method is implemented in the routine ell a pr. 31 (135) 14 Generate an ellipsoid which does not cover any speci ed points We are given a point c, a set of P points p(j) (j = 1 : P ) and a positive length rmax . The problem is to generate an ellipsoid E which 1. is centered at c 2. has principal semi-axes no larger than rmax 3. does not cover any of the P points 4. and is as large as possible (in an unde ned sense). The algorithm used to solve this problem is implemented in the routine ell pts uncover. It involves a user-speci ed parameter (0 < 1) which a ects the shape of the resulting ellipsoid, E. The algorithm has two phases. In the rst phase there are n stages which generate a succession of ellipsoids E1 , E2 , . . ., En . In the second phase, E is formed by shrinking En uniformly and minimally so that none of the points is covered. In the rst phase, a principle axis is determined on each stage. The ellipsoid Ek is determined on the kth stage, and it has the following properties: 1. Ek is centered at c 2. for 1 < k, the th principle axis of Ek is the same as that of E (previously determined on stage ) 3. the kth principle axis (of half-length rk rmax ) is determined on the kth stage 4. for k < n, the th principle axis of Ek is of half-length rk . Note that E1 is a ball of radius r1 . An orthonormal basis is developed with basis vectors e1 , e2 , ..., en . On the kth stage e ( k) is modi ed, but subsequently ek is not altered. At the end of the kth stage, the basis vectors are principal axes of Ek . The vectors y(j) (j = 1 : P ) store the coordinates of the points (relative to c) in the current basis. For the jth particle we de ne 32 n hj = i=1 (yi /ri )2 . (j) (136) At the end of the kth stage (since r = rk for k) we can decompose hj as 2 hj = fj + gj /rk , (137) where fj = k 1 i=1 (yi /ri )2 , (j) (j) (138) and gj = n i=1 (yi )2 . (139) The ellipsoid Ek does not cover the point j if hj is greater than unity. This condition (hj > 1) can be re-expressed as 2 gj /rk > (1 fj ). (140) Points with fj > 1 cannot be covered by Ek regardless of how large rk is. Such points are excluded. Points with 2 < fj 1, (141) are partially excluded. Such points cannot be covered by Ek shrunk by a factor of . The remaining points (i.e., with fj 2 ) are included. We de ne rk to be the minimum value over the included points of gj /(1 2 fj ), and denote by the index of a point that achieves this minimum. The j signi cance of rk is that if rk is set to rk , then the point is on the boundary j of Ek , but no points are in the interior of Ek . If rk is greater than rmax (or if there are no included points), then we get r = rmax for all k, and omit the remaining stages of the rst phase. Otherwise rk is set to rk , and the ( j) basis vectors e (k n) are re-de ned so that y = 0 for > k. In this way, in subsequent stages, the ellipsoid can expand in directions orthogonal to y(j) , with y(j) remaining on the boundary. At the end of the rst phase there are no included points, and En does not cover any excluded points. However En may cover one or more partially excluded points. Consequently, the nal result E is obtained by shrinking En uniformly, as little as possible so that none of the partially included points is covered. 33 (a) E1 (b) O H xh v y1 yh E2 y2 E2 E1 H Figure 13: (a) Sketch of the non-intersecting ellipsoids E1 and E2 and a separating hyperplane H (b) Corresponding sketch in the transformed space in which E1 is the unit ball. 15 Separating hyperplane of two ellipsoids Given two ellipsoids, E1 and E2 , the task is to determine if they intersect; and, if they do not intersect, to determine a separating hyperplane, H. This is performed by the routine ell pair separate. The ellipsoid E1 is given by E1 {x | LT (x c1 ) 1}, 1 or, equivalently, E1 {x | y 1, y LT (x c1 )}. 1 With the transformation y LT (x c1 ), 1 (143) E1 is transformed to the unit ball at the origin (denoted by E1 ) and E2 is transformed to E2 , see Fig. 13. The point y2 in E2 which is closest to the origin is determined by the closest point algorithm (see Sec. 8). If y2 is less than unity, then E1 and E2 intersect. For the case in which E1 and E2 do not intersect, we de ne y1 = v = y2 / y2 , 34 (144) (142) so that y1 and y2 are the pair of closest points in E1 and E2 . Note that v is a unit vector. Then we de ne 1 yh 2 (y1 + y2 ), (145) and a separating hyperplane is de ned by H {y | vT (y yh ) = 0}. (146) Inverting the transformation Eq.(143), we obtain the separating hyperplane in the original space: H {x | uT (x xh ) = 0}, where u= and xh = c1 + L T yh . 1 (149) L1 v , L1 v (147) (148) For a hyperplane H given by Eq.(147) (for some u and xh ), we de ne the quality q as follows. Let x1 be the point in E1 closest to H, and similarly let x2 be the point in E2 closest to H. The quality q is de ned as q uT (x2 x1 ) 1. |x2 x1 | (150) This is the distance between the supporting hyperplanes at x1 and x2 relative to their separation. The maximum possible distance between the supporting hyperplanes is achieved for q = 1, and this occurs when x1 and x2 are the mutually closest points in E1 and E2 . In the routine ell pair separate there is an option to iteratively improve the quality of the hyperplane. Initially x1 and x2 are set to the points in E1 and E2 corresponding to y1 and y2 (in E1 and E2 ). Then, successively, x1 is replaced by the closest point in E1 to x2 ; and then x2 is replaced by the closest point in E2 to x1 . The hyperplane is then taken as the perpendicular bisector of the line x1 x2 . It is not guaranteed that this hyperplane is separating, nor that the quality increases with the iterations, although it generally does. Consequently, H is taken as the hyperplane with the greatest quality encountered initially or during the iterations. 35 (Note that the separating hyperplane of quality q = 1 could alternatively be obtained by solving the quadratic programming problem, of determining the mutually closest points x1 and x2 in E1 and E2 . This has not been implemented.) 16 Pair covering query Given a pair of ellipsoids E1 and E2 (centered at c1 and c2 and with Cholesky factors L1 and L2 ), the problem is to determine whether E1 covers E2 . This query is answered by the routine ell pair cover query by the following algorithm. With the same transformation as used in Sec. 15, E1 is transformed to the unit ball, at the origin, and E2 to E2 (see Eq.(143) and Fig.13). The algorithm described in Sec. 9 is then used to determine the distance s from the origin to the furthest point in E2 . The ellipsoid E1 covers E2 if and only if s is less than or equal to unity. 17 Shrink ellipsoid so that it is covered by a concentric ellipsoid Given two concentric ellipsoids, E1 , and E2 , the task is to form the maximalvolume ellipsoid, E, which is covered by both E1 and E2 . This is performed by the routine ell pair shrink. Clearly E is concentric with E1 and E2 , and so, without loss of generality, we take the origin at the mutual center. Then, E1 is given by E1 {x | LT x 1}, 1 and similarly for E2 and E (in terms of L2 and L, respectively). The transformation y = LT x, 1 (151) (152) maps E1 to the unit ball, E1 , and it maps E2 to the ellipsoid E2 with Cholesky factor L2 = L 1 L2 , (153) 1 see Fig. 14. Note that E1 shares the principal axes of E2 . Hence E (i.e., the mapped covered ellipsoid E) has the same principal directions as E2 ; and 36 (a) E2 E E1 (b) E1 E E2 Figure 14: (a) Sketch of the concentric ellipsoids E1 and E2 and the maximalvolume ellipsoid E which is covered by them (b) Corresponding ellipsoids in the transformed space in which E1 is the unit ball. the lengths of its principal axes are the lesser of those of E1 (all of which are unity) and those of E2 . Thus if L2 = U VT , is the SVD of L2 , so that L2 L2T = U 2 UT , then the Cholesky matrix of E is given by L L T = U UT , where the singular values ii = i are given by i = max( i , 1), (157) 2 (154) (155) (156) where i ii . The Cholesky matrix of E is obtained by inverting the transformation: L = L1 L . (158) The same algorithm can be used to determine the minimal-volume ellipsoid which covers E1 and E2 . In that case, in contrast to Eq.(157), the appropriate singular values are: i = min( i , 1). 37 (159) 18 Ellipsoid that covers two given ellipsoids Given two ellipsoids E1 and E2 , the task is to determine a third ellipsoid E that covers both E1 and E2 . Ideally E is of minimal volume. It is a problem of convex optimization to determine the minimum-volume covering ellipsoid (see, e.g., Boyd and Vandenberghe 2004). An algorithm is provided by Yildirim (2006). It appears that the solution to this convex optimization problem is computationally expensive. Instead, in the subsections, we describe heuristic algorithms with determine ellipsoids E (not of minimal volume) which cover E1 and E2 . These methods are implemented in the routine ell pair cover, with the parameter algorithm determining the particular algorithm to be used. 18.1 Spheroid algorithm The ellipsoids E1 and E2 have centers c1 and c2 , and outer radii rout,1 and rout,2 . Thus the ball B1 centered at c1 of radius rout,1 covers E1 , and similarly we de ne the ball B2 which covers E2 . In this spheroid algorithm (algorithm = 1) we take E to be the minimum-volume ball B which covers B1 and B2 . The center and radius of B are determined as follows. Consider the line of centers with distance s measured from c1 towards c2 . Let the distance between the centers be c | c2 c1 |. There are four intersections between the ball B1 and B2 and the line. The outermost of these corresponds to smax = max(rout,1 , c + rout,2 ), and smin = min( rout,1 , c rout,2 ). Thus the center of B is 1 c = c1 + (smin + smax )(c2 c1 )/ c, 2 and its radius is 1 r = (smax smin ). 2 (162) (161) (160) (163) 38 A variant is the spheroid algorithm with shrinking (algorithm = 4), in which r is decreased to r to yield the ball B of minimum volume centered at c which covers E1 and E2 . The radius r is determined as the greatest distance from c to any point in E1 and E2 (which is determined by ell pt near far). 18.2 Covariance algorithm E1 {x | (x c1 )T A1 (x c1 ) 1}, (164) The ellipsoid E1 is de ned by and E2 and E are similarly de ned by c2 and A2 , and by c and A, respectively. In the covariance algorithm described in this section, we de ne the covering ellipsoid E by 1 (165) c 2 (c1 + c2 ), A A0 , and A0 A 1 + A 1 + 1 [c1 c2 ][c1 c2 ]T 1 2 4 1 (166) , (167) where is a positive parameter to be determined. To determine , we consider the ellipsoid E0 de ned by c and A0 . With L0 being the Cholesky factor L0 LT = A0 , 0 we perform the linear transformation y = c + L T x. 0 (169) (168) As depicted in Fig. 15, this transforms E0 to the unit ball, and E1 and E2 to ellipsoids denoted by E1 and E2 . Using the furthest-point algorithm (see Sec. 9), we determine y1 and y2 , de ned as the points on E1 and E2 , respectively, which are furthest from the origin (see Fig. 15). Clearly, the ball E of radius r max( y1 , y2 ) (170) covers E1 and E2 . This corresponds to E (the transformation of E) with 39 E E1 E0 y1 E2 y2 Figure 15: In the transformed space, the ellipsoid E0 which is the unit ball, the ellipsoids E1 and E2 , and the covering ellipsoid E . = r 2 . (171) In summary the ellipsoid E which covers E1 and E2 is de ned by c Eq.(165) and A Eq.(166), where A0 and are given by Eqs. (167), (170), and (171). The routine ell pair cover includes two implementations of this covariance algorithm. The rst (algorithm = 5) is a direct implementation of the above equations in which A0 is formed and then L0 is obtained by Cholesky decomposition. The second implementation (algorithm = 2), now described, is preferred, since it avoids the formation of A0 , and hence is signi cantly more accurate. With L1 and L2 being the Cholesky factors of A1 and A2 , and with d 1 (c1 c2 ), we de ne the (2n + 1) n matrix B by 2 BT = [L T 1 L T 2 d], (172) so that the inverse of Eq.(167) can be written A 1 = L T L 1 = BT B. 0 0 0 Now let the QL factorization of B be B=Q so that 40 0 L , (174) (173) BT B = LT L. (175) By comparing the above two equations, we see that L0 is the inverse of L. Thus, in this preferred QL implementation of the covariance algorithm, the Cholesky matrix L0 required in Eq.(169) is obtained as L 1 , where L is obtained from the QL factorization of B, de ned by Eq.(172). 18.3 Iterative algorithm The algorithm described here (which is also implemented in the routine ell pair cover with algorithm= 2), is more elaborate and expensive than the covariance algorithm described in the previous subsection, but in most circumstances it generates a covering ellipsoid E of smaller volume. Given E1 and E2 , the algorithm proceeds through six stages (described below) to generate the covering ellipsoid E. Stages 1, 2 and 6 are trivial, but are retained for consistency with the implementation in ell pair cover. In Stages 3 and 4 the shape of E (but not its center and size) are determined. The center of E is taken to be on the line of centers of E1 and E2 . Its location is determined iteratively in Stage 5 so as to minimize the volume of E. Various spaces are considered, and are referred to as x-space, y-space, z-space and -space. The ellipsoids E1 and E2 are given in x-space, and the covering ellipsoid E is to be determined in this space. The other spaces are obtained by successive linear transformations; and E2 (z), for example, denotes the ellipsoid E2 viewed in z-space, which has center c2 (z) and Cholesky triangle L2 (z). 18.3.1 Stage 1 The two ellipsoids E1 and E2 are given in x-space (in terms of c1 , L1 , c2 and L2 ): see Fig.16(a). 18.3.2 Stage 2 0 0 We denote by E1 and E2 the two given ellipsoids shifted to the origin: see 0 Fig.16(b). (In general, Em denotes Em shifted to the origin.) 18.3.3 Stage 3 The transformed variable y is de ned by 41 (a) E1 c1 (b) O E2 c2 0 E1 0 E2 Figure 16: In x-space, sketches of (a) the ellipsoids E1 and E2 in stage 1 (b) the 0 0 ellipsoids E1 and E2 in stage 2. E3 E1 O 0 E2 E2 c2 (y) Figure 17: In y-space, sketches of the various ellipsoids in stage 3. y L 1 (x c1 ), 1 (176) so that E1 (y) is the unit ball at the origin. At this stage a test is performed to determine if E1 covers E2 ; if it does, then E = E1 is the minimal volume covering ellipsoid. Similarly, a test is performed to determine if E2 covers E1 . The way in which this testing is performed is described in Section 18.3.7. If neither E1 nor E2 covers the other, then the algorithm proceeds. The ellipsoid E3 is de ned to be the minimal-volume ellipsoid which cov0 0 ers E1 and E2 (see Fig.17). This is readily determined: E3 (y) has the same 0 principal directions as E2 , and the lengths of its principal semi-axes are the 0 0 greater of those of E2 (y) and E1 (y) (which are unity). 18.3.4 Stage 4 A transformation (to z-space) is performed which consists of a rotation such that the principal axes of E1 , E2 and E3 are aligned with the coordinate 42 smax E2 B2 w E1 smin E3 = B1 L E4 Figure 18: In z-space, sketch of ellipsoids in stage 4 showing the construction of E4 based on the extreme intersections (smin and smax ) of the bounding balls B1 and B2 with the line of centers L. axes, followed by a stretching in the coordinate directions to make E3 (z) the unit ball at the origin (see Fig.18). This transformation is readily determined from the SVD of L2 (y). Note that the principal semi-axes of E1 (z) and E2 (z) are aligned with the coordinate directions; their lengths are at most unity; and in each direction the length of the principal semi-axis of E1 (z) and/or E2 (z) equals unity. The line of centers L between E1 (z) and E2 (z) can be written L {z | z = ws}, where w = c2 (z)/|c2 (z)| is the unit vector from the origin to the center of E2 (z), and s is the distance along the line. Let B1 denote the bounding ball of E1 (z). This is centered at the origin and has radius less than or equal to unity. There are two intersections of B1 and L which occur at distances s1 and s1+ along L. Similarly the bounding ball B2 of E2 (z) intersects L at distances s2 and s2+ . We de ne the extrema 43 of these intersections by smin min(s1 , s1+ , s2 , s2+ ), and smax max(s1 , s1+ , s2 , s2+ ). The ellipsoid E4 is now de ned to be centered at s0 1 (smin + smax ), 2 (179) (178) (177) to have a principal semi-axis of length 1 (smax smin ) in the direction w, and 2 to have all other principal axes unity. In other words, E4 (z) is formed from the unit ball at s0 by stretching it in the w directions so that it intersects the extrema smin and smax . 18.3.5 Stage 5 A transformation is performed to -space such that E4 ( ) is the unit ball at the origin, and the line of centers L( ) is in the rst coordinate direction: L( ) = { | i = t i1 }, (180) where t measures the distance along the line (see Fig.19). The covering ellipsoid being constructed E( ) is a ball of radius r0 centered at a distance t0 along L. It remains to determine r0 and t0 . Consider a ball centered at a distance t along L. Let r1 (t) denote the distance from the ball s center to the furthest point in E1 ( ). This can be determined by the algorithm described in Section 9. Similarly, let r2 (t) denote the distance to the furthest point in E2 ( ); and we de ne r(t) max(r1 (t), r2 (t)). (181) Thus the ball centered at a distance t along L and of radius r(t) covers both E1 ( ), and E2 ( ). In the de nition of E( ), we take t0 to be (an approximation to) the value of t at which r(t) is minimum, and then de ne r0 r(t0 ) As one moves along the line L from the center of E1 ( ) to the center of E2 ( ), r1 (t) continually increases and r2 (t) continually decreases. Except in unusual circumstances, the minimum of r(t) occurs between the centers, where r1 (t) equals r2 (t). A simple iterative procedure usually determines the location of the minimum (to reasonable accuracy) in two or three iterations. 44 r1 (t) t r2 (t) E 2 O L E1 E4 Figure 19: In -space, sketch showing ellipsoids in stage 5, the general point t on the line of centers L, and the distances r1 (t) and r1 (t) to the furthest points in E1 and E2 , respectively. 18.3.6 Stage 6 The covering ellipsoid E( ) obtained in Stage 5 is transformed to the original coordinate system to yield the required ellipsoid E(x) which covers E1 and E2 . 18.3.7 Mutual covering In Stage 3 it is required to determine if E1 (y) covers E2 (y) and vice versa. The procedure used is now described. Let be the distance from the origin to the furthest point in E2 (y). Since E1 (y) is the unit ball at the origin, E1 covers E2 if, and only if, is less than or equal to unity. While can be evaluated using the algorithm described in Section 9, this involves quadratic minimization and hence is somewhat expensive. As now described, there are easily computed bounds on so that its exact evaluation can often be avoided. Let d be the distance between the centers of E1 (y) and E2 (y), and let r2,min and r2,max be the lengths of the smallest and largest principal semi-axes of E2 (y). From the triangle inequality we then have 45 d + r2,max , and d + r2,min . We also have r2,max . (182) (183) (184) Thus, if d + r2,max 1, then 1, and so E1 covers E2 . On the other hand, if d + r2,min > 1 or r2,max > 1, then > 1, and so E1 does not cover E2 . In the remaining cases, E1 may cover E2 , and so is evaluated. 18.3.8 Discussion While this algorithm is more elaborate and expensive than that given in Section 18.2, it generally yields a covering ellipsoid E of smaller volume. In particular it yields the ellipsoid of minimal volume if E1 and E2 are concentric or if one covers the other. The principal computational expenses are: 1. The SVD of L2 (y) performed in Stage 3 2. The quadratic minimization involved in the furthest-point algorithm (Section 9) which is invoked 0, 1 or 2 times (in Stage 3) to determine if E1 and E2 cover each other 3. The quadratic minimization involved in the furthest-point algorithm which is invoked twice per iteration (to evaluate r1 (t) and r2 (t)) in Stage 5. 19 Conclusions Algorithms have been described for performing some basic geometric operations on ellipsoids. A Fortran implementation of these algorithms is provided by the Ell LIB library. 20 Acknowledgments I am particularly grateful to Professor Charles Van Loan for his help on numerous occasion with issues of numerical linear algebra. 46 References Averick, B. M., R. Carter, and J. More (1993). MINPACK-2. http://wwwfp.mcs.anl.gov/OTC/minpack/summary.html. Boyd, S. and L. Vandenberghe (2004). Convex Optimization. Cambridge: Cambridge University Press. Golub, G. H. and C. F. Van Loan (1996). Matrix Computations (3rd ed.). Baltimore: Johns Hopkins University Press. Pope, S. B. (1997). Computationally e cient implementation of combustion chemistry using in situ adaptive tabulation. Combust. Theory Modelling 1, 41 63. Yildirim, E. A. (2006). On the minimum volume covering ellipsoid of ellipsoids. SIAM Journal on Optimization 17, 621 641. 47

Find millions of documents here - Study Guides, Homework Solutions, Papers, Exam Answer Keys and more. Course Hero has millions of course related materials that will enable you to learn better, faster and get an A in all your courses.
Below is a small sample set of documents: