Hilbert geometry of the Siegel disk: The Siegel-Klein disk model

We introduce and study the Hilbert geometry induced by the Siegel disk, an open bounded convex set of complex matrices. This Hilbert geometry naturally yields a generalization of the Klein disk model of hyperbolic geometry, which we term the Siegel-Klein model to differentiate it with the usual Siegel upper plane and Siegel disk domains. In the Siegel-Klein disk, geodesics are by construction always straight, allowing one to build efficient geometric algorithms and data-structures from computational geometry. For example, we show how to approximate the Smallest Enclosing Ball (SEB) of a set of complex matrices in the Siegel domains: We compare two implementations of a generalization of the iterative algorithm of [Badoiu and Clarkson, 2003] in the Siegel-Poincar\'e disk and in the Siegel-Klein disk. We demonstrate that geometric computing in the Siegel-Klein disk allows one (i) to bypass the time-costly recentering operations to the origin (Siegel translations) required at each iteration of the SEB algorithm in the Siegel-Poincar\'e disk model, and (ii) to approximate numerically fast the Siegel distance with guaranteed lower and upper bounds.


Introduction
German mathematician Carl Ludwig Siegel [89] (1896-1981) and Chinese mathematician Loo-Keng Hua [45] (1910-1985) have introduced independently the symplectic geometry in the 1940's (with a preliminary work of Siegel [88] released in German in 1939). The adjective symplectic stems from the greek "complex", meaning mathematically the number field C instead of the ordinary real field R. Symplectic geometry 1 was originally motivated by the study of complex multivariate functions in the two landmark papers of Siegel [89] and Hua [45]. We refer the reader to the thesis [37,51] for an overview of Siegel domains. More generally, the Siegel domains have been studied and classified in the most general setting of bounded symmetric irreducible homogeneous domains of 6 types by E. Cartan [25] in 1935 (see also [50,14]).
The Siegel half-space and the Siegel disk provide generalizations of the complex Poincaré upper plane and the complex Poincaré disk to spaces of symmetric square complex matrices. We shall term them the Siegel-Poincaré upper plane and the Siegel-Poincaré disk in the remainder. The Siegel upper space includes the cone of real symmetric positive-definite (PD) matrices [35] (SPD manifold), and the well-known affine-invariant PD Riemannian metric [40] can be recovered as a special case of the Siegel metric.
In this paper, we extend the Klein disk model [83] of the hyperbolic geometry to the Siegel disk by considering the Hilbert geometry [43] induced by the open bounded convex Siegel disk [80,59]. We term this model the Klein-Siegel model for short to contrast with the Poincaré-Siegel upper plane and disk models. The main advantages of using the Klein-Siegel disk model instead of the usual Siegel upper plane or the Siegel-Poincaré disk are that the geodesics are always straight and therefore this Siegel-Klein disk model is very well-suited for algorithms and data-structures. Moreover, in the Siegel-Klein disk model, we have an efficient method to approximate with guarantees the calculation of the Siegel distance (specially useful when handling high-dimensional matrices). The algorithmic advantage was already observed for real hyperbolic geometry (included as a special case of the Siegel-Klein model): For example, the hyperbolic Voronoi diagrams can be efficiently computed as an affine power diagram clipped with the boundary ball [71,70,73,72]. To demonstrate the advantage of the Siegel-Klein disk over the Siegel-Poincaré disk, we consider approximating the smallest encloding ball of the a set of matrices in the disk. This problem has potential applications in image morphology [5,56] or anomaly detection of covariance matrices [95,29]. We state the problem as follows: Problem 1 (Smallest-radius Enclosing Ball (SEB)) Given a metric space (X, ρ) and a finite set {p 1 , . . . , p n } of n points, find the smallest-radius enclosing ball with center c * minimizing the following objective function: min c∈X max i∈{1,...,n} ρ(c, p i ).
In general, the SEBs may not be unique in a metric space: For example, the SEBs are not unique in a discrete Hamming metric space 2 [61]. The SEB is proven unique in the Euclidean metric space [98], the hyperbolic geometry [69], the Riemannian positive-definite matrix manifold 3 [55,68], and more generally in any Cartan-Hadamard manifold [6]. A fast (1 + )-approximation algorithm which requires 1 2 iterations was reported in [7,6]. Since the approximation factor does not depend on the dimension, this SEB approximation algorithm had many applications in machine learning [96] (e.g., in Reproducing Kernel Hilbert Spaces [81], RKHS).

Paper outline and contributions
In Section 2, we concisely recall the usual models of the hyperbolic complex plane: The Poincaré upper plane, and the Poincaré and Klein disk models. We then briefly review the Siegel upper plane in §3 and the Siegel disk in §4. Section 5 introduces the novel Siegel-Klein model using the Hilbert geometry. To demonstrate the algorithmic advantage of using the Siegel-Klein disk model, we compare two implementations of the Badoiu and Clarkson's approximation algorithm [7] extended to Siegel spaces in 6. Finally, we conclude this work in §7.
We list our main contributions as follows: • First, we formulate a generalization of the Klein disk model of hyperbolic geometry to the Siegel disk in Definition 2. We report the formula of the Siegel-Klein distance to the origin in Theorem 1 (and more generally the Siegel-Klein distance between two points whose line passes through the origin), describe how to convert the Siegel-Poincaré disk to the Siegel-Klein disk and vice versa in Proposition 2, report an exact algorithm to calculate the Siegel-Klein distance for diagonal matrices in Theorem 4. In practice, we show how to obtain a fast guaranteed approximation of the Siegel-Klein distance using bisection searches with guaranteed lower and upper bounds detailed in Theorem 5.
• Second, we report the exact solution to a geodesic cut problem in the Siegel-Poincaré/Siegel-Klein disks in Proposition 179. This result yields an explicit equation for the geodesic linking the origin of the Siegel disk to any other matrix point (Proposition 3 and Proposition 4). We then report an implementation of the Badoiu and Clarkson's iterative algorithm [7] for approximating the smallest enclosing ball tailored to the Siegel disk domains. In particular, we show in §6 that the implementation in the Siegel-Klein model yields a fast algorithm which bypasses the recentering operations required in the Siegel-Poincaré model.

Matrix spaces and matrix norms
Let F be a number field considered in the remainder to be either the real number field R or the complex number field C. For a complex number z = a + ib ∈ C, we denote z = a − ib its conjugate, and |z| = √ zz = √ a 2 + b 2 its modulus. Let Re(z) = a and Im(z) = b denote the real part and the imaginary part of the complex z, respectively.
Let M (d, F) be the space of d × d matrices with coefficients in F, and GL(d, F) denotes its subspace of invertible matrices. Let Sym(d, F) denote the space of d × d symmetric matrices with coefficients in F. The identity matrix is denoted by I (or I d when we want to emphasize its dimension d). A real matrix M ∈ M (d, R) is said positive-definite (PD) iff x M x > 0 for all x ∈ R d with x = 0. This positive-definiteness property is written M 0, where is the partial Löwner ordering [74]. Let PD(d, R) = {P 0 : P ∈ Sym(d, R)} be the space of real symmetric positivedefinite matrices [35,63,55,65] of dimension d × d. This space is a cone, i.e., if P 1 , P 2 ∈ PD(d, R) then P 1 +λP 2 ∈ PD(d, R) for all λ > 0. The boundary of the cone consists of rank-deficient positive semi-definite matrices.
The eigenvalues of a square complex matrix M are ordered such that |λ 1 (M )| ≥ . . . ≥ |λ d (M )|, where |·| denotes the complex modulus. The spectrum λ(M ) of a matrix M is its set of eigenvalues: λ(M ) = {λ 1 (M ), . . . , λ d (M )}. In general, real matrices may have complex eigenvalues but symmetric matrices (including PD matrices) have always real eigenvalues. The singular values σ i (M ) of M are always real: 4 Also denoted by the star operator (i.e., M * ) or the dagger (i.e., M † ) in the literature.
The induced Fröbenius distance between two complex matrices C 1 and C 2 is The operator norm or spectral norm of a matrix M is: = σ max (M ). To calculate the largest singular value, we may use a faster numerical power method of Lanczos [54] or its optimized variants [30]. 2 Hyperbolic geometry in the complex plane: The Poincaré and Klein models We concisely review the Poincaré upper plane, the Poincaré disk, and the Klein disk models of the hyperbolic plane [24,39]. In information geometry, the Fisher-Rao geometry of location-scale families amount to hyperbolic geometry [66].

Poincaré complex upper plane
The Poincaré upper plane domain is The Hermitian metric tensor is 4 or equivalently the Riemannian metric tensor is: Geodesics between z 1 and z 2 are arcs of semi-circles whose centers are on the real axis and orthogonal to the real axis, or vertical line segments when Im(z 1 ) = Im(z 2 ).
The geodesic distance is or equivalently where Equivalent formula can be obtained by using the identity log(x) = arcosh where By interpreting a complex number z = x + iy as a 2D point with Cartesian coordinates (x, y), the metric can be rewritten as where ds 2 E = dx 2 + dy 2 is the Euclidean metric. That is, the Poincaré upper plane metric can be rewritten as a conformal factor times the Euclidean metric. Thus the metric of Eq. 16 shows that the Poincaré upper plane model is a conformal model: That is, the Euclidean angle measurements in the (x, y) chart coincides with the underlying hyperbolic angles.
The group of orientation-preserving isometries (i.e., without reflections) is the real projective special group PSL(2, R) = SL(2, R)/{±I}, where SL(2, R) is the special linear group of matrices with unit determinant: The left group action is a fractional linear transformation (also called a Möbius transformation): The condition ab − cd = 0 is to ensure that the Möbius transformation is not constant. The set of Möbius transformations form a group Moeb(R, 2). The elements of the Möbius group can be represented by corresponding matrices of PSL(2, R):

5
The neutral element e is encoded by the identity matrix. The fractional linear transformations are the analytic mappings C ∪ {∞} → C ∪ {∞} of the Poincaré upper plane onto itself. The action is transitive (i.e., ∀z 1 , z 2 ∈ H, ∃g such that g.z 1 = z 2 ) and faithful (i.e., if g.z = z∀z then g = e). The stabilizer of i is the rotation group: The unit speed geodesic anchored at i and going up (geodesic with initial condition) is: Since the other geodesics are obtained by the action of PSL(2, R), it follows that the geodesics in H are parameterized by:

Poincaré disk
The Poincaré unit disk is The metric tensor is Since ds 2 D = 2 1− x 2 2 ds 2 E , we deduce that the metric is conformal. The geodesics between w 1 and w 2 are either arcs of circles intersecting orthogonally the disk boundary ∂D, or straight lines passing through the origin 0 of the disk and clipped to the disk domain.
The geodesic distance in the Poincaré disk is = 2arctanh The group of orientation preserving isometry is the complex projective special group PSL(2, C) = SL(2, C)/{±I} where SL(2, C) is the group of 2 × 2 complex matrices with unit determinant.
In the Poincaré disk model, the transformation corresponds to a hyperbolic motion (a Möbius transformation [82]) which moves point z 0 to the origin 0, and then makes a rotation of angle θ. The group of such transformations is the automorphism group of the disk, Aut(D), and transformation T z 0 ,θ is called a biholomorphic automorphism (a one-to-one conformal mapping of the disk onto itself).

Klein disk
The Klein disk model [24,83] (or Klein-Beltrami model) is defined on the unit disk as the Poincaré disk model. The metric is It is not a conformal metric (except at the origin), and therefore the Euclidean angles in the (x, y) chart do not correspond to the underlying hyperbolic angles. The Klein distance between two points k 1 = (x 1 , y 1 ) and k 2 = (x 2 , y 2 ) is (An equivalent formula will be reported in page 21 in a more general Theorem 4.) The advantage of the Klein disk over the Poincaré disk is that geodesics are straight Euclidean lines clipped to the unit disk domain. Therefore this model is well adapted to implement computational geometric algorithms and data structures, see for example [71,48].
The group of isometries in the Klein model are projective maps RP 2 preserving the disc.

Poincaré and Klein distances to the disk origin and conversions
In the Poincaré disk, the distance of a point w to the origin 0 is Since the Poincaré disk model is conformal (and Möbius transformations are conformal maps), Eq. 31 shows that Poincaré disks have Euclidean disk shapes (with displaced centers).
In the Klein disk, the distance of a point k to the origin is Observe the multiplicative factor of 1 2 in Eq. 32. Thus we can easily convert a point p ∈ C in the Poincaré disk to a point k ∈ C in the Klein disk, and vice-versa as follows: Let C K→D (k) and C D→K (w) denote these conversion functions with We can write C K→D (k) = α(k)k and C D→K (w) = β(w)w, so that α(k) > 1 is an expansion factor, and β(w) < 1 is a contraction factor.
The conversion functions are Möbius transformations represented by the matrices: For sanity check, let w = r + 0i be a point in the Poincaré disk with equivalent point k = 2 1+r 2 r + 0i in Poincaré disk. Then we have: We can convert a point z in the Poincaré upper plane to a corresponding point w in the Poincaré disk, or vice versa, using the following Möbius transformations: Notice that we compose Möbius transformations by multiplying their matrix representations.

The Siegel upper space
The Siegel upper space [88,89,64,37] SH(d) is defined as the space of symmetric complex matrices of size d × d which have positive-definite imaginary part: The space SH(d) is a tube domain of dimension d(d + 1). We can extract the components X and Y from Z as X = 1 2 (Z +Z) and . The pair (X, Y ) belongs to the Cartesian product of a vector space with the SPD cone: (X, Y ) ∈ Sym(d, R) × PD(d, R). When d = 1, the Siegel upper space coincides with the Poincaré upper plane: SH(1) = H. The geometry of the Siegel upper space was studied independently by Siegel [89] and Hua [45] from different viewpoints in the late 1930's-1940's. Historically, these class of complex matrices Z ∈ SH(d) were first studied by Riemann [85], and later eponymously named Riemann matrices. Riemann matrices are used to define Riemann theta functions [84,92,2,1].
The Siegel distance in the upper plane is induced by the following metric tensor: The formula for the Siegel upper distance between Z 1 and Z 2 ∈ SH(d) was calculated in Siegel's paper [89] as follows: where with R(Z 1 , Z 2 ) denoting the matrix generalization [19] of the cross-ratio 5 : and λ i (M ) denotes the i-th largest (real) eigenvalue of (complex) matrix M . This Siegel distance in the upper plane is a smooth spectral distance function: That is, is the eigenvalue map, and f is the following symmetric function (i.e., invariant under parameter permutations): A remarkable property is that all eigenvalues of R(Z 1 , Z 2 ) are positive (see [89]) although R may not necessarily be a Hermitian matrix. 6 Thus calculating the Siegel distance on the upper plane requires cubic time, i.e., cost of computing the eigenvalue decomposition. This Siegel distance in the upper plane SH(d) generalizes several distances: • When Z 1 = iY 1 and Z 2 = iY 2 , we have the Riemannian distance between Y 1 and Y 2 on the symmetric positive-definite manifold [35,63]: In that case, the Siegel upper metric for Z = iY becomes the affine-invariant metric: • In 1D, the Siegel upper distance ρ U (Z 1 , Z 2 ) between Z 1 = [z 1 ] and Z 2 = [z 2 ] (with z 1 and z 2 in C) amounts to the hyperbolic distance on the Poincaré upper plane H: where • The Siegel distance between two diagonal matrices Z = diag(z 1 , . . . , z d ) and Z = diag(z 1 , . . . , z d ) is: Observe that the Siegel distance is a non-separable metric distance, but its squared distance is separable when the matrices are diagonal: The Siegel metric in the upper plane is invariant by generalized matrix Möbius transformations (also called linear fractional transformations or rational transformation): where S ∈ M (2d, R) is the following 2d × 2d block matrix: which satisfies AB = BA , CD = DC , AD − BC = I.
The map φ S (·) = φ(S, ·) is called a symplectic map. The set of matrices S encoding the symplectic maps forms a group called the real symplectic group Sp(d, R) [37] (the group of Siegel motion): It can be shown that symplectic matrices have unit determinant [60,86], and therefore Sp(d, R) is a subgroup of SL(2d, R), the special group of real invertible matrices with unit determinant. We Matrix S denotes the representation of the group element g S . The symplectic group operation corresponds to matrix multiplications of their representations, the neutral element is encoded by Here, we use the parenthesis notation S (−1) to indicate that it is the group inverse and not the usual matrix inverse S −1 . The symplectic group is a Lie group of dimension d(2d + 1).
The action of the group is transitive: That is, for any Z = A + iB and S(Z) = we have φ S(Z) (iI) = Z. Therefore, by taking the group inverse we get The action φ S (Z) can be interpreted as a "Siegel translation" moving matrix iI to matrix Z, and conversely the action φ S (−1) (Z) as moving Z to iI.
The stabilizer group of Z = iI (also called isotropy group, the set of group elements S ∈ Sp(d, R) whose action fixes Z) is the subgroup of symplectic orthogonal matrices SpO(2d, R): We have SpO(2d, The geodesics in the Siegel upper space can be obtained by applying symplectic transformations to the geodesics of the positive-definite manifold (also known as SPD manifold) which is a totally geodesic submanifold of SU(d). Let Z 1 = iP 1 and Z 2 = iP 2 . Then the geodesic Z 12 (t) with Z 12 (0) = Z 1 and Z 12 (1) = Z 2 is expressed as: where Exp(M ) denotes the matrix exponential: and Log(M ) is the principal matrix logarithm, unique when M has all positive eigenvalues. The equation of the geodesic emanating from P with tangent vector S ∈ T p (symmetric matrix) on the SPD manifold is: Both the exponential and the principal logarithm of a matrix M can be calculated in cubic time when the matrices are diagonalizable: Let V denote the matrix of eigenvectors so that where λ 1 , . . . , λ d are the corresponding eigenvalues. Then for a scalar function f (e.g., f (u) = exp(u) or f (u) = log u), we define the corresponding matrix function f (M ) as

Siegel disk
The Siegel disk 7 [89] is an open convex complex domain defined by The Siegel disk can be written equivalently as SD(d) The boundary of the Siegel disk is called the Shilov boundary [28,37,36]). The Shilov boundary is a stratified manifold where each stratum is defined as a space of constant rank-deficient matrices [12].
The metric in the Siegel disk is: When d = 1, we recover ds 2 D = 1 (1−|w| 2 ) 2 dwdw which is the usual metric in the Poincaré disk (up to a missing factor of 4, see Eq. 25).
This Siegel metric induces a Kähler geometry [9] with Kähler potential: The distance between W 1 and W 2 in SD(d) is calculated as follows: where is a Siegel translation which moves W 1 to the origin O of the disk: We have Φ W (W ) = 0. Notice that the Siegel disk distance, although a spectral distance function via the operator norm, is not smooth because of it uses the maximum singular value (recall that the Siegel upper plane distance uses all eigenvalues of a matrix cross-ratio R).
It follows that the cost of calculating a Siegel distance in the Siegel disk is cubic: We require to compute a symmetric matrix square root [91] in Eq. 77, and then compute the largest singular value for the operator norm in Eq. 76.
Notice that when d = 1, the "1d" scalar matrices commute, and we have: This corresponds to a translation of w 1 to 0 (see Eq. 28). Let us observe the following special cases of the Siegel-Poincaré distance: • Distance to the origin: When W 1 = 0 and W 2 = W , we have Φ 0 (W ) = W , and therefore the Siegel distance in the disk between a matrix W and the origin 0 is: In particular, when d = 1, we recover the formula of Eq. 31: ρ D (0, w) = log 1+|w| 1−|w| .
Thus the diagonal matrices belong to the polydisk domain. Then we have Notice that the polydisk domain is a Cartesian product of 1D complex disk domains, but it is not the unit d-dimensional complex ball {z ∈ C d : We can convert a matrix Z in the Siegel upper space to an equivalent matrix W in the Siegel disk by using the following matrix Cayley transformation for Z ∈ SH d : Notice that the imaginary positive-definite matrices iP of the upper plane (vertical axis) are mapped to i.e., the real symmetric matrices belonging to the horizontal-axis of the disk. The inverse transformation for a matrix W in the Siegel disk is a matrix in the Siegel upper space.

13
With those mappings, the origin of the disk 0 ∈ SD(d) coincides with iI ∈ SH(d).
A key property is that the geodesics passing through the matrix origin 0 are expressed by straight line segments in the Siegel disk. We can check that for any α ∈ [0, 1].
To describe the geodesics between W 1 and W 2 , we first move W 1 to 0 and W 2 to Φ W 1 (W 2 ). Then the geodesic between 0 and Φ W 1 (W 2 ) is a straight line segment, and we map back this geodesic via Φ −1 W 1 (·). The inverse of a symplectic map is a symplectic map which corresponds to the action of an element of the complex symplectic group.
The complex symplectic group is with for the d × d identity matrix I. Notice that the condition M JM = J amounts to check that The conversions between the Siegel upper plan to the Siegel disk (and vice versa) can be expressed using complex symplectic transformations associated to the matrices: It can be shown that with and the left action of g ∈ Sp(d, C) is iI The isotropy group at the origin 0 is where U (d) is the unitary group: Thus we can "rotate" a matrix W with respect to the origin so that its imaginary part becomes 0: There exists A such that Re(AW W −1Ā−1 ) = 0.
More generally, we can define a Siegel rotation [62] in the disk with respect to a center W 0 ∈ SD(d) as follows: whereĀ In 1D, the Poincaré disk can be embedded non-diagonally onto the Siegel upper plane [87].
wherep andq are the unique two intersection points of the line (pq) with the boundary ∂Ω of the domain Ω as depicted in Figure 2, and CR denotes the cross-ratio of four points (a projective invariant): When p = q, we have: The Hilbert distance is a metric distance which does not depend on the underlying norm of the vector space: Proposition 1 (Formula of Hilbert distance) The Hilbert distance between two points p and q of an open bounded convex domain Ω is Proof: For distinct points p and q of Ω, let α + > 1 be such thatq = p + α + (q − p), and α − < 0 We may also write the source points p and q as linear interpolations of the extremal pointsp andq on the boundary: p = (1 − α p )p + α pq and q = (1 − α q )p + α qq with 0 < α p < α q < 1. In that case, the Hilbert distance can be written as The space (Ω, H Ω ) is a metric space. Notice that the above formula has demonstrated that H Ω,κ (p, q) = H Ω∩(pq),κ (p, q). That is, the Hilbert distance between two points of a d-dimensional domain is equivalent to the Hilbert distance between the two points on the 1D domain defined by Ω restricted to the line (pq) passing through the points p and q.
Notice that the boundary ∂Ω of the domain may not be smooth (e.g., Ω may be a simplex or a polytope). The Hilbert geometry for a unit disk with κ = 1 2 yields the Klein model [49] (or Klein-Beltrami model [13]) of hyperbolic geometry. The Hilbert geometry for an ellipsoid yields the Cayley-Klein hyperbolic model [26,83,70] generalizing the Klein model. The Hilbert geometry for a simplicial polytope is isometric to a normed vector space [32,76]. We refer to the handbook [80] for a survey with recent results on Hilbert geometry. The Hilbert geometry of the elliptope (i.e., space of correlation matrices) was studied in [76]. Hilbert geometry may be studied from the viewpoint of Finslerian geometry which is Riemannian if and only if the domain Ω is an ellipsoid (i.e., Klein and Cayley-Klein hyperbolic geometries). When d = 1, the Siegel-Klein disk is the Klein disk model of hyperbolic geometry, and the Klein distance [71] between two any points k 1 ∈ C and k 2 ∈ C restricted to the unit disk is

Hilbert geometry of the Siegel disk domain
where This formula can be derived from the Hilbert distance induced by the Klein unit disk [83].

Calculating the Siegel-Klein distance
The Siegel disk domain SD(d) = W ∈ Sym(d, C) : I − W W 0 can be rewritten using the operator norm as Let {K 1 + α(K 2 − K 1 ), α ∈ R} denote the line passing through (matrix) points K 1 and K 2 . That line intersects the Shilov boundary when When K 1 = K 2 , there are two unique solutions 9 : Let one solution be α + with α + > 1, and the other solution be α − with α − < 0. The Siegel-Klein distance is then defined as LetK 1 = K 1 + α − (K 2 − K 1 ) andK 2 = K 1 + α + (K 2 − K 1 ) be the extremal matrices (rank deficient belonging to the Shilov boundary).
In practice, we may perform a bisection search on the line (K 1 K 2 ) to approximate these two extremal pointsK 1 andK 2 (such that these matrices are ordered along the line asK 1 , K 1 , K 2 , K 2 ). We may find a lower bound for α − and a upper bound for α + as follows: We seek α on the line (K 1 K 2 ) such that K 1 + α(K 2 − K 1 ) falls outside the Siegel disk: Since · O is a matrix norm, we have Thus we deduce that

Siegel-Klein distance to the origin
When K 1 = 0 (the 0 matrix denoting the origin of the disk), and K 2 = K ∈ SD(d), it is easy to solve: We have |α| = 1 K O , that is, In that case, the Siegel-Klein distance of Eq. 112 is expressed as: where ρ D (0, W ) is defined in Eq. 80.

Theorem 1 (Siegel-Klein distance to the origin) The Siegel-Klein distance of matrix
The constant κ = 1 2 is chosen in order to ensure that when d = 1 the corresponding Klein disk has negative unit curvature.
The result can be easily extended to the case of the Siegel-Klein distance between K 1 and K 2 where the origin O belongs to the line (K 1 K 2 ). In that case, K 2 = λK 1 for some λ ∈ R (e.g., λ = tr(K 2 ) tr(K 1 ) ). It follows that Thus we get the two values defining the intersection of (K 1 K 2 ) with the Shilov boundary: We then apply formula Eq. 112: Theorem 2 The Siegel-Klein distance between two points K 1 = 0 and K 2 on a line (K 1 K 2 ) passing through the origin is where λ == tr(K 2 ) tr(K 1 ) .

Converting Siegel-Poincaré matrices ⇔ Siegel-Klein matrices
From Eq. 122, we deduce that we can convert a matrix K in the Siegel-Klein disk to a corresponding matrix W in the Siegel-Poincaré disk, and vice versa, as follows: • Converting K to W : We convert a matrix K in the Siegel-Klein model to an equivalent matrix W in the Siegel-Poincaré model as follows: • Converting W to K: We convert a matrix W in the Siegel-Poincaré model to an equivalent matrix K in the Siegel-Klein model as follows: (130)

Proposition 2 (Conversions Siegel-Poincaré⇔Siegel-Klein disk)
The conversion of a matrix K of the Siegel-Klein model to its equivalent matrix W in the Siegel-Poincaré model, and vice versa, is done by the following radial expansion and contraction functions: By virtue of the cross-ratio property 10 , the (pre)geodesics 11 in the Hilbert-Klein disk are Euclidean straight. Thus we can write the pregeodesics as: Another way to get a generic closed-form formula for the Siegel-Klein distance is by using the formula for the Siegel-Poincaré disk after converting the matrices to their equivalent matrices in the Siegel-Poincaré disk. We get the following expression: Theorem 3 (Formula for the Siegel-Klein distance) The Siegel-Klein distance between K 1 and K 2 in the Siegel disk is ρ K (K 1 , K 2 ) = 1 2 log The isometries in Hilbert geometry have been studied in [90]. 10 The cross-ratio (p, q; P, Q) = p−P q−Q p−Q q−P of four collinear points on a line is such that (p, q; P, Q) = (p, r; P, Q) × (r, q; P, Q) whenever r belongs to that line. 11 Geodesics are paths which minimize locally the distance and are parameterized proportionally to the arc-length. A pregeodesic is a path which minimizes locally the distance but is not necessarily parameterized proportionally to the arc-length. 20

Siegel-Klein distance between diagonal matrices
Solving for the general case: Let K α = K 1 + αK 21 with K 21 = K 2 − K 1 . We seek for the extremal values of α such that This last equation is reminiscent to a Linear Matrix Inequality [33] (LMI, i.e., i y i S i 0 with y i ∈ R and S i ∈ Sym(d, R) where the coefficients y i are however linked between them).
Let us consider the special case of diagonal matrices corresponding to the polydisk domain: K = diag(k 1 , . . . , k d ) and K = diag(k 1 , . . . , k d ) of the disk.
First, let us start with the simple case d = 1, i.e., the Siegel disk SD(1) which is the complex open unit disk {k ∈ C :kk < 1}. Let k α = (1 − α)k 1 + αk 2 = k 1 + αk 21 with k 21 = k 2 − k 1 . We havek α k α = aα 2 + bα + c with a =k 2 1k 21 , b =k 1 k 21 +k 21 k 1 and c =k 1 k 1 . To find the two intersection points of line (k 1 k 2 ) with the boundary of SD(1), we need to solvek α k α = 1. This amounts to solve an ordinary quadratic equation since all coefficients a, b, and c are provably reals. Let ∆ = b 2 −4ac be the discriminant (∆ > 0 when k 1 = k 2 ). We get the two solutions , and apply the 1D formula for the Hilbert distance: Doing so, we obtain a formula equivalent to Eq. 30. For diagonal matrices with d > 1, we get the following system of d inequalities: For each inequality, we solve the quadratic equation as in the 1d case above, yielding two solutions α − i and α + i . Then we satisfy all those constraints by setting and we compute the Hilbert distance: Theorem 4 (Siegel-Klein distance for diagonal matrices) The Siegel-Klein distance between two diagonal matrices in the Siegel-Klein disk can be calculated exactly in linear time.
Notice that the proof extends to triangular matrices as well.
When the matrices are non-diagonal, we have to solve analytically the equation: such that with the following symmetric matrices: When S 0 , S 1 and S 2 are simultaneously diagonalizable via congruence [21], the optimization problem becomes: such that where D i = P S i P for some P ∈ GL(d, C), and we apply Theorem 4. The same result applies for simultaneously diagonalizable matrices S 0 , S 1 and S 2 via similarity: Notice that the Hilbert distance (or its squared distance) is not a separable distance, even in the case of diagonal matrices. (But recall that the squared Siegel-Poincaré distance in the upper plane is separable for diagonal matrices.)

A fast guaranteed approximation of the Siegel-Klein distance
In the general case, we use the bisection approximation algorithm which is a geometric approximation technique that requires to only calculate operator norms (and not the square root matrices required in the functions Φ · (·)).
We have the following important property of the Hilbert distance: Property 1 (Bounding Hilbert distance) Let Ω + ⊂ Ω ⊂ Ω − be strictly nested open convex bounded domains. Then we have the following inequality for the corresponding Hilbert distances: See Figure 5 for a visual proof. Figure 4 illustrates the Property 1 of Hilbert distances corresponding to nested domains. Notice that when Ω − is a large enclosing ball of Ω with radius increasing to infinity, we have α − α + , and therefore the Hilbert distance tends to zero.
Therefore the bisection search for finding the values of α − and α + yields lower and upper bounds on the exact Siegel-Klein distance as follows: Let α − ∈ (l − , u − ) and α + ∈ (l + , u + ) where l − , u − , l + , u + are reals defining the extremities of the intervals. Using Property 1, we get the following theorem:

Theorem 5 (Lower and upper bounds on the Siegel-Klein distance) The
Siegel-Klein distance between two matrices K 1 and K 2 of the Siegel disk is bounded as follows: where Figure 6 depicts the guaranteed lower and upper bounds obtained by performing the bisection search for approximating the pointK 1 ∈ (K 1 ,K 1 ) and the pointsK 2 ∈ (K 2 ,K 2 ).
We have: and hence H Ω , 1 Notice that the approximation of the Siegel-Klein distance by line bisection requires only to calculate an operator norm M O at each step: This involves calculating the smallest and largest eigenvalues of M , or the largest eigenvalue of MM . To get a (1 + )-approximation, we need to perform O(log 1 ) dichotomic steps. This yields a fast method to approximate the Siegel-Klein distance compared with the exact calculation of the Siegel-Klein distance of Eq. 132 which requires to calculate Φ · (·) functions: This involves the calculation of a square root of a complex matrix. Notice that the operator norm can be numerically approximated using a Lanczos's power iteration scheme [53,42] (see also [57]).

Hilbert-Fröbenius distances and fast simple bounds on the Siegel-Klein distance
Let us notice that although the Hilbert distance does not depend on the chosen norm in the vector space, but the Siegel complex ball SD(d) is defined according to the operator norm. In a finitedimensional vector space, all norms are equivalent: That is, given two norms · a and · b of Figure 6: Guaranteed lower and upper bounds for the Siegel-Klein distance.
vector space X, there exists positive constants c 1 and c 2 such that In particular, this property holds for the operator norm and Fröbenius norm of complex matrices with positive constants c d , C D , c d and C d depending on the dimension d of the square matrices: Therefore we have where H FD(d,r), 1 2 denotes the Fröbenius-Klein distance, i.e., the Hilbert distance induced by the Fröbenius balls FD(d, r) with constant κ = 1 2 . Now, we can calculate in closed-form the Fröbenius-Klein distance by computing the two intersection points of the line (K 1 K 2 ) with the Fröbenius ball FD(d, r). This amounts to solve an ordinary quadratic equation K 1 + α(K 2 − K 1 ) 2 F = r for parameter α: where K i,j denotes the coefficient of matrix K at row i and column j. Notice that is a real. Once α − and α + are found, we apply the 1D formula of the Hilbert distance of Eq. 107.
We summarize the result as follows: 6 The smallest enclosing ball in the SPD manifold and in the Siegel spaces The goal of this section is to compare two implementations of a generalization of the Badoiu and Clarkson's algorithm [7] to approximate the smallest enclosing ball of a set of complex matrices either in the Siegel-Poincaré disk or in the Siegel-Klein disk.
In general, we may encode a pair of features (S, P ) ∈ Sym(d, R) × P ++ (d, R) in applications as a Riemann matrix Z(S, P ) := S + iP , and consider the underlying symplectic geometry. For example, anomaly detection of time-series maybe considered by considering (Σ(t), Σ(t)) where Σ(t) is the covariance matrix at time t andΣ(t) 1 dt (Σ(t + dt) − Σ(t)) is the approximation of the derivative of the covariance matrix (a symmetric matrix).
The generic Badoiu and Clarkson's algorithm [7] (BC algorithm) for a set {p 1 , . . . , p n } of n p points in a metric space (X, ρ) is described as follows: • Initialization: Let c 1 = p 1 and l = 1 • Repeat L times: -Calculate the farthest point: f l = arg min i∈[d] ρ(c l , p i ).
-Geodesic cut: Let c l+1 = c l # t l f t , where p# t l q is the point which satisfies ρ(p, p# X t l q) = t l ρ(p, q).
This elementary SEB approximation algorithm has been instantiated to various metric spaces with proofs of convergence according to the sequence 12 {t l } l : see [69] for the case of hyperbolic geometry, [6] for Riemannian geometry with bounded sectional curvatures, [77] for dually flat spaces (a non-metric space equipped with a Bregman divergences), etc. The number of iterations L to get a (1+ )-approximation of the SEB depends on the underlying geometry and the sequence {t l } l . For example, in Euclidean geometry, setting t l = 1 l+1 with L = 1 2 steps yield a (1 + )-approximation of the SEB [7].
We start by recalling the Riemannian generalization of the BC algorithm, and then consider the Siegel spaces.

Approximating the SEB in Riemannian spaces
We first instantiate a particular example of Riemannian space, the space of Symmetric Positive-Definite matrix manifold (PD/SPD manifold), and then consider the general case on a Riemannian manifold (M, g).

Approximating the SEB on the SPD manifold
Given n positive-definite matrices [17,34] P 1 , . . . , P n of size d × d, we ask to calculate the SEB with circumcenter P * minimizing the following objective function: This is a minimax optimization problem. The SPD cone is not a complete metric space with respect to the Fröbenius distance, but is a complete metric space with respect to the natural Riemannian distance. When the minimization is performed with respect to the Fröbenius distance, we can solve this problem using techniques of Euclidean computational geometry [16,7] by vectorizing the PSD matrices P i into corresponding vectors v i = vec(P i ) of R d×d such that P − P F = vec(P ) − vec(P ) 2 , where vec(·) : Sym(d, R) → R d×d vectorizes a matrix by stacking its column vectors. In fact, since the matrices are symmetric, it is enough to half-vectorize the matrices: , see [74].

Property 2
The smallest enclosing ball of a finite set of positive-definite matrices is unique.
• Another proof of the uniqueness of the SEB on a SPD manifold consists in noticing that the SPD manifold is a Cartan-Hadamard manifold [6], and the SEB on Cartan-Hadamard manifolds are guaranteed to be unique.
We shall use the invariance property of the Riemannian distance by congruence: In particular, choosing C = P − 1 2 1 , we get The geodesic from I to P is γ I,P (α) = Exp(αLogP ) = P α . The set {λ i (P α )} of the d eigenvalues of P α coincide with the set {λ i (P ) α } of eigenvalues of P raised to the power α (up to a permutation). Thus to cut the geodesic I# PD t P , we have to solve the following problem: ρ PD (I, P α ) = t × ρ PD (I, P ).

26
That is The solution is α = t. Thus I# PD t P = P t . For arbitrary P 1 and P 2 , we first apply the congruence transformation with C = P Theorem 7 (Geodesic cut on the SPD manifold) For any t ∈ (0, 1), we have the closed-form expression of the geodesic cut on the manifold of positive-definite matrices: The matrix P can be rewritten using the orthogonal eigendecomposition as U DU , where D is the diagonal matrix of generalized eigenvalues. Thus the PD geodesic can be rewritten as We instantiate the generic algorithm to positive-definite matrices as follows: Algorithm ApproximatePDSEB({P 1 , . . . , P n }, L): • Initialization: Let C 1 = P 1 and l = 1 • Repeat L times:l -Calculate the index of the farthest matrix: f l = arg min i∈{1,...,d} ρ PD (C t , P i ).
-Geodesic walk: The complexity of the algorithm is in O(d 3 nT ) where T is the number of iterations, d the row dimension of the square matrices P i and n the number of matrices.
Observe that the solution corresponds to the arc-length parameterization of the geodesic with boundary values on the SPD manifold: In fact, we have shown the following property: Property 3 (Riemannian geodesic cut) Let γ p,q (t) denote the Riemannian geodesic linking p and q on a Riemannian manifold (M, g) (i.e., parameterized proportionally to the arc-length and with respect to the Levi-Civita connection induced by the metric tensor g). Then we have Thus it follows the following generic Riemannian algorithm: Algorithm ApproximateRieSEB({p 1 , . . . , p n }, g, L): • Initialization: Let c 1 = p 1 and l = 1 • Repeat L times: -Calculate the index of the farthest point: f l = arg min i∈{1,...,d} ρ g (c l , p i ).
Theorem 1 of [6] guarantees the convergence of the ApproximateRieSEB algorithm provided that we have a lower bound and an upper bound on the sectional curvatures of the manifold (M, g). The sectional curvatures of the PD manifold have been proven to be negative [41]. The SPD manifold is a Cartan-Hadamard manifold with scalar curvature 1 8 d(d + 1)(d + 2) [3] depending on the dimension d of the matrices. Notice that we can identify P ∈ PD(d) with an element of the quotient space GL(d, R)/O(d) since O(d) is the isotropy subgroup of the GL(d, R) for the action P → C T P C (i.e., I → C T IC = I when C ∈ O(d)). Thus we have PD(d) ∼ = GL(d, R)/O(d).

Implementation in the Siegel-Poincaré disk
Given n d × d complex matrices W 1 , . . . , W n ∈ SD(d), we ask to find the smallest-radius enclosing ball with center W * minimizing the following objective function: This problem may have potential applications in image morphology [4] or anomaly detection 13 of covariance matrices [95].
The Siegel-Poincaré upper plane and disk are not Bruhat-Tits space 14 , but spaces of non-positive curvatures [31].
Notice that when d = 1, the hyperbolic ball in the Poincaré disk have Euclidean shape. This is not true anymore when d > 1: Indeed, the equation of the ball centered at the origin 0: amounts to When d = 1, W O = |w| = (Re(w), Im(w)) 2 , and Poincaré balls have Euclidean shapes. Otherwise, when d > 1, W O = σ max (W ) and σ max (W ) ≤ e r −1 e r +1 is not a complex Fröbenius ball. In order to apply the generic algorithm, we need to implement the geodesic cut operation W 1 # t W 2 . We consider the complex symplectic map Φ W 1 (W ) in the Siegel disk that maps W 1 to 0 and W 2 to W 2 = Φ W 1 (W 2 ). Then the geodesic between 0 and W 2 is a straight line.
We need to find α(t)W = 0# SD t W (with α(t) > 0) such that ρ D (0, α(t)W ) = tρ D (0, W ). That is, we shall solve the following equation: We find the exact solution as Proposition 3 (Siegel-Poincaré geodesics from the origin) The geodesic in the Siegel disk is Thus the midpoint W 1 # SD W 2 := W 1 # SD 1 2 W 2 of W 1 and W 2 can be found as follows: where To summarize, the algorithm recenters at every step the current center C t to the Siegel disk origin 0: Algorithm ApproximateSiegelSEB({W 1 , . . . , W n }): • Initialization: Let C 1 = 0 and l = 1.
• Repeat L times: -Calculate the index of the farthest point: F l = arg min i∈[d] ρ D (0, W i ).
• Let the approximate circumcenter be mapped back to be consistent with the input: This amounts to calculate the symplectic map associated to the matrix S = C The farthest point to the current approximation of the circumcenter can be calculated using the data-structure of the Vantage Point Tree (VPT), see [75].
The Riemannian curvature tensor of the Siegel space is non-positive [89,44] and the sectional curvatures are non-positive [31] and bounded above by a negative constant. In our implementation, we chose the step sizes t l = 1 l+1 . Barbaresco [11] also adopted this iterative recentering operation for calculating the median in the Siegel disk. However at the end of his algorithm, he does not map back the median among the source matrix set. Recentering is costly because we need to calculate a square root matrix to calculate Φ C (W ). A great advantage of Siegel-Klein space is that we have straight geodesics anywhere in the disk so we do not need to perform recentering.

Fast implementation in the Siegel-Klein disk
The main advantage of implementing the Badoiu and Clarkson's algorithm [7] in the Siegel-Klein disk is to avoid to perform the costly recentering operations (which require calculation of square root matrices). Moreover, we do not have to roll back our approximate circumcenter at the end of the algorithm.
First, we state the following expression of the geodesics in the Siegel disk: Proposition 4 (Siegel-Klein geodesics from the origin) The geodesic from the origin in the Siegel-Klein disk is expressed γ SK 0,K (t) = α(t)K with The proof follows straightforwardly from Proposition 3 because we have ρ K (0, K) = 1 2 ρ D (0, K).

30
In this paper, we have generalized the Klein Since the geodesics in the Siegel-Klein disks are by construction straight, this model is wellsuited to implement techniques of computational geometry [16]. Furthermore, the calculation of the Siegel distance in the Siegel-Klein disk does not require to recenter one of its arguments to the disk origin, a costly operation. We reported a linear-time algorithm for computing the exact Siegel-Klein distance between diagonal matrices of the disk (Theorem 4), and a fast way to numerically approximate the Siegel distance by bisection searches with guaranteed lower and upper bounds (Theorem 5). Finally, we demonstrated the algorithmic advantage of using the Siegel-Klein disk model instead of the Siegel-Poincaré disk model for approximating the smallest-radius enclosing ball of a finite set of matrices in the Siegel disk. In future work, we shall consider more generally the Hilbert geometry of homogeneous complex domains applications and applications of the Siegel-Klein geometry for radar processing [11], image morphology [56], computer vision, and machine learning [52].

Snippet code
We implemented our software library and smallest enclosing ball algorithms in Java TM . The code below is a snippet written in Maxima: A computer algebra system, freely downloadable at http://maxima.sourceforge.net/