Computing Hulls And Centerpoints In Positive Definite Space

Computing Hulls And Centerp oints In P ositive Definite Space ∗ P . Thomas Fletcher fletcher@sci.utah.edu John Moeller moeller@cs.utah.edu Jef f M. Phillips jeffp@cs.utah.edu Suresh V enkatasubramanian suresh@cs.utah.edu Abstract In this paper , we present algorithms for computing approximate hulls and centerpoints for collections of matrices in positiv e definite space. There are many applications where the data under consideration, rather than being points in a Euclidean space, are positiv e definite (p.d.) matrices. These applications include dif fusion tensor imaging in the brain, elasticity analysis in mechanical engineering, and the theory of kernel maps in machine learning. Our work centers around the notion of a horoball : the limit of a ball fixed at one point whose radius goes to infinity . Horoballs possess many (though not all) of the properties of halfspaces; in particular , they lack a strong separation theorem where two horoballs can completely partition the space. In spite of this, we show that we can compute an approximate “horoball hull” that strictly contains the actual con ve x hull. This approximate hull also preserves geodesic extents, which is a result of independent value: an immediate corollary is that we can approximately solve problems like the diameter and width in positiv e definite space. W e also use horoballs to show existence of and compute approximate robust centerpoints in positiv e definite space, via the horoball- equiv alent of the notion of depth. ∗ This research was supported in part by NSF SGER-0841185 and a subaward to the Univ ersity of Utah under NSF award 0937060 to the Computing Research Association. 1 Intro duction There are many application areas where the basic objects of interest, rather than points in Euclidean space, are symmetric positive-definite n × n matrices (denoted by P ( n ) ). In diffusion tensor imaging [3], matrices in P (3) model the flo w of water at each v oxel of a brain scan. In mechanical engineering [11], stress tensors are modeled as elements of P (6) . Kernel matrices in machine learning are elements of P ( n ) [25]. In all these areas, a problem of great interest is the analysis [13, 14] of collections of such matrices (finding central points, clustering, doing regression). For all of these problems, we need the same kinds of geometric tools a v ailable to us in Euclidean space, including basic structures like halfspac es, con vex hulls, V oronoi diagrams, v arious notions of centers, and the like. P ( n ) is non-Euclidean; in particular , it is negati vely (and variably) curved, which poses fundamental problems for the design of geometric algorithms. This is in contrast to hyperbolic space (which has constant curv ature of − 1 ), in which many standard geometric algorithms carry ov er . In this paper , we develop a number of basic tools for manipulating positiv e definite space, with a focus on appli- cations in data analysis. 1.1 Our W ork Ho roballs. A main technical contrib ution of this work is the use of hor oballs as generalization of halfspaces. Suppose we allo w a ball to gro w to infinite radius while al ways touching a fix ed point on its boundary . In Euclidean space, this construction yields a halfspace; in general Cartan-Hadamard manifolds (of which P ( n ) is a special case), this construction yields a horoball. Because of the curvature of space, horoballs are not flat and the complement of a horoball is not a horoball. Howe ver , we sho w that these objects can be effecti vely used as proxies for halfspaces, allo wing us to define a number of different geometric structures in P ( n ) . Ball Hulls. The first structure we study is the conv ex hull. Apart from its importance as a fundamental primitive in computational geometry , the con vex hull also provides a compact description of the boundary of a data set, can be used to define the center of a data set (via the notion of con vex hull peeling depth [23, 2]), and also captures e xtremal properties of a data set like its diameter , width and bounding volume (e ven in its approximate form [1]). The con vex hull of a set of points in P ( n ) can be naturally defined as the intersection of all conv ex sets containing the points. Alternativ ely , it can be defined as the set of all points that are “con vex combinations” (in a geodesic sense) of the input points. A significant obstacle to the con vex hull in P ( n ) is that it is not even kno wn whether the con ve x hull of a finite collection of points in P ( n ) can be represented finitely [4]. Another approach to defining the conv ex hull is via halfspaces: we can define the conv ex hull in Euclidean space as the intersection of all halfspaces that contain all the points. Unfortunately , ev en this notion f ails to generalize: the rele vant structures are called totally geodesic submanifolds , and we cannot guarantee that any set of d + 1 points admits such a submanifold passing through them. Our main technical contribution here is a generalization of the con ve x hull called the ball hull that is based on the relationship between horoballs and halfplanes. The ball hull is the intersection of all horoballs that contain the input points. Although the ball hull itself might require an infinite number of balls to define it, it is closed, it can be approximated ef ficiently , it is identical to the con vex hull in Euclidean space, and it al ways contains the con ve x hull in P ( n ) . In the process of proving this result, we also develop a generalized notion of extent [1] in positive definite space that might be of independent interest for other analysis problems. Centerp oints. One important motiv ation for studying collections of points in positiv e definite space is to compute measures of centrality (or mean shapes ) [14]. A robust centerpoint can be obtained by finding a point of maximum (halfspace) depth among a collection of points. W e first prov e, using a generalization of Helly’ s theorem to negati vely curved spaces, that for any set of points in P ( n ) , there e xists a point of lar ge depth, where depth is defined in terms of horoballs. W e then dev elop an algorithm to compute an approximation to such a point, using an LP-type framework. 1 The point we compute is a geometric approximation: it does not approximate the depth of the optimal point, but is guaranteed to be close to such a point. 1.2 Related W ork The mathematics of Riemannian manifolds, Cartan-Hadamard manifolds and P ( n ) is well-understood: the book by Bridson and Haefliger [6] is an in valuable reference on metric spaces of nonpositiv e curvature, and Bhatia [5] provides a detailed study of P ( n ) in particular . Ho wev er , there are many fewer algorithmic results for problems in these spaces. T o the best of our kno wledge, the only prior work on algorithms for positive definite space are the work by Moakher [21] on mean shapes in positiv e definite space, and papers by Fletcher and Joshi [13] on doing principal geodesic analysis in symmetric spaces, and the robust median algorithms of Fletcher et al [14] for general manifolds (including P ( n ) and SO ( n ) ). Geometric algorithms in hyperbolic space are much more tractable. The Poincar ´ e and Klein models of hyperbolic space preserve different properties of Euclidean space, and many algorithm carry over directly with no modifica- tions. Leibon and Letscher [18] were the first to study basic geometric primitives in general Riemannian manifolds, constructing V oronoi diagrams and Delaunay triangulations for sufficiently dense point sets in these spaces. Epp- stein [12] described hierarchical clustering algorithms in hyperbolic space. Krauthgamer and Lee [16] studied the nearest neighbor problem for points in δ -hyperbolic space; these spaces are a combinatorial generalization of neg- ati vely curved space and are characterized by global, rather than local, definitions of curv ature. Chepoi et al [8, 9] adv anced this line of research, providing algorithms for computing the diameter and minimum enclosing ball of collections of points in δ -hyperbolic space. 2 Prelimina ries P ( n ) is the set of symmetric positiv e-definite real matrices. It is a Riemannian metric space with tangent space at point p equal to S ( n ) , the vector space of symmetric matrices with inner product h A, B i p = tr( p − 1 Ap − 1 B ) . The exp map, exp p : S ( n ) → P ( n ) is defined exp p ( tA ) = c ( t ) = pe tpA , where c ( t ) is the geodesic with unit tangent A and c (0) = p . For simplicity , we often assume that p = I so exp I ( tA ) = e tA . The log map, log p : P ( n ) → S ( n ) , indicates direction and distance and is the in verse of exp p . The metric d ( p, q ) = k log p ( q ) k = p tr(log( p − 1 q ) 2 ) . Convex Hulls in P ( n ) . P ( n ) is an e xample of a proper CA T (0) space [6, II.10], and as such admits a well-defined notion of con vexity , in which metric balls are con ve x. W e can define the conv ex hull C ( X ) of a set of points X as the smallest con ve x set that contains the points. This hull can be realized as the limit of an iterativ e procedure where we draw all geodesics between data points, add all the ne w points to the set, and repeat. Lemma 2.1 ([5]) . If X 0 = X and X i +1 = S a,b ∈ X i [ a, b ] , then C ( X ) = S ∞ i =0 X i . Pr oof. W e will use the notation X ∞ = S ∞ i =0 X i . It is easy to demonstrate by straightforward induction that X ∞ is contained in any con vex set that contains X . Therefore C ( X ) ⊇ X ∞ . W e also kno w that if p, q ∈ X ∞ there must be some m for which p, q ∈ X m , since X ∞ is the nested union of X i . Then [ p, q ] ⊂ X m +1 ⊂ X ∞ . This means that X ∞ is con ve x, so C ( X ) ⊆ X ∞ . Berger [4] notes that it is unknown whether the conv ex hull of three points is in general closed, and the standing conjecture is that it is not. The abov e lemma bears this out, as it is an infinite union of closed sets, which in general is not closed. These facts present a significant barrier to the computation of con ve x hulls on general manifolds. 2.1 Busemann F unctions In R d , the con vex hull of a finite set can be described by a finite number of hyperplanes each supported by d points from the set. A hyperplane through a point may also be thought of as the limiting case of a sphere whose center has 2 been mov ed away to infinity while a point at its surface remains fixed. W e generalize this notion with the definition of a Busemann function. For this notion to work, we must restrict ourselves to a class of spaces called CA T (0) spaces. The y are metric spaces with non-positi ve curvature. Additionally , they must be complete ; that is, Cauchy sequences in the space must conv erge to a point in the space. Euclidean space, hyperbolic space, and P ( n ) are all examples of complete CA T (0) spaces. T o talk about “sending a point away to infinity , ” we must provide a rigorous definition of what we mean by infinity in a complete CA T (0) space. T wo geodesic rays c 1 , c 2 : R + → M in a complete CA T (0) space M are asymptotic if lim t →∞ d ( c 1 ( t ) , c 2 ( t )) < K for some K ∈ R + . If c 1 and c 2 are asymptotic, then we say c 1 ∼ c 2 . This forms an equiv alence relation ∼ so that [ c ] describes the set of all geodesics c 0 such that c ∼ c 0 . Let ξ = [ c ] where ξ is identified with the limit of any geodesic ray asymptotic to c . W e say that ξ is a point at infinity . Moreov er , for any point x ∈ M we can find a member of [ c ] that issues from x [6, II.8]. Definition 2.1. F or a complete CAT (0) space M , given a geodesic ray c ( t ) : R + → M , a Busemann function b c : M → R is defined b c ( p ) = lim t →∞ d ( p, c ( t )) − t. It should be noted that if we construct a Busemann function from an y geodesic ray in [ c ] , it is the same function up to addition by a constant [6, II.8]. It’ s con venient then to normalize a Busemann function by assuring that b c ( I ) = 0 . A Busemann function is an e xample of a hor ofunction [6, II.8]. A hor ospher e S r ( h ) ⊂ M is a lev el set of a horofunction h ; that is, S r ( h ) = h − 1 ( r ) , where r ∈ R . A hor oball B r ( h ) ⊂ M is a sublev el set of h ; that is, B r ( h ) = h − 1 (( −∞ , r ]) . Horofunctions are con ve x [6, II.8], so any sublev el set of a horofunction is conv ex, and therefore any horoball is con vex. Example: Busemann functions in R n . As an illustration, we can easily compute the Busemann function in Euclidean space associated with a ray c ( t ) = t u , where u is a unit vector . Since lim t →∞ 1 2 t ( k p − t u k + t ) = 1 , b c ( p ) = lim t →∞ ( k p − t u k − t ) = lim t →∞ 1 2 t ( k p − t u k 2 − t 2 ) = lim t →∞ 1 2 t ( k p k 2 − 2 h p, t u i + k t u k 2 − t 2 ) = lim t →∞ k p k 2 2 t − h p, u i = − h p, u i . Horospheres in Euclidean space are then just hyperplanes, and horoballs are halfspaces. 2.1.1 Decomp osing P ( n ) In order to construct Busemann functions in P ( n ) it is necessary to decompose the space into simpler components. The notion of a hor ospherical pr ojection will be very useful. The horospherical group. There is a subgroup of GL ( n ) , N ξ (the hor ospherical gr oup ), that leav es the Buse- mann function b c in v ariant [6, II.10]. That is, gi ven p ∈ P ( n ) , and ν ∈ N ξ , b c ( ν pν T ) = b c ( p ) . Let A be diagonal, where A ii > A j j , ∀ i > j . Let c ( t ) = e tA , and ξ = c ( ∞ ) . Then ν ∈ N ξ if and only if ν is a upper-triangular matrix with ones on the diagonal 1 . If A ∈ S ( n ) is not sorted-diagonal, we may still use this characterization of N ξ without loss of generality , since we may compute an appropriate diagonalization A = QA 0 Q T , QQ T = I , then apply the isometry Q T pQ to any element p ∈ P ( n ) . 1 For simplicity , we consider only those rays with unique diagonal entries, but this definition may be e xtended to those with multiplicity . 3 horospheres Figure 1: Left: projection of X ⊂ P (2) onto det( x ) = 1 . Right: X ⊂ P (2) . Two ho rospheres are drawn in b oth views. Flats. Let A ∈ S ( n ) and c ( t ) = e tA as above. If we consider all elements f ∈ P ( n ) that share eigen vectors Q with e A , then all such elements commute with each other and f e A = e A f . W e call this space F , the n -flat containing c . Since we may assume that Q ∈ SO ( n ) , ev ery flat F corresponds to an element of SO ( n ) . Moreov er , since members of F commute, log( ab ) = log a + log b for all a, b ∈ F . So if u and v are in F , then the distance between them is p tr(log( u − 1 v ) 2 ) = p tr((log( v ) − log( u )) 2 ) . Since p tr(( · ) 2 ) is a Euclidean norm on log ( F ) , we have that F is isometric to R n with a Euclidean metric under log( · ) . Ho rospherical p rojection. Gi ven p ∈ P ( n ) , there is a unique decomposition p = ν f ν T where ( ν , f ) ∈ N ξ × F [6, II.10]. Let p ∈ P ( n ) and ( ν, f ) ∈ N ξ × F . If p = ν f ν T , then define the hor ospherical pr ojection function π F : P ( n ) → F as π F ( p ) = ν − 1 pν − T = f . 2.1.2 Busemann functions in P ( n ) . W e can no w giv e an explicit expression for a Busemann function in P ( n ) . For geodesic c ( t ) = e tA , where A ∈ S ( n ) , the Busemann function b c : P ( n ) → R is b c ( p ) = − tr( A log ( π F ( p ))) , where π F is defined as abov e [6, II.10]. In P (2) it is conv enient to visualize Busemann functions through horospheres. W e can embed P (2) in R ¯ 3 where the log of the determinant of elements grows along one axis. The orthogonal planes contain a model of hyperbolic space called the P oincar ´ e disk that is modeled as a unit disk, with boundary at infinity represented by the unit circle. Thus the entire space can be seen as a cylinder , as sho wn in Figure 1. W ithin each cross section with constant determinant, the horoballs are disks tangent to the boundary at infinity . 3 Ball Hulls W e now introduce our variant of the con vex hull in P ( n ) , which we call the ball hull. For a subset X ⊂ P ( n ) , the ball hull B ( X ) is the intersection of all horoballs that also contain X : B ( X ) = \ b c ,r B r ( b c ) , X ⊂ B r ( b c ) . Note that the ball hull can be seen as an alternate generalization of the Euclidean con ve x hull (i.e. via intersection of halfspaces) to P ( n ) . Furthermore, since it is the intersection of closed sets, it is itself guaranteed to be closed. 3.1 Prop erties Of The Ball Hull W e know that an y horoball is con ve x. Because the ball hull is the intersection of conv ex sets, it is itself con ve x (and therefore C ( X ) ⊆ B ( X ) ). W e can also show that it shares critical parts of its boundary with the conv ex hull (Theorem 3.1), but unfortunately , we cannot represent it as a finite intersection of horoballs (Theorem 3.2). 4 Theorem 3.1. Every x ∈ X ( X finite) on the boundary of B ( X ) is also on the boundary of C ( X ) (i.e., X ∩ ∂ B ( X ) ⊆ X ∩ ∂ C ( X ) ). Pr oof. Since X ⊂ C ( X ) , either x ∈ ∂ C ( X ) or x ∈ int C ( X ) . Assume that x ∈ int C ( X ) . Then there is a neighborhood U of x contained wholly in C ( X ) . Because x ∈ ∂ B ( X ) , there is a horofunction h such that X ⊂ B r ( h ) and h ( x ) = r . Since B r ( h ) is conv ex, U ⊂ C ( X ) ⊆ B r ( h ) . This implies that h ( x ) < r , but h ( x ) = r , a contradiction. Thus x ∈ ∂ C ( X ) . Theorem 3.2. In general, the ball hull cannot be described as the intersection of a finite set of hor oballs. Pr oof. W e construct an example in P (4) with a point set X = { x 1 , x 2 } of size 2 where the ball hull cannot be described as the intersection of a finite number of horoballs. Let X ⊂ H 3 , the three dimensional hyperbolic space, as embedded in P (4) . In particular , let the geodesic that contains x 1 and x 2 also contain I , the identity matrix, at their midpoint on the geodesic. Consider the family of horofunctions H such that for h ∈ H , h ( x 1 ) = h ( x 2 ) . By construction, h ( x 1 ) = h ( x 2 ) = r for some r ∈ R . In the Poincar ´ e ball model of H 3 , the horospheres S r ( h ) are literally spheres that are tangent to the boundary at infinity , and touch x 1 and x 2 . So constructed, any horosphere in H will contact the boundary at infinity on a great circle that is equidistant from both points. The ball hull B ( X ) is defined { T B r ( h ) | h ∈ H } , and is a “spindle” (a three dimensional lune) with tips at x 1 and x 2 and bulges out from the geodesic segment between them. Any finite family of horofunctions will intersect in a region strictly larger than B ( X ) , so ev ery horoball generated by a member in H is necessary for B ( X ) , and so there is no finite set of horoballs that describe B ( X ) . 4 The ε -Ball Hull Theorem 3.2 indicates that we cannot maintain a finite representation of a ball hull. Howe ver , as we shall show in this section, we can maintain a finite-sized appr oximation to the ball hull. Our approximation will be in terms of extents : intuitively , we say that a set of horoballs approximates the ball hull if a geodesic traveling in any direction trav erses approximately the same distance inside the ball hull as it does inside the approximate hull. ho ro extent In Euclidean space, we can capture extent by measuring the distance between two parallel hy- perplanes that sandwich the set. Since measuring extent by recording the distance between two parallel planes does not ha ve a direct analogue in P ( n ) , we define a notion we call a hor oextent . Let c ( t ) = q e tq − 1 A be a geodesic, and X ⊂ P ( n ) . The hor oe xtent E c ( X ) with respect to c is defined as: E c ( X ) =     max p ∈ X b c + ( p ) + max p ∈ X b c − ( p )     , where b c + is the Busemann function created when we allow t to approach positi ve infinity as normal, while b c − is the Busemann function created when we allo w the limit to go the other direction; that is: b c + ( p ) = lim t → + ∞ ( d ( c ( t ) , p ) − t ) , b c − ( p ) = lim t →−∞ ( d ( c ( t ) , p ) + t ) . Observe that for any c , E c ( X ) = E c ( C ( X )) = E c ( B ( X )) . (Note that we cannot simply substitute min b c + for max b c − ; the horoballs are generated by Busemann functions tied to opposite points at infinity .) If we were to compute E c ( X ) in a Euclidean space, it would be apparent that the extent would be the distance between two parallel planes. Because the minimum distance between two Euclidean horoballs is a constant, no matter where we place the base point q , we can measure horoextent simply by measuring width in a particular direction. Ho we ver , this is not true in general. For instance in P ( n ) , horofunctions are nonlinear , so the distance between opposing horoballs is not constant. The width of the intersection of the opposing horoballs is taken along the geodesic c , and a geodesic is described by a point q and a direction A . W e fix the point q so that we need only choose a uniform grid of directions A for our approximation. 5 Definition 4.1. An intersection of hor oballs is called an ε -ball hull with origin q ( B ε,q ( X ) ) if for all g eodesic rays c such that c (0) = q , | E c ( B ε,q ( X )) − E c ( X ) | ≤ ε . When q is clear , we will refer to B ε,q ( X ) as just an ε -ball hull. Shifting the o rigin. Let the geodesic anisotr opy [21] of a point p ∈ P ( n ) be defined as GA( p ) = d ( n p det( p ) I , p ) (so if det( p ) = 1 then GA( p ) = d ( I , p ) ). Let d X = max p ∈ X d ( p, I ) ≥ max p ∈ X GA( p ) . The size of the ε -ball hulls we construct will depend d X , but this is not an intrinsic parameter of the data, since we can change it merely by isometrically translating the point set. For some point q ∈ C ( X ) , if we could translate the data set so that q was at the origin I , then d X ≤ diam( X ) = max p,q ∈ X d ( p, q ) . W e no w prov e that such a translation is always possible. Lemma 4.1. F or a point q ∈ P ( n ) , a geodesic c such that c (0) = q , and a point set X , if we define an operation ˆ S on a set S ∈ P ( n ) such that ˆ p = q − 1 2 pq − 1 2 for any p ∈ S , then E c ( X ) = E ˆ c ( ˆ X ) . Pr oof. Let c ( t ) = q e tq − 1 A . Then ˆ c ( t ) = q − 1 2 ( q e tq − 1 A ) q − 1 2 = q 1 2 e tq − 1 A q − 1 2 = e tq − 1 2 Aq − 1 2 = e t ˆ A . Note that k A k q = q tr(( q − 1 2 Aq − 1 2 ) 2 ) = k ˆ A k I , so ˆ c is a geodesic such that ˆ c (0) = I with the same speed as c . Let b c be the Busemann function of c , so b c ( p ) = lim t →∞ ( d ( c ( t ) , p ) − t ) . Since conjugation by q − 1 2 is an isometry on P ( n ) , b c ( p ) = lim t →∞ ( d ( c ( t ) , p ) − t ) = lim t →∞ ( d (ˆ c ( t ) , ˆ p ) − t ) = b ˆ c ( ˆ p ) , and therefore E c ( X ) = E ˆ c ( ˆ X ) . For con venience, we will now assume that our data has been shifted into a reasonable frame where some point q ∈ C ( X ) is the base point of our horofunction. That is in this shifted frame I ∈ C ( X ) and, we can bound, d X ≤ diam( X ) . Main result. The main result of this section is a construction of a finite-sized ε -ball hull. Theorem 4.1. F or a set X ⊂ P ( n ) of size N (for constant n ), we can construct an ε -ball hull of size O ((sinh( d X ) /ε ) n − 1 · N b n/ 2 c ) in time O ((sinh( d X ) /ε ) n − 1 ( N b n/ 2 c + N log N )) . Pro of Overview. W e make much use of the structure of flats in our proof, so it is helpful to describe some con ventions. Consider the set of unit-length tangent vectors at I , part of the tangent space S ( n ) ; in other words, the set of “directions” from I . If we choose a flat F to work in, then the tangent space of F contains a subset of those directions. All these directions, though, share the rotation Q identified with F . So in much of the rest of the paper , we refer to this rotation Q as a “direction, ” even though it is not a member of the tangent space. Our proof uses tw o k ey ideas. First, we show that within a flat F (i.e., given a direction Q ∈ SO ( n ) ) we can find a finite set of minimal horoballs exactly . This is done by showing an equiv alence between halfspaces in F and horoballs in P ( n ) in Section 4.1. The result implies that computing minimal horoballs with respect to a direction Q is equi valent to computing a con vex hull in Euclidean space. Second, we show that instead of searching ov er the entire space of directions SO ( n ) , we can discretize it into a finite set of directions such that when we calculate the horoballs with respect to each of these directions, the horoextents of the resulting ε -ball hull are not too far from those of the ball hull. In order to do this, we prove a Lipschitz bound for horofunctions (and hence horoextents) on the space of directions. Since any two flats F and F 0 are identified with rotations Q and Q 0 , we can move a point from F to F 0 simply by applying the rotation Q T Q 0 , and measure the angle θ between the flats. If we consider a geodesic c ⊂ F such that c (0) = I , we can apply Q T Q 0 to c to get c 0 , then for any point p ∈ P ( n ) we bound | b c ( p ) − b c 0 ( p ) | as a function of θ . 6 Proving this theorem is quite technical. W e first prov e a Lipschitz bound in P (2) , where the space of directions is a circle (as in the left part of Figure 1). After providing a bound in P (2) we decompose the distance between two directions in SO ( n ) into  n 2  angles defined by 2 × 2 submatrices in an n × n matrix. In this setting it is possible to apply the P (2) Lipschitz bound  n 2  times to get the full bound. The proof for P (2) is presented in Section 4.2, and the generalization to P ( n ) is presented in Section 4.3. Finally , we combine these results in an algorithm in Section 4.4. The follo wing lemma describes how geodesics (and horofunctions) are transformed by a rotation. Lemma 4.2. F or a point p ∈ P ( n ) , a r otation matrix Q , geodesics c ( t ) = e tA and c 0 ( t ) = e tQAQ T , b c 0 ( p ) = b c ( Q T pQ ) . Pr oof. If A 0 = QAQ T is the tangent vector of c 0 , and F 0 is the flat containing c 0 , b c 0 ( p ) = − tr( A 0 log( π F 0 ( p ))) = − tr(( QAQ T ) log (( Qν − 1 Q T ) p ( Qν − 1 Q T ) T )) = − tr( AQ T log( Qν − 1 ( Q T pQ ) ν − T Q T ) Q ) = − tr( A log ( ν − 1 ( Q T pQ ) ν − T )) = b c ( Q T pQ ) . In particular, this allo ws us to pick a flat where computation of b c is conv enient, and rotate the point set by Q to compute b c instead of attempting computation of b c 0 directly , which may be more cumbersome; we will utilize this idea later . 4.1 Projection to k -flat For the first part of our proof, we establish an equiv alence between horospheres and halfspaces. That is, after we compute the projection of our point set, we can say that the point set lies inside a horoball B r ( b c ) if and only if its projection lies inside a halfspace H r of F (recall that F is isometric to a Euclidean space under log ). Lemma 4.3. F or any hor oball B r ( b c ) , ther e is a halfspace H r ⊂ log ( F ) ⊂ S ( n ) suc h that log ( π F ( B r ( b c ))) = H r . Pr oof. If b c ( p ) ≤ r , p ∈ P ( n ) , and c ( t ) = e tA , then − tr( A log ( π F ( p ))) ≤ r . Since π F ( p ) is positiv e-definite, log( π F ( p )) is symmetric. But tr(( · )( · )) defines an inner product on the Euclidean space of symmetric n × n matrices. Then the set of all Y such that − tr( AY ) ≤ r defines a halfspace whose boundary is perpendicular to A . This gives us a means to compute horoballs by using π F to project our point set onto F , and lev erage a familiar Euclidean en vironment. 4.2 A Lipschitz b ound in P (2) 4.2.1 Rotations in P (2) W e start with some technical lemmas that describe the locus of rotating points in P (2) . Lemma 4.4. Given a r otation matrix Q ∈ SO (2) corr esponding to an angle of θ / 2 , Q acts on a point p ∈ P (2) via QpQ T as a r otation by θ about the (geodesic) axis e tI = e t I . Pr oof. If p = e t I , t ∈ R , then QpQ T = e t QI Q T = e t I , so e t I is in variant under the action of Q . Any action GpG T where G ∈ GL ( n ) is an isometry on P ( n ) , so the distance from p to the axis e t I remains fixed [6, II.10]. Computing QpQ T as a function of θ , we get:  cos θ 2 − sin θ 2 sin θ 2 cos θ 2   u w w v   cos θ 2 sin θ 2 − sin θ 2 cos θ 2  =  u + v 2 + u − v 2 cos θ − w sin θ u − v 2 sin θ + w cos θ u − v 2 sin θ + w cos θ u + v 2 − u − v 2 cos θ + w sin θ  , which is 2 π -periodic. (In fact, it is easy to see a rotation in the “coordinates” u − v 2 and w .) 7 By Lemma 4.4, we know that as we apply a rotation to p , it moves in a circle. Because any rotation Q has determinant 1 , det( QpQ T ) = det( p ) . This leads to the following corollary: Corollary 4.1. In P (2) , the radius of the circle that p travels on is GA( p ) . Such a cir cle lies entirely within a submanifold of constant determinant. In fact, any submanifold P (2) r of points with determinant equal to some r ∈ R + is isometric to any other such submanifold P (2) s for s ∈ R + . This is seen very easily by considering the distance function tr(log( p − 1 q )) — the determinants of p and q will cancel. W e pick a natural representativ e of these submanifolds, P (2) 1 . This submanifold forms a complete metric space of its o wn that has special structure: Lemma 4.5. P (2) 1 has constant sectional curvatur e − 1 2 . Pr oof. Let p ∈ P (2) 1 . Then if p =  x w w y  , det( p ) = xy − w 2 = 1 . Let u = x + y 2 and v = x − y 2 so that x = u + v and y = u − v . Then det( p ) = u 2 − v 2 − w 2 = 1 . This describes a hyperboloid of two sheets, and restricting u > 0 , is a model for hyperbolic space H 2 . T o analyze the metrics between the two spaces, we may consider our other point to be the identity matrix, since in P (2) 1 , d ( p, q ) = d ( q − 1 / 2 pq − 1 / 2 , I ) . The equi valent point on the hyperboloid to I is ( u, v , w ) = (1 , 0 , 0) . The distance between the two points in the hyperbolic metric is d H 2 ( p, I ) = cosh − 1 ( u 1 u 2 − v 1 v 2 − w 1 w 2 ) = cosh − 1  x + y 2  = ln   x + y 2 + s  x + y 2  2 − 1   . And in the metric of P ( n ) : d P (2) ( p, I ) = q tr(log 2 ( p )) = q ln 2 λ 1 + ln 2 λ 2 = √ 2 ln λ 1 = √ 2 ln   x + y 2 + s  x + y 2  2 − 1   . Since this is a constant multiple of the hyperbolic metric, P (2) 1 is a complete metric space of constant negati ve sectional curv ature. W e can find the curv ature κ by solving 1 / √ − κ = √ 2 to get κ = − 1 / 2 [6, I.2]. 4.2.2 Bounding k∇ b c k T o bound the error incurred by discretizing the space of directions, we need to understand the behavior of b c as a function of a rotation Q . W e show that the deri v ativ e of a geodesic is constant on P ( n ) . Lemma 4.6. F or a geodesic ray c ( t ) = e tA , k A k = 1 , then k∇ b c k = 1 at any point p ∈ P ( n ) . Pr oof. Let u = b c ( p ) . If we define the map γ p to map p to its projection onto any horoball B u − t ( b c ) , where t > 0 , then γ p is a geodesic ray [6, II.8]. Since any γ p ( t ) is the projection onto the horoball B u − t ( b c ) , the geodesic segment [ p, γ p ( t )] is perpendicular to B u − t ( b c ) at γ p ( t ) , and therefore the tangent v ector of γ p points directly opposite to ∇ b c at γ p ( t ) . Because γ p intersects B u − t ( b c ) at t , b c ( γ p ( t )) must change at a rate opposite to γ p ( t ) , along γ p , and since ∇ b c points in the opposite direction as γ 0 p ( t ) at γ p ( t ) , ∇ b c = − γ 0 p ( t ) . Also, since d ( p, B u − t ( b c )) = t for any p ∈ B u ( b c ) , for any u , and k γ 0 p k is constant along γ p , k∇ b c k is the same anywhere in P ( n ) . Since c ( t ) is the projection of c (0) onto B − t ( b c ) by construction, and geodesics are unique in P ( n ) , k∇ b c k = k A k = 1 . 8 4.2.3 A Lipschitz condition on Busemann functions in P (2) W e are no w ready to prove the main Lipschitz result in P (2) . W e start with a more specific bound that depends on the geodesic anisotropy of a point: Lemma 4.7. Given a point p ∈ P (2) , a r otation matrix Q corr esponding to an angle of θ / 2 , geodesics c ( t ) = e tA and c 0 ( t ) = e tQAQ T , | b c ( p ) − b c 0 ( p ) | ≤ | θ | · √ 2 sinh  GA( p ) √ 2  . Pr oof. The deriv ati ve of a function f along a curve γ ( t ) has the form  ∇ f | γ ( t ) , γ 0 ( t )  , and has greatest magnitude when the tangent vector γ 0 ( t ) to the curve and the gradient ∇ f | γ ( t ) are parallel. When this happens, the deriv ati ve reaches its maximum at k∇ f | γ ( t ) k · k γ 0 ( t ) k . Since k∇ b c k = 1 anywhere by Lemma 4.6, the deri vati ve of b c along γ at γ ( t ) is bounded by k γ 0 ( t ) k . W e are interested in the case where γ ( θ ) is the circle in P (2) defined by tracing Q ( θ / 2) pQ ( θ / 2) T for all − π < θ ≤ π . By Corollary 4.1, we know that this circle has radius GA( p ) and lies entirely within a submanifold of constant determinant, which by Lemma 4.5 also has constant curv ature κ = − 1 / 2 . This implies that k γ 0 ( θ ) k = 1 √ − κ sinh( √ − κ r ) = √ 2 sinh  GA( p ) √ 2  . for any v alue of θ ∈ ( − π , π ] [6, I.6]. Then | b c ( p ) − b c 0 ( p ) | = | b c ( p ) − b c ( Q T pQ ) | ≤ | θ | · √ 2 sinh  GA( p ) √ 2  . W e can no w state our main Lipschitz result in P (2) . Theorem 4.2. F or Q ∈ SO (2) corr esponding to θ / 2 , and let γ ( t ) = e tA and γ 0 ( t ) = e tQ T AQ . Then for any p ∈ X | b γ ( p ) − b γ 0 ( p ) | ≤ | θ | · √ 2 sinh  d X √ 2  . 4.3 Generalizing to P ( n ) No w to generalize to P ( n ) we need to decompose the projection operation π F ( · ) and the rotation matrix Q . W e can compute π F recursi vely , and it turns out that this f act helps us to break do wn the analysis of rotations. Since we can decompose any rotation into a series of 2 × 2 rotation matrices, decomposing the computation of π F in a similar manner lets us build a Lipschitz condition for P ( n ) . Lemma 4.8. Let F ⊂ P ( n ) be the flat containing diagonal matrices, let r + s = n , and let π F,r and π F,s be the pr ojection operations for r × r and s × s flats of diagonal matrices, r espectively . Then π F ( p ) =  π F,r ( h r ) π F,s ( p s )  , wher e p s is the lower right s × s block of p , and h r is the Schur complement of p s . Pr oof. π F is computed by decomposing p into ν − 1 f ν − T , where f is diagonal (and positi ve-definite) and ν is upper - triangular , with ones on the diagonal. This can be done by computing the Schur complement of some lower right square block of p , and putting the complement in the upper right block of a new matrix, with the lower block in the 9 other corner . The original matrix can be reconstructed by conjugating by an upper-triangular matrix of appropriate form: p =  p r a a T p s  =  I r ap − 1 s I s   h r p s   I r p − 1 s a T I s  . Performing this process recursiv ely yields a product of upper -triangular matrices with ones on the diagonal, which is again an upper -triangular matrix with ones on the diagonal, yielding ν , and a diagonal matrix, f = π F ( p ) . If we wish to compute π F ( p ) for other flats, then we can apply the rotation of F to p , compute using the recursi ve formula described above, and then apply the opposite rotation to the resulting diagonal matrix. Most of the time, ho wev er , it is most con venient to condition our data so that π F is computed for a diagonal flat F . W e no w wish to analyze a simpler form of rotation, one that can be broken into rotations on separate axes. Lemma 4.9. Given a point p ∈ P ( n ) , a r otation matrix Q =  Q r Q s  , such that r + s = n and Q r , Q s ar e r × r , s × s r otation matrices, r espectively , and a geodesic c ( t ) = e tQAQ T with A =  A r A s  sorted-diagonal, b c ( p ) = − tr( A r log( π F,r ( Q T r h r Q r ))) − tr( A s log( π F,s ( Q T s p s Q s ))) . Pr oof. From Lemma 4.2 we kno w that b c ( p ) = − tr( A log ( π F ( Q T pQ ))) . Compute π F ( Q T pQ ) by first decompos- ing p 7→  p r a a T p s  : Q T pQ =  Q T r Q T s   p r a a T p s   Q r Q s  =  Q T r p r Q r Q T r aQ s Q T s a T Q r Q T s p s Q s  . No w compute the Schur complement of Q T s p s Q s : Q T r p r Q r − Q T r aQ s ( Q T s p s Q s ) − 1 Q T s a T Q r = Q T r p r Q r − Q T r aQ s Q T s p − 1 s Q s Q T s a T Q r = Q T r p r Q r − Q T r ap − 1 s a T Q r = Q T r ( p r − ap − 1 s a T ) Q r . But h r = p r − ap − 1 s a T is just the Schur complement of p s , so b c ( p ) = − tr( A log ( π F ( Q T pQ ))) = − tr  A r log( π F,r ( Q T r h r Q r )) A s log( π F,s ( Q T s p s Q s ))  = − tr( A r log( π F,r ( Q T r h r Q r ))) − tr( A s log( π F,s ( Q T s p s Q s ))) . This allows us to break the Lipschitz bound into smaller pieces that we can analyze individually . The following two corollaries gi ve us a way to analyze the ef fects of 2 × 2 rotation matrices: Corollary 4.2. Given a point p ∈ P ( n ) , a r otation matrix Q =   I r Q 0 I s   , wher e r + s + 2 = n , Q 0 is a 2 × 2 r otation matrix corr esponding to an angle of θ / 2 , g eodesics c ( t ) = e tA and c 0 ( t ) = e tQAQ T , then | b c ( p ) − b c 0 ( p ) | is bounded as in Lemma 4.7. Pr oof. This is easily seen after observing that I r and I s are also rotation matrices, so Lemma 4.9 can be applied twice. 10 Corollary 4.3. The r esults of Cor ollary 4.2 extend to r otations between any two coordinates, that is, wher e Q 0 is of the form   cos( θ / 2) − sin( θ / 2) I t sin( θ / 2) cos( θ / 2)   . Pr oof. First observe that a rotation matrix Q i,j that rotates between axes i and j is equal to a matrix E i +1 ,j Q i,i +1 E T i +1 ,j , where E i,j is a permutation that mov es row i to row j and shifts the intervening rows up. W e assume E = E i +1 ,j , Q 0 = Q i,i +1 , and Q = E Q 0 E T from here on. Assuming that A is sorted-diagonal, we can compute b c 0 ( p ) as: b c 0 ( p ) = − tr(( E Q 0 E T ) A ( E Q 0 E T ) T log((( E Q 0 E T ) ν − 1 ( E Q 0 E T ) T ) p (( E Q 0 E T ) ν − 1 ( E Q 0 E T ) T ) T )) = − tr(( E T AE ) log (( E T ν − 1 E ) Q 0 T ( E T pE ) Q 0 ( E T ν − 1 E ) T )) = − tr( ˆ A log ( ˆ ν − 1 Q 0 T ˆ pQ 0 ˆ ν − T )) = − tr( ˆ A log ( π ˆ F ( Q 0 T ˆ pQ 0 ))) , which can be computed as abov e; some care must be taken, howe ver , since the order of elements of ˆ A is dif ferent than that of A . That is, in certain places, the Schur complement of the upper corner must be taken to compute π ˆ F , rather than that of the lo wer corner . 4.3.1 A Lipschitz condition on Busemann functions in P ( n ) W e are now ready to prove the main Lipschitz result in P ( n ) . W e start with a more specific bound that depends on the distance from a point p to I : Lemma 4.10. Given a point p ∈ P ( n ) , a r otation matrix Q ∈ SO ( n ) corresponding to an angle of θ / 2 , geodesics c ( t ) = e tA and c 0 ( t ) = e tQAQ T , | b c ( p ) − b c 0 ( p ) | ≤ | θ | ·  n 2  · √ 2 sinh  d ( p, I ) √ 2  . Pr oof. Every rotation Q may be decomposed into a product of rotations Q = Q 1 Q 2 . . . Q k where k =  n 2  and Q i is a 2 × 2 sub-block rotation corresponding to an angle of θ i / 2 with | θ i | ≤ | θ | . Then | b c ( p ) − b c 0 ( p ) | =      k X i =1 ( b i − 1 c 0 ( p ) − b i c 0 ( p ))      ≤ k X i =1 | b i − 1 c 0 ( p ) − b i c 0 ( p ) | , where b 0 c 0 ( p ) = b c ( p ) and b i c 0 ( p ) is b c ( p ) with the first i rotations successiv ely applied, so if Q 0 i = Q i j =1 Q j , b i c 0 ( p ) = b c (( Q 0 i ) T p ( Q 0 i )) . But then | b i − 1 c 0 ( p ) − b i c 0 ( p ) | ≤ | θ i | · √ 2 sinh  d ( p, I ) √ 2  , and therefore | b c ( p ) − b c 0 ( p ) | ≤ k X i =1 | θ i | ! · √ 2 sinh  d ( p, I ) √ 2  ≤ | θ | ·  n 2  · √ 2 sinh  d ( p, I ) √ 2  , since for all i we have | θ i | ≤ | θ | . 11 W e can no w state our main Lipschitz result in P ( n ) . Theorem 4.3 (Lipschitz condition on Busemann functions in P ( n ) ) . Consider a set X ⊂ P ( n ) , a r otation matrix Q ∈ SO ( n ) corresponding to an angle θ/ 2 , geodesics c ( t ) = e tA and c 0 ( t ) = e tQAQ T . Then for any p ∈ X | b c ( p ) − b c 0 ( p ) | ≤ | θ | ·  n 2  · √ 2 sinh  d X √ 2  . 4.4 Algo rithm For X ⊂ P ( n ) we can construct ε -ball hull as follows. W e place a grid G ε on SO ( n ) so that for any Q 0 ∈ SO ( n ) , there is another Q ∈ G ε such that the angle between Q and Q 0 is at most ( ε/ 2) / (2  n 2  √ 2 sinh( d X / √ 2)) . For each Q ∈ G ε , we consider π F ( X ) , the projection of X into the associated n -flat F associated with Q . W ithin F , we construct a con vex hull of π F ( X ) , and return the horoball associated with each hyperplane passing through each facet of the con vex hull, as in Lemma 4.3. T o analyze this algorithm we can now consider an y direction Q 0 ∈ SO ( n ) and a horofunction b c 0 that lies in the associated flat F 0 . There must be another direction Q ∈ G ε such that the angle between Q and Q 0 is at most ( ε/ 2) / (2  n 2  √ 2 sinh( d X / √ 2)) . Let b c be the similar horofunction to b c 0 , except it lies in the flat F associated with Q . This ensures that for any point p ∈ X , we have | b c 0 ( p ) − b c ( p ) | ≤ ε/ 2 . Since E c 0 ( X ) depends on two points in X , and each point changes at most ε/ 2 from b c 0 to b c we can argue that | E c 0 ( X ) − E c ( X ) | ≤ ε . Since this holds for any direction Q 0 ∈ SO ( n ) and for Q ∈ G ε the function E c ( X ) is exact, the returned set of horoballs defines an ε -ball hull. Since (with constant n ) the grid G ε is of size O ((sinh( d X ) /ε ) n − 1 ) and computing the con vex hull in each flat takes O ( N b n/ 2 c + N log N ) time this proves Theorem 4.1. 5 Center P oints In Euclidean space a center point p of a set X ⊂ R ¯ d of size N has the property that any halfspace that contains p also contains at least N/ ( d + 1) points from X . Center points alw ays exist [22] and there exists sev eral algorithms for computing them exactly [15] and approximately [10, 20]. W e cannot directly replicate the notion of center points in P ( n ) with horoballs. Instead we replace it with a slightly weaker notion, which is equiv alent in Euclidean space. A hor o-center point p of a set X ⊂ P ( n ) (or R ¯ d ) of size N has the property that any horoball that contains more than N d/ ( d + 1) points must contain p , where we define d = n ( n + 1) / 2 so that P ( n ) is a d -dimensional manifold. Construction fo r no center p oint in P ( n ) . Analogous to Euclidean center points, a center point p of X ⊂ P ( n ) of N points has the property that any horoball that contains p must also contain at least N / ( d + 1) points from X . Theorem 5.1. F or a set X ⊂ P ( n ) ther e may be no center point. Pr oof. Consider a set of distinct points X ∈ P ( n ) such that all points X lie on a single geodesic α between x 1 and x N where x 1 , x N ∈ X . Furthermore, let the points lie in a hyperbolic submanifold of P ( n ) . No w , for any point p not on α there is a horoball that contains p but contains none of X . So if there is a center point, it must lie on α . Ho wev er , also for any point p ∈ α there is a horoball that intersects α at only p , since the cross-section of any horoball in the hyperbolic submanifold will be strictly con vex (it can be represented as a hyperball in the Poincar ´ e model, and geodesics are circular arcs, so there is a horoball tangent to the geodesic at one point). Thus for any possible center point p there is a horoball that contains at most 1 point of X . Hence, there can be no center point. 12 Ho ro-center p oints in P ( n ) . The “simple” proof of the existence of center points [19] uses Helly’ s theorem to sho w that a horo-center point always exist, and then in Euclidean space a halfspace separation theorem can be used to show that a horo-center point is also a center point. W e replicate the first part in P ( n ) , but cannot replicate the second part because horoballs do not hav e the proper separation properties when not defined in R ¯ d . Theorem 5.2. Any set X ⊂ P ( n ) has a hor o-center point. Pr oof. W e use the follo wing Helly Theorem on Cartan-Hadamard manifolds (which include P ( n ) ) of dimension d [17]. For a family F of closed con ve x sets, if any set of d + 1 sets from F contain a common point, then the intersection of all sets in F contain a common point. In P ( n ) we consider the f amily F of closed con ve x sets defined as follo ws. A set F ∈ F is defined by a horofunc- tion b c and a subset X 0 ⊂ X of size greater than N d/ ( d + 1) such that X 0 is the intersection of X and a horoball B r ( b c ) . Then F = B ( X 0 ) is the ball hull of X 0 , so F ⊂ B r ( b c ) and F is compact. W e can ar gue that any set of d + 1 sets from F must intersect. W e can count the number of point s not in any d + 1 sets as S < d +1 X i =1 ( N − N d/ ( d + 1)) = d +1 X i =1 ( N (1 / ( d + 1)) = N . So there must be at least one point in X that is in all of the d + 1 sets. Then by the Helly-type theorem there exists a point p such that p ∈ F for any F ∈ F . W e can no w sho w that this point p must be a horo-center point. Any horoball that contains more than N d/ ( d + 1) points from X contains an element of F , thus it must also contain p . 5.1 Algo rithms for Ho ro-Center P oints W e provide justification for why it appears difficult to describe an exact algorithm for constructing horocenter points in P ( n ) and then provide an algorithm for an approximate horocenter point in P ( n ) . Before we begin we need a useful definition of a family of problems. An LP-type Pr oblem [24] takes as input a set of constraints H and a function ω : 2 H → R ¯ that we seek to minimize, and it has the follo wing two properties. M O N O T O N I C I T Y : For an y F ⊆ G ⊆ H , ω ( F ) ≤ ω ( G ) . L O C A L I T Y : For an y F ⊆ G ⊆ H with ω ( F ) = ω ( G ) and an h ∈ H such that ω ( G ∪ h ) > ω ( G ) implies that ω ( F ∪ h ) > ω ( F ) . A basis for an LP-type problem is a subset B ⊂ H such that ω ( B 0 ) < ω ( B ) for all proper subsets B 0 of B . And we say that B is a basis for a subset G ⊆ H if ω ( B ) = ω ( G ) and B is a basis. The cardinality of the largest basis is the combinatorial dimension of the LP-type problem. LP-type problems with constant combinatorial dimensions can be solved in time linear in the number of constraints [7]. Lemma 5.1. A set H of hor oballs in P ( n ) , and a function ω ( G ) = min p ∈ T H det( p ) is an LP-type pr oblem with constant combinatorial dimension. Pr oof. Monotonicity holds since in adding more horoballs to the a set F ⊂ H to get a set G ⊂ H (i.e. so F ⊂ G ) we hav e T G ⊆ T F . T o sho w locality , we consider subsets F ⊆ G ⊆ H such that ω ( F ) = ω ( G ) . Let P be the set of points { p ∈ P ( n ) | ω ( p ) = min q ∈ T G ω ( q ) } . Adding a constraint (a horoball) h to G only causes ω ( G ∪ h ) > ω ( G ) if P ∩ h = ∅ and thus P ∩ ( T G ∪ h ) = ∅ . Since T G ⊂ T F , then also P ∩ ( T F ∪ h ) = ∅ and ω ( F ∪ h ) > ω ( F ) . W e no w sho w that our problem has combinatorial dimension d = n ( n + 1) / 2 . P ( n ) is a d -dimensional manifold. Each constraint (a horosphere) is a ( d − 1) -dimensional sub-manifold of P ( n ) . Thus let p ∗ = arg min p ∈ T F ω ( p ) . If a constraint h lies in the the basis B ⊂ F , then p ∗ must lie on the corresponding horosphere. Hence, this reduces the problem by 1 dimension. And each subsequent horosphere h 0 we add to the basis, must also include p ∗ , so it must intersect h (and all other horospheres in the basis) transversally , reducing the dimension by 1 . (If h, h 0 ∈ B do not intersect transv ersally , then we can remo ve either one from the basis without changing p .) This process can only add d horospheres to B because P ( n ) is d -dimensional, thus the maximum basis size is d . 13 This lemma suggests the follo wing algorithm for constructing a horo-center point. Consider all subsets X 0 ⊂ X of N d/ ( d + 1) points, find the minimal horoball(s) which contain X 0 . Each of these horoballs can then be seen as a constraint for the LP-type problem. Then we solve the LP-type problem, returning a horo-center point. Unfortunately , there is no finite bound on the number of horoballs defined by a subset X 0 . Theorem 3.2 indicates that it could be infinite. Thus there are an infinite number of constraints that may need to be considered. In order to approximate the horocenter point, we use a similar approach as we did to approximate the ball hull. W e discretize the set of directions, and create a finite family of constraints for each direction. Then we can solve the associated LP-type problem to find a horo-center point. More formally , we place a grid G ε on SO ( n ) so that for any Q 0 ∈ SO ( n ) there is another Q ∈ G ε such that the angle between Q 0 and Q is at most ε/ (  n 2  2 √ 2 sinh( d X / √ 2)) . For each c corresponding to Q ∈ G ε , we consider π F ( X ) , the projection of X onto the ( d − 1) -flat F corresponding to Q . W ithin F , we can consider all subsets X 0 ⊂ X of N d/ ( d + 1) points and find the hyperplanes defining the con vex hull of π F ( X 0 ) . This finite set of hyperplanes corresponds to a finite set of horoballs which serve as constraints for the LP-type problem in P ( n ) . W e say a point ˆ p is an ε -appr oximate hor o-center point if there is a horo-center point p such that for any horo- function b c we hav e | b c ( p ) − b c ( ˆ p ) | ≤ ε . Lemma 5.2. A point ˆ p that satisfies all of the constraints defined by X and G ε is an ε -approximate hor o-center point of X . Pr oof. Let C ( X ) ⊂ P ( n ) be the set of horo-center points. Assume that ˆ p / ∈ C ( X ) , otherwise let p = ˆ p and we are done. W e sho w the existence of a specific nearby horo-center point p with the follo wing property . Let c be the geodesic connecting ˆ p and p , and let p = max p 0 ∈ C ( X ) b c ( p 0 ) . For any q ∈ ∂ C ( X ) let α be the geodesic connecting q and ˆ p . If q = arg max p 0 ∈ C ( X ) b α ( p 0 ) , we are done, if not, let q α ∈ ∂ C ( X ) such that q α = arg max p 0 ∈ C ( X ) b α ( p 0 ) . Then the geodesic ray on C ( X ) that goes from q to q α describes a flo w on C ( X ) . W e can see that the fixed point of this flo w is p , because as we move to ¯ q in the direction of this flo w , b α gets closer to is maximum, and the geodesic on P ( n ) from ˆ p to ¯ q is closer to direction defining the horofunction that ¯ q maximizes. No w , we can show that b c ( ˆ p ) − b c ( p ) = δ ≤ ε . This follo ws by Theorem 4.3 since there must be another direction Q ∈ G ε where the corresponding flat contains a geodesic c 0 such that there are more than N d/ ( d + 1) points x ∈ X such that b c 0 ( x ) < b c 0 ( ˆ p ) and | b c ( x ) − b c 0 ( x ) | ≤ ε . Thus there are more than N d/ ( d + 1) points x ∈ X such that b c ( x ) − ε < b c ( ˆ p ) , and there must be exactly b N d/ ( d + 1) + 1 c points x ∈ X such that b c ( x ) < b c ( p ) . Since ˆ p and p lie on the geodesic c , we have δ ≤ ε . W e can now use Lemma 4.6 to bound the difference in values for any horofunction. ||∇ b c || is constant for any b c , hence for any other horofunction b c 0 we hav e | b c 0 ( ˆ p ) − b c 0 ( p ) | ≤ | b c ( ˆ p ) − b c ( p ) | ≤ ε (where they are only equal when c and c 0 asymptote at opposite points). Since p is a horo-center point, this concludes the proof. When constructing the set of constraints in each flat F corresponding to a direction in G ε we do not need to explicitly consider all  N N d/ ( d +1)  subsets of X . Each constraint only depends on n points in X , so we can instead consider  N n  = O ( N n ) subsets of X of size n , and then check if either of the halfspaces it defines in F contain at least N d/ ( d + 1) points from X in O ( N ) time. Only these constraints need to be considered in the definition of ˆ p . Theorem 5.3. Given a set X ⊂ P ( n ) of size N , (for n constant) we can construct an ε -appr oximate hor o-center point in time O ((sinh( d X ) /ε ) n − 1 N n +1 ) time. Pr oof. As per the abov e construction, for each Q ∈ G ε we only need to consider O ( N n ) potential constraints, and each takes O ( N ) time to e valuate. By Lemma 4.7 G ε is of size O ((sinh( d X ) /ε ) n − 1 ) so the total number of constraints we need to consider in the LP-type problem is O ((sinh( d X ) /ε ) n − 1 N n ) . The total runtime is thus dominated by constructing the constraints and takes O ((sinh( d X ) /ε ) n − 1 N n +1 ) time. 14 References [1] A G A RW A L , P . K . , H A R - P E L E D , S . , A N D V A R A D A R A JA N , K . R . Approximating extent measures of points. J A CM 51 (2004). [2] B A R N E T T , V . The ordering of multiv ariate data. Journal of the Royal Statistical Society . Series A 139 (1976), 318 – 355. [3] B A S S E R , P . J . , M A T T I E L L O , J . , A N D L E B I H A N , D . MR diffusion tensor spectroscopy and imaging. Biophys. J . 66 , 1 (1994), 259–267. [4] B E R G E R , M . A P anoramic V iew of Riemannian Geometry . Springer , 2007. [5] B H A T I A , R . P ositive Definite Matrices . Princeton Uni versity Press, 2006. [6] B R I D S O N , M . R . , A N D H A E FL I G E R , A . Metric Spaces of Non-P ositive Curvatur e . Springer, 2009. [7] C H A Z E L L E , B . , A N D M A T O U ˇ S E K , J . On linear-time deterministic algorithms for optimization problems in fixed dimension. J ournal of Algorithms 21 (1996), 579–597. [8] C H E P O I , V . , D R AG A N , F . , E S T E L L O N , B . , H A B I B , M . , A N D V A X ` E S , Y . Diameters, centers, and approxi- mating trees of delta-hyperbolicgeodesic spaces and graphs. Annual Symposium on Computational Geometry (2008), 9. [9] C H E P O I , V . , A N D E S T E L L O N , B . Packing and Cov ering δ -Hyperbolic Spaces by Balls. Lecture Notes In Computer Science; V ol. 4627 (2007). [10] C L A R K S O N , K . L . , E P P S T E I N , D . , M I L L E R , G . L . , S T U RT I V A N T , C . , A N D T E N G , S . - H . Approximating center points with iterative radon points. International Journal of Computational Geometry and Applications 6 (1996), 357–377. [11] C OW I N , S . The structure of the linear anisotropic elastic symmetries. J ournal of the Mechanics and Physics of Solids 40 (1992), 1459–1471. [12] E P P S T E I N , D . Squarepants in a tree:Sum of subtree clustering and hyperbolic pants decomposition. A CM T ransactions on Algorithms (T ALG) 5 (2009). [13] F L E T C H E R , P . T . , A N D J O S H I , S . Principal geodesic analysis on symmetric spaces: Statistics of dif fusion tensors. In Computer V ision and Mathematical Methods in Medical and Biomedical Image Analysis (2004), pp. 87 – 98. [14] F L E T C H E R , P . T. , V E N K A TA S U B R A M A N I A N , S . , A N D J O S H I , S . The geometric median on riemannian manifolds with application to robust atlas estimation. Neur oImage 45 (2009), S143–52. [15] J A D H A V , S . , A N D M U K H O PA D H Y AY , A . Computing a center point of a finite planar set of points in linear time. Discr ete and Computational Geometry 12 (1994), 291–312. [16] K R AU T H G A M E R , R . , A N D L E E , J . R . Algorithms on negati vely curved spaces. FOCS (2006). [17] L E DY A E V , Y . S . , T R E I M A N , J . S . , A N D Z H U , Q . J . Helly’ s intersection theorem on manifolds of nonpositi ve curv ature. Journal of Comple x Analysis (2006). [18] L E I B O N , G . , A N D L E T S C H E R , D . Delaunay triangulations and V oronoi diagrams for Riemannian manifolds. In A CM Symposium on Computational Geometry (2000). [19] M A T O U ˇ S E K , J . Lectur es on Discr ete Geometry . Springer , 2002. 15 [20] M I L L E R , G . L . , A N D S H E E H Y , D . R . Approximate center points with proofs. In 25th Annual Symposium on Computational Geometry (2009). [21] M OA K H E R , M . , A N D B A T C H E L O R , P . G . Symmetric positi ve-definite matrices: From geometry to applica- tions and visualization. In V isualization and Pr ocessing of T ensor F ields . Springer , 2006. [22] R A D O , R . A theorem on general measure. J ournal of London Mathematical Society 21 (1947), 291–300. [23] S H A M O S , M . I . Geometry and statistics: Pr oblems at the interface . Computer Science Department, Carnegie- Mellon Uni versity , 1976. [24] S H A R I R , M . , A N D W E L Z L , E . A combinatorial bound for linear programming and related problems. In Pr oceedings 9th Annual Symposium on Theor etical Aspects of Computer Science (1992). [25] S H AW E - T A Y L O R , J . , A N D C R I S T I A N I N I , N . Kernel Methods for P attern Analysis . Cambridge Uni versity Press, 2004. 16