HoroPCA: Hyperbolic Dimensionality Reduction via
Horospherical Projections

Ines Chami Equal contribution. Albert Gu^∗ Dat Nguyen^∗ Christopher Ré

Abstract

This paper studies Principal Component Analysis (PCA) for data lying in hyperbolic spaces. Given directions, PCA relies on: (1) a parameterization of subspaces spanned by these directions, (2) a method of projection onto subspaces that preserves information in these directions, and (3) an objective to optimize, namely the variance explained by projections. We generalize each of these concepts to the hyperbolic space and propose HoroPCA, a method for hyperbolic dimensionality reduction. By focusing on the core problem of extracting principal directions, HoroPCA theoretically better preserves information in the original data such as distances, compared to previous generalizations of PCA. Empirically, we validate that HoroPCA outperforms existing dimensionality reduction methods, significantly reducing error in distance preservation. As a data whitening method, it improves downstream classification by up to 3.9% compared to methods that don’t use whitening. Finally, we show that HoroPCA can be used to visualize hyperbolic data in two dimensions.

1 Introduction

Learning representations of data in hyperbolic spaces has recently attracted important interest in Machine Learning (ML) [25, 27] due to their ability to represent hierarchical data with high fidelity in low dimensions [28]. Many real-world datasets exhibit hierarchical structures, and hyperbolic embeddings have led to state-of-the-art results in applications such as question answering [30], node classification [5], link prediction [1, 7] and word embeddings [32]. These developments motivate the need for algorithms that operate in hyperbolic spaces such as nearest neighbor search [21, 35], hierarchical clustering [24, 6], or dimensionality reduction which is the focus of this work.

Euclidean Principal Component Analysis (PCA) is a fundamental dimensionality reduction technique which seeks directions that best explain the original data. PCA is an important primitive in data analysis and has many important uses such as (i) dimensionality reduction (e.g. for memory efficiency), (ii) data whitening and pre-processing for downstream tasks, and (iii) data visualization.

Here, we seek a generalization of PCA to hyperbolic geometry. Given a core notion of directions, PCA involves the following ingredients:

1.

A nested sequence of affine subspaces (flags) spanned by a set of directions.
2.

A projection method which maps points to these subspaces while preserving information (e.g. dot-product) along each direction.
3.

A variance objective to help choose the best directions.

Refer to caption — (a) Euclidean projections.

The PCA algorithm is then defined as combining these primitives: given a dataset, it chooses directions that maximize the variance of projections onto a subspace spanned by those directions, so that the resulting sequence of directions optimally explains the data. Crucially, the algorithm only depends on the directions of the affine subspaces and not their locations in space (Fig. 1(a)). Thus, in practice we can assume that they all go through the origin (and hence become linear subspaces), which greatly simplifies computations.

Generalizing PCA to manifolds is a challenging problem that has been studied for decades, starting with Principal Geodesic Analysis (PGA) [10] which parameterizes subspaces using tangent vectors at the mean of the data, and maximizes distances from projections to the mean to find optimal directions. More recently, the Barycentric Subspace Analysis (BSA) method [26] was introduced. It finds a more general parameterization of nested sequences of submanifolds by minimizing the unexplained variance. However, both PGA and BSA map points onto submanifolds using closest-point or geodesic projections, which do not attempt to preserve information along principal directions; for example, they cannot isometrically preserve hyperbolic distances between any points and shrink all path lengths by an exponential factor (3.5).

Fundamentally, all previous methods only look for subspaces rather than directions that explain the data, and can perhaps be better understood as principal subspace analysis rather than principal component analysis. Like with PCA, they assume that all optimal subspaces go through a chosen base point, but unlike in the Euclidean setting, this assumption is now unjustified: “translating” the submanifolds does not preserve distances between projections (Fig. 1(b)). Furthermore, the dependence on the base point is sensitive: as noted above, the shrink factor of the projection depends exponentially on the distances between the subspaces and the data. Thus, having to choose a base point increases the number of necessary parameters and reduces stability.

Here, we propose HoroPCA, a dimensionality reduction method for data defined in hyperbolic spaces which better preserves the properties of Euclidean PCA. We show how to interpret directions using points at infinity (or ideal points), which then allows us to generalize core properties of PCA:

1.

To generalize the notion of affine subspace, we propose parameterizing geodesic subspaces as the sets spanned by these ideal points. This yields multiple viable nested subspaces (flags) (Section 3.1).
2.

To maximally preserve information in the original data, we propose a new projection method that uses horospheres, a generalization of complementary directions for hyperbolic space. In contrast with geodesic projections, these projections exactly preserve information – specifically, distances to ideal points – along each direction. Consequently, they preserve distances between points much better than geodesic projections (Section 3.2).
3.

Finally, we introduce a simple generalization of explained variance that is a function of distances only and can be computed in hyperbolic space (Section 3.3).

Combining these notions, we propose an algorithm that seeks a sequence of principal components that best explain variations in hyperbolic data. We show that this formulation retains the location-independence property of PCA: translating target submanifolds along orthogonal directions (horospheres) preserves projected distances (Fig. 1(c)). In particular, the algorithm’s objective depends only on the directions and not locations of the submanifolds (Section 4).

We empirically validate HoroPCA on real datasets and for three standard PCA applications. First, (i) we show that it yields much lower distortion and higher explained variance than existing methods, reducing average distortion by up to 77%. Second, (ii) we validate that it can be used for data pre-processing, improving downstream classification by up to 3.8% in Average Precision score compared to methods that don’t use whitening. Finally, (iii) we show that the low-dimensional representations learned by HoroPCA can be visualized to qualitatively interpret hyperbolic data.

2 Background

We first review some basic notions from hyperbolic geometry; a more in-depth treatment is available in standard texts [22]. We discuss the generalization of coordinates and directions in hyperbolic space and then review geodesic projections. We finally describe generalizations of the notion of mean and variance to non-Euclidean spaces.

2.1 The Poincaré Model of Hyperbolic Space

Hyperbolic geometry is a Riemannian geometry with constant negative curvature $-1$ , where curvature measures deviation from flat Euclidean geometry. For easier visualizations, we work with the $d$ -dimensional Poincaré model of hyperbolic space: $\mathbb{H}^{d}=\{{x}\in\mathbb{R}^{d}:\|x\|<1\}$ , where $\|\cdot\|$ is the Euclidean norm. In this model, the Riemannian distance can be computed in cartesian coordinates by:

d_{\mathbb{H}}(x,y)=\operatorname{cosh^{-1}}\left(1+2\frac{\|x-y\|^{2}}{(1-\|x\|^{2})(1-\|y\|^{2})}\right).

(1)

Geodesics

Shortest paths in hyperbolic space are called geodesics. In the Poincaré model, they are represented by straight segments going through the origin and circular arcs perpendicular to the boundary of the unit ball (Fig. 2).

Geodesic submanifolds

A submanifold $M\subset\mathbb{H}^{d}$ is called (totally) geodesic if for every $x,y\in M$ , the geodesic line connecting $x$ and $y$ belongs to $M$ . This generalizes the notion of affine subspaces in Euclidean spaces. In the Poincaré model, geodesic submanifolds are represented by linear subspaces going through the origin and spherical caps perpendicular to the boundary of the unit ball.

	Euclidean	Hyperbolic
Component	Unit vector $w$	Ideal point $p$
Coordinate	Dot product $x\cdot w$	Busemann func. $B_{p}(x)$

Table 1: Analogies of components and their corresponding coordinates, in both Euclidean and hyperbolic space.

2.2 Directions in Hyperbolic space

The notions of directions, and coordinates in a given direction can be generalized to hyperbolic spaces as follows.

Ideal points

As with parallel rays in Euclidean spaces, geodesic rays in $\mathbb{H}^{d}$ that stay close to each other can be viewed as sharing a common endpoint at infinity, also called an ideal point. Intuitively, ideal points represent directions along which points in $\mathbb{H}^{d}$ can move toward infinity. The set of ideal points $\mathbb{S}_{\infty}^{d-1}$ , called the boundary at infinity of $\mathbb{H}^{d}$ , is represented by the unit sphere $\mathbb{S}_{\infty}^{d-1}=\{\|x\|=1\}$ in the Poincaré model. We abuse notations and say that a geodesic submanifold $M\subset\mathbb{H}^{d}$ contains an ideal point $p$ if the boundary of $M$ in $\mathbb{S}_{\infty}^{d-1}$ contains $p$ .

Busemann coordinates

In Euclidean spaces, each direction can be represented by a unit vector $w$ . The coordinate of a point $x$ in the direction of $w$ is simply the dot product $w\cdot x$ . In hyperbolic geometry, directions can be represented by ideal points but dot products are not well-defined. Instead, we take a ray-based perspective: note that in Euclidean spaces, if we shoot a ray in the direction of $w$ from the origin, the coordinate $w\cdot x$ can be viewed as the normalized distance to infinity in the direction of that ray. In other words, as a point $y=tw$ , ( $t>0$ ) moves toward infinity in the direction of $w$ :

w\cdot x=\lim_{t\to\infty}\left(d(0,tw)-d(x,tw)\right).

This approach generalizes to other geometries: given a unit-speed geodesic ray $\gamma(t)$ , the Busemann function $B_{\gamma}(x)$ of $\gamma$ is defined as:¹¹1Note that compared to the above formula, the sign convention is flipped due to historical reasons.

B_{\gamma}(x)=\lim_{t\to\infty}\left(d(x,\gamma(t))-t\right).

Up to an additive constant, this function only depends on the endpoint at infinity of the geodesic ray, and not the starting point $\gamma(0)$ . Thus, given an ideal point $p$ , we define the Busemann function $B_{p}(x)$ of $p$ to be the Busemann function of the geodesic ray that starts from the origin of the unit ball model and has endpoint $p$ . Intuitively, it represents the coordinates of $x$ in the direction of $p$ . In the Poincaré model, there is a closed formula:

B_{p}(x)=\ln\frac{\|p-x\|^{2}}{1-\|x\|^{2}}.

Horospheres

The level sets of Busemann functions $B_{p}(x)$ are called horospheres centered at $p$ . In this sense, they resemble spheres, which are level sets of distance functions. However, intrinsically as Riemannian manifolds, horospheres have curvature zero and thus also exhibit many properties of planes in Euclidean spaces.

Every geodesic with endpoint $p$ is orthogonal to every horosphere centered at $p$ . Given two horospheres with the same center, every orthogonal geodesic segment connecting them has the same length. In this sense, concentric horospheres resemble parallel planes in Euclidean spaces. In the Poincaré model, horospheres are Euclidean spheres that touch the boundary sphere $\mathbb{S}_{\infty}^{d-1}$ at their ideal centers (Fig. 2). Given an ideal point $p$ and a point $x$ in $\mathbb{H}^{d}$ , there is a unique horosphere $S(p,x)$ passing through $x$ and centered at $p$ .

2.3 Geodesic Projections

PCA uses orthogonal projections to project data onto subspaces. Orthogonal projections are usually generalized to other geometries as closest-point projections. Given a target submanifold $M$ , each point $x$ in the ambient space is mapped to the closest-point to it in $M$ :

\pi^{\operatorname{G}}_{M}(x)=\operatorname*{argmin}_{y\in M}d_{M}(x,y).

One could view $\pi^{\operatorname{G}}_{M}(\cdot)$ as the map that pushes each point $x$ along an orthogonal geodesic until it hits $M$ . For this reason, it is also called geodesic projection. In the Poincaré model, these can be computed in closed-form (see Appendix C).

2.4 Manifold Statistics

PCA relies on data statistics which do not generalize easily to hyperbolic geometry. One approach to generalize the arithmetic mean is to notice that it is the minimizer of the sum of squared distances to the inputs. Motivated by this, the Fréchet mean [11] of a set of points $S$ in a Riemannian manifold $(M,d_{M})$ is defined as:

\mu_{M}(S)\coloneqq\operatorname*{argmin}_{y\in M}\sum_{x\in S}d_{M}(x,y)^{2}.

This definition only depends on the intrinsic distance of the manifold. For hyperbolic spaces, since squared distance functions are convex, $\mu(S)$ always exists and is unique.²²2For more general geometries, existence and uniqueness hold if the data is well-localized [20]. Analogously, the Fréchet variance is defined as:

\sigma^{2}_{M}(S)\coloneqq\frac{1}{|S|}\sum_{x\in S}^{N}d_{M}(x,\mu(S))^{2}.

(2)

We refer to [16] for a discussion on different intrinsic statistics in non-Euclidean spaces, and a study of their asymptotic properties.

3 Generalizing PCA to the Hyperbolic Space

We now describe our approach to generalize PCA to hyperbolic spaces. The starting point of HoroPCA is to pick $1\leq K\leq d$ ideal points $p_{1},\dots,p_{K}\in\mathbb{S}_{\infty}^{d-1}$ to represent $K$ directions in hyperbolic spaces (Section 2.2). Then, we generalize the core concepts of Euclidean PCA. In Section 3.1, we show how to generalize flags. In Section 3.2, we show how to project points onto the lower-dimensional submanifold spanned by a given set of directions, while preserving information along each direction. In Section 3.3, we introducing a variance objective to optimize and show that it is a function of the directions only.

3.1 Hyperbolic Flags

In Euclidean spaces, one can take the linear spans of more and more components to define a nested sequence of linear subspaces, called a flag. To generalize this to hyperbolic spaces, we first need to adapt the notion of linear/affine spans. Recall that geodesic submanifolds are generalizations of affine subspaces in Euclidean spaces.

Definition 3.1.

Given a set of points $S$ (that could be inside $\mathbb{H}^{d}$ or on the boundary sphere $\mathbb{S}_{\infty}^{d-1}$ ), the smallest geodesic submanifold of $\mathbb{H}^{d}$ that contains $S$ is called the geodesic hull of $S$ and denoted by $\operatorname{GH}(S)$ .

Thus, given $K$ ideal points $p_{1},p_{2},\dots,p_{K}$ and a base point $b\in\mathbb{H}^{d}$ , we can define a nested sequence of geodesic submanifolds $\operatorname{GH}(b,p_{1})\subset\operatorname{GH}(b,p_{1},p_{2})\subset\dots\subset\operatorname{GH}(b,p_{1},\dots,p_{K})$ . This will be our notion of flags.

Remark 3.2.

The base point $b$ is only needed here for technical reasons, just like an origin $\mathbf{o}$ is needed to define linear spans in Euclidean spaces. We will see next that it does not affect the projection results or objectives (Theorem 4.1).

Remark 3.3.

We assume that none of $b,p_{1},\dots,p_{K}$ are in the geodesic hull of the other $K$ points. This is analogous to being linearly independent in Euclidean spaces.

3.2 Projections via Horospheres

In Euclidean PCA, points are projected to the subspaces spanned by the given directions in a way that preserves coordinates in those directions. We seek a projection method in hyperbolic spaces with a similar property.

Recall that coordinates are generalized by Busemann functions (Table 1), and that horospheres are level sets of Busemann functions. Thus, we propose a projection that preserves coordinates by moving points along horospheres. It turns out that this projection method also preserves distances better than the traditional geodesic projection.

As a toy example, we first show how the projection is defined in the $K=1$ case (i.e. projecting onto a geodesic) and why it tends to preserve distances well. We will then show how to use $K\geq 1$ ideal points simultaneously.

3.2.1 Projecting onto $K=1$ Directions

For $K=1$ , we have one ideal point $p$ and base point $b$ , and the geodesic hull $\operatorname{GH}(b,p)$ is just a geodesic $\gamma$ . Our goal is to map every $x\in\mathbb{H}^{d}$ to a point $\pi^{\operatorname{H}}_{b,p}(x)$ on $\gamma$ that has the same Busemann coordinate in the direction of $p$ :

B_{p}(x)=B_{p}(\pi^{\operatorname{H}}_{b,p}(x)).

Since level sets of $B_{p}(x)$ are horospheres centered at $p$ , the above equation simply says that $\pi^{\operatorname{H}}_{b,p}(x)$ belongs to the horosphere $S(p,x)$ centered at $p$ and passing through $x$ . Thus, we define:

\pi^{\operatorname{H}}_{b,p}(x)\coloneqq\gamma\cap S(p,x).

(3)

Another important property that $\pi^{\operatorname{H}}_{b,p}(\cdot)$ shares with orthogonal projections in Euclidean spaces is that it preserves distances along a direction – lengths of geodesic segments that point to $p$ are preserved after projection (Fig. 3):

Proposition 3.4.

For any $x\in\mathbb{H}^{d}$ , if $y\in\operatorname{GH}(x,p)$ then:

d_{\mathbb{H}}(\pi^{\operatorname{H}}_{b,p}(x),\pi^{\operatorname{H}}_{b,p}(y))=d_{\mathbb{H}}(x,y).

Proof.

This follows from the remark in Section 2.2 about horospheres: every geodesic going through $p$ is orthogonal to every horosphere centered at $p$ , and every orthogonal geodesic segment connecting concentric horospheres has the same length (Fig. 2). In this case, the segments from $x$ to $y$ and from $\pi^{\operatorname{H}}_{b,p}(x)$ to $\pi^{\operatorname{H}}_{b,p}(y)$ are two such segments, connecting $S(p,x)$ and $S(p,y)$ . ∎

3.2.2 Projecting onto $K>1$ Directions

We now generalize the above construction to projections onto higher-dimensional submanifolds. We describe the main ideas here; Appendix A contains more details, including an illustration in the case $K=2$ (Fig. 5).

Fix a base point $b\in\mathbb{H}^{d}$ and $K>1$ ideal points $\{p_{1},\dots,p_{K}\}$ . We want to define a map from $\mathbb{H}^{d}$ to $M\coloneqq\operatorname{GH}(b,p_{1},\dots,p_{K})$ that preserves the Busemann coordinates in the directions of $p_{1},\dots,p_{K}$ , i.e.:

B_{p_{j}}(x)=B_{p_{j}}\left(\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}(x)\right)\text{ for every }j=1,\dots,K.

As before, the idea is to take the intersection with the horospheres centered at $p_{j}$ ’s and passing through $x$ :

	$\displaystyle\pi^{\text{H}}_{b,p_{1},\dots,p_{K}}:\mathbb{H}^{d}$	$\displaystyle\to M$
	$\displaystyle x$	$\displaystyle\mapsto M\cap S(p_{1},x)\cap\dots\cap S(p_{K},x).$

It turns out that this intersection generally consists of two points instead of one. When that happens, one of them will be strictly closer to the base point $b$ , and we define $\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}(x)$ to be that point.

As with 3.4, $\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}(\cdot)$ preserves distances along $K$ -dimensional manifolds (A.10). In contrast, geodesic projections in hyperbolic spaces never preserve distances (except between points already in the target):

Proposition 3.5.

Let $M\subset\mathbb{H}^{d}$ be a geodesic submanifold. Then every geodesic segment of distance at least $r$ from $M$ gets at least $\cosh(r)$ times shorter under the geodesic projection $\pi^{\operatorname{G}}_{M}(\cdot)$ to $M$ :

\operatorname{length}(\pi^{\operatorname{G}}_{M}(I))\leq\frac{1}{\cosh(r)}\operatorname{length}(I).

In particular, the shrink factor grows exponentially as the segment $I$ moves away from $M$ . ∎

The proof is in Appendix B.

Computation

Interestingly, horosphere projections can be computed without actually computing the horospheres. The key idea is that if we let $P=\operatorname{GH}(p_{1},\dots,p_{K})$ be the geodesic hull of the horospheres’ centers, then the intersection $S(p_{1},x)\cap\dots\cap S(p_{K},x)$ is simply the orbit of $x$ under the rotations around $P$ . (This is true for the same reason that spheres whose centers lie on the same axis must intersect along a circle around that axis). Thus, $\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}(\cdot)$ can be viewed as the map that rotates $x$ around until it hits $M$ . It follows that it can be computed by:

1.

Find the geodesic projection $c=\pi^{\operatorname{G}}_{P}(x)$ of $x$ onto $P$ .
2.

Find the geodesic $\alpha$ on $M$ that is orthogonal to $P$ at $c$ .
3.

Among the two points on $\alpha$ whose distance to $c$ equals $d_{\mathbb{H}}(x,c)$ , returns the one closer to $b$ .

The detailed computations and proof that this recovers horospherical projections are provided in Appendix A.

3.3 Intrinsic Variance Objective

In Euclidean PCA, directions are chosen to maximally preserve information from the original data. In particular, PCA chooses directions that maximize the Euclidean variance of projected data. To generalize this to hyperbolic geometry, we define an analog of variance that is intrinsic, i.e. dependent only on the distances between data points. As we will see in Section 4, having an intrinsic objective helps make our algorithm location (or base point) independent.

The usual notion of Euclidean variance is the squared sum of distances to the mean of the projected datapoints. Generalizing this is challenging because non-Euclidean spaces do not have a canonical choice of mean. Previous works have generalized variance either by using the unexplained variance or Fréchet variance. The former is the squared sum of residual distances to the projections, and thus avoids computing a mean. However, it is not intrinsic. The latter is intrinsic [10] but involves finding the Fréchet mean, which is not necessarily a canonical notion of mean and can only be computed by gradient descent.

Our approach uses the observation that in Euclidean space:

\sigma^{2}(S)=\frac{1}{n}\sum_{x\in S}\left\|x-\mu(S)\right\|^{2}=\frac{1}{n^{2}}\sum_{x,y\in S}\|x-y\|^{2}.

Thus, we propose the following generalization of variance:

\sigma_{\mathbb{H}}^{2}(S)=\frac{1}{n^{2}}\sum_{x,y\in S}d_{\mathbb{H}}(x,y)^{2}.

(4)

This function agrees with the usual variance in Euclidean space, while being a function of distances only. Thus it is well defined in non-Euclidean space, is easily computed, and, as we will show next, has the desired invariance due to isometry properties of horospherical projections.

4 HoroPCA

Section 3 formulated several simple primitives – including directions, flags, projections, and variance – in a way that is generalizable to hyperbolic geometry. We now revisit standard PCA, showing how it has a simple definition that combines these primitives using optimization. This directly leads to the full HoroPCA algorithm by simply using the hyperbolic analogs of these primitives.

Euclidean PCA

Given a dataset $S$ and a target dimension $K$ , Euclidean PCA greedily finds a sequence of principal components $p_{1},\dots,p_{K}$ that maximizes the variance of orthogonal projections $\pi^{\operatorname{E}}_{\mathbf{o},p_{1},\dots,p_{k}}(\cdot)$ onto the linear ³³3Here $\mathbf{o}$ denotes the origin, and $\pi^{\operatorname{E}}_{\mathbf{o},p_{1},\dots,p_{k}}(\cdot)$ denotes the projection onto the affine span of $\{\mathbf{o},p_{1},\dots,p_{k}\}$ , which is equivalent to the linear span of $\{p_{1},\dots,p_{k}\}$ . subspaces spanned by these components:

	$\displaystyle p_{1}$	$\displaystyle=\underset{\|\|p\|\|=1}{\mathrm{argmax}}\ \sigma^{2}(\pi^{\operatorname{E}}_{\mathbf{o},p}(S))$
	$\displaystyle\text{and}\ p_{k+1}$	$\displaystyle=\underset{\|\|p\|\|=1}{\mathrm{argmax}}\ \sigma^{2}(\pi^{\operatorname{E}}_{\mathbf{o},p_{1},\dots,p_{k},p}(S)).$

Thus, for each $1\leq k\leq K$ , $\{p_{1},\dots,p_{k}\}$ is the optimal set of directions of dimension $k$ .

HoroPCA

Because we have generalized the notions of flag, projection, and variance to hyperbolic geometry, the HoroPCA algorithm can be defined in the same fashion. Given a dataset $S$ in $\mathbb{H}^{d}$ and a base point $b\in\mathbb{H}^{d}$ , we seek a sequence of $K$ directions that maximizes the variance of horosphere-projected data:

\begin{split}p_{1}&=\operatorname*{argmax}_{p\in\mathbb{S}_{\infty}^{d-1}}\sigma^{2}_{\mathbb{H}}(\pi^{\operatorname{H}}_{b,p}(S))\\ \text{and}\ p_{k+1}&=\operatorname*{argmax}_{p\in\mathbb{S}_{\infty}^{d-1}}\sigma^{2}_{\mathbb{H}}(\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{k},p}(S)).\end{split}

(5)

Base point independence

Finally, we show that algorithm (5) always returns the same results regardless of the choice of a base point $b\in\mathbb{H}^{d}$ . Since our variance objective only depends on the distances between projected data points, it suffices to show that these distances do not depend on $b$ .

Theorem 4.1.

For any $b,b^{\prime}$ and any $x,y\in\mathbb{H}^{d}$ , the two projected distances $d_{\mathbb{H}}(\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}(x),\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}(y))$ and $d_{\mathbb{H}}(\pi^{\operatorname{H}}_{b^{\prime},p_{1},\dots,p_{K}}(x),\pi^{\operatorname{H}}_{b^{\prime},p_{1},\dots,p_{K}}(y))$ are equal.

The proof is included in Appendix A. Thus, HoroPCA retains the location-independence property of Euclidean PCA: only the directions of target subspaces matter; their exact locations do not (Fig. 1). Therefore, just like in the Euclidean setting, we can assume without loss of generality that $b$ is the origin $\mathbf{o}$ of the Poincaré model. This:

1.

alleviates the need to use $d$ extra parameters to search for an appropriate base point, and
2.

simplifies computations, since in the Poincaré model, geodesics submanifolds that go through the origin are simply linear subspaces, which are easier to deal with.

After computing the principal directions which span the target $M=\operatorname{GH}(\mathbf{o},p_{1},\dots,p_{K})$ , the reduced dimensionality data can be found by applying an Euclidean rotation that sends $M$ to $\mathbb{H}^{K}$ , which also preserves hyperbolic distances.

5 Experiments

	Balanced Tree		Phylo Tree		Diseases		CS Ph.D.
	distortion ( $\downarrow$ )	variance ( $\uparrow$ )	distortion ( $\downarrow$ )	variance ( $\uparrow$ )	distortion ( $\downarrow$ )	variance ( $\uparrow$ )	distortion ( $\downarrow$ )	variance ( $\uparrow$ )
PCA	0.84	0.34	0.94	0.40	0.90	0.26	0.84	1.68
tPCA	0.70	1.16	0.63	14.34	0.63	3.92	0.56	11.09
PGA	0.63 $\pm$ 0.07	2.11 $\pm$ 0.47	0.64 $\pm$ 0.01	15.29 $\pm$ 0.51	0.66 $\pm$ 0.02	3.16 $\pm$ 0.39	0.73 $\pm$ 0.02	6.14 $\pm$ 0.60
PGA-Noise	0.87 $\pm$ 0.08	0.29 $\pm$ 0.30	0.64 $\pm$ 0.02	15.08 $\pm$ 0.77	0.88 $\pm$ 0.04	0.53 $\pm$ 0.19	0.79 $\pm$ 0.03	4.58 $\pm$ 0.64
BSA	0.50 $\pm$ 0.00	3.02 $\pm$ 0.01	0.61 $\pm$ 0.03	18.60 $\pm$ 1.16	0.52 $\pm$ 0.02	5.95 $\pm$ 0.25	0.70 $\pm$ 0.01	8.15 $\pm$ 0.96
BSA-Noise	0.74 $\pm$ 0.12	1.06 $\pm$ 0.67	0.68 $\pm$ 0.02	13.71 $\pm$ 0.72	0.80 $\pm$ 0.11	1.62 $\pm$ 1.30	0.79 $\pm$ 0.02	4.41 $\pm$ 0.59
hAE	0.26 $\pm$ 0.00	6.91 $\pm$ 0.00	0.32 $\pm$ 0.04	45.87 $\pm$ 3.52	0.18 $\pm$ 0.00	14.23 $\pm$ 0.06	0.37 $\pm$ 0.02	22.12 $\pm$ 2.47
hMDS	0.22	7.54	0.74	40.51	0.21	15.05	0.83	19.93
HoroPCA	0.19 $\pm$ 0.00	7.15 $\pm$ 0.00	0.13 $\pm$ 0.01	69.16 $\pm$ 1.96	0.15 $\pm$ 0.01	15.46 $\pm$ 0.19	0.16 $\pm$ 0.02	36.79 $\pm$ 0.70

Table 2: Dimensionality reduction results on 10-dimensional hyperbolic embeddings reduced to two dimensions. Results are averaged over 5 runs for non-deterministic methods. Best in bold and second best underlined.

We now validate the empirical benefits of HoroPCA on three PCA uses. First, for dimensionality reduction, HoroPCA preserves information (distances and variance) better than previous methods which are sensitive to base point choices and distort distances more (Section 5.2). Next, we validate that our notion of hyperbolic coordinates captures variation in the data and can be used for whitening in classification tasks (Section 5.3). Finally, we visualize the representations learned by HoroPCA in two dimensions (Section 5.4).

5.1 Experimental Setup

Baselines

We compare HoroPCA to several dimensionality reduction methods, including: (1) Euclidean PCA, which should perform poorly on hyperbolic data, (2) Exact PGA, (3) Tangent PCA (tPCA), which approximates PGA by moving the data in the tangent space of the Fréchet mean and then solves Euclidean PCA, (4) BSA, (5) Hyperbolic Multi-dimensional Scaling (hMDS) [27], which takes a distance matrix as input and recovers a configuration of points that best approximates these distances, (6) Hyperbolic autoencoder (hAE) trained with gradient descent [12, 15]. To demonstrate their dependence on base points, we also include two baselines that perturb the base point in PGA and BSA. We open-source our implementation⁴⁴4https://github.com/HazyResearch/HoroPCA and refer to Appendix C for implementation details on how we implemented all baselines and HoroPCA.

Datasets

For dimensionality reduction experiments, we consider standard hierarchical datasets previously used to evaluate the benefits of hyperbolic embeddings. More specifically, we use the datasets in [27] including a fully balanced tree, a phylogenetic tree, a biological graph comprising of diseases’ relationships and a graph of Computer Science (CS) Ph.D. advisor-advisee relationships. These datasets have respectively 40, 344, 516 and 1025 nodes, and we use the code from [13] to embed them in the Poincaré ball. For data whitening experiments, we reproduce the experimental setup from [8] and use the Polbooks, Football and Polblogs datasets which have 105, 115 and 1224 nodes each. These real-world networks are embedded in two-dimensions using Chamberlain et al. [4]’s embedding method.

Evaluation metrics

To measure distance-preservation after projection, we use average distortion. If $\pi(\cdot)$ denotes a mapping from high- to low-dimensional representations, the average distortion of a dataset $S$ is computed as:

\frac{1}{\binom{|S|}{2}}\sum_{x\neq y\in S}\frac{|d_{\mathbb{H}}(\pi(x),\pi(y))-d_{\mathbb{H}}(x,y)|}{d_{\mathbb{H}}(x,y)}.

We also measure the Fréchet variance in Eq. 2, which is the analogue of the objective that Euclidean PCA optimizes⁵⁵5All mentioned PCA methods, including HoroPCA, optimize for some forms of variance but not Fréchet variance or distortion.. Note that the mean in Eq. 2 cannot be computed in closed-form and we therefore compute it with gradient-descent.

5.2 Dimensionality Reduction

We report metrics for the reduction of 10-dimensional embeddings to two dimensions in Table 2, and refer to Appendix C for additional results, such as more component and dimension configurations. All results suggest that HoroPCA better preserves information contained in the high-dimensional representations.

On distance preservation, HoroPCA outperforms all methods with significant improvements on larger datasets. This supports our theoretical result that horospherical projections better preserve distances than geodesic projections. Furthermore, HoroPCA also outperforms existing methods on the explained Fréchet variance metric on all but one dataset. This suggests that our distance-based formulation of the variance (Eq. 4) effectively captures variations in the data. We also note that as expected, both PGA and BSA are sensitive to base point choices: adding Gaussian noise to the base point leads to significant drops in performance. In contrast, HoroPCA is by construction base-point independent.

5.3 Hyperbolic Data Whitening

An important use of PCA is for data whitening, as it allows practitioners to remove noise and decorrelate the data, which can improve downstream tasks such as regression or classification. Recall that standard PCA data whitening consists of (i) finding principal directions that explain the data, (ii) calculating the coordinates of each data point along these directions, and (iii) normalizing the coordinates for each direction (to have zero mean and unit variance).

Because of the close analogy between HoroPCA and Euclidean PCA, these steps can easily map to the hyperbolic case, where we (i) use HoroPCA to find principal directions (ideal points), (ii) calculate the Busemann coordinates along these directions, and (iii) normalize them as Euclidean coordinates. Note that this yields Euclidean representations, which allow leveraging powerful tools developed specifically for learning on Euclidean data.

We evaluate the benefit of this whitening step on a simple classification task. We compare to directly classifying the data with Euclidean Support Vector Machine (eSVM) or its hyperbolic counterpart (hSVM), and also to whitening with tPCA. Note that most baselines in Section 5.1 are incompatible with data whitening: hMDS does not learn a transformation that can be applied to unseen test data, while methods like PGA and BSA do not naturally return Euclidean coordinates for us to normalize. To obtain another baseline, we use a logarithmic map to extract Euclidean coordinates from PGA.

	Polbooks	Football	Polblogs
eSVM	69.9 $\pm$ 1.2	20.7 $\pm$ 3.0	92.3 $\pm$ 1.5
hSVM	68.3 $\pm$ 0.6	20.9 $\pm$ 2.5	92.2 $\pm$ 1.6
tPCA+eSVM	68.5 $\pm$ 0.9	21.2 $\pm$ 2.2	92.4 $\pm$ 1.5
PGA+eSVM	64.4 $\pm$ 4.1	21.7 $\pm$ 2.2	82.3 $\pm$ 1.2
HoroPCA+eSVM	72.2 $\pm$ 2.8	25.0 $\pm$ 1.0	92.8 $\pm$ 0.9

Table 3: Data whitening experiments. We report classification accuracy averaged over 5 embedding configurations. Best in bold and second best underlined.

We reproduce the experimental setup from [8] who split the datasets in 50% train and 50% test sets, run classification on 2-dimensional embeddings and average results over 5 different embedding configurations as was done in the original paper (Table 3). ⁶⁶6Note that the results slightly differ from [8] which could be because of different implementations or data splits. HoroPCA whitening improves downstream classification on all datasets compared to eSVM and hSVM or tPCA and PGA whitening. This suggests that HoroPCA can be leveraged for hyperbolic data whitening. Further, this confirms that Busemann coordinates do capture variations in the original data.

5.4 Visualizations

When learning embeddings for ML applications (e.g. classification), increasing the dimensionality can significantly improve the embeddings’ quality. To effectively work with these higher-dimensional embeddings, it is useful to visualize their structure and organization, which often requires reducing their representations to two or three dimensions. Here, we consider embeddings of the mammals subtree of the Wordnet noun hierarchy learned with the algorithm from [25]. We reduce embeddings to two dimensions using PGA and HoroPCA and show the results in Fig. 4. We also include more visualizations for PCA and BSA in Fig. 8 in the Appendix. As we can see, the reduced representations obtained with HoroPCA yield better visualizations. For instance, we can see some hierarchical patterns such as “feline hypernym of cat” or “cat hypernym of burmese cat”. These patterns are harder to visualize for other methods, since these do not preserve distances as well as HoroPCA, e.g. PGA has 0.534 average distortion on this dataset compared to 0.078 for HoroPCA.

6 Related Work

PCA methods in Riemannian manifolds

We first review some approaches to extending PCA to general Riemannian geometries, of which hyperbolic geometry is a special case. For a more detailed discussion, see Pennec [26]. The simplest such approach is tangent PCA (tPCA), which maps the data to the tangent space at the Fréchet mean $\mu$ using the logarithm map, then applies Euclidean PCA. A similar approach, Principal Geodesic Analysis (PGA) [10], seeks geodesic subspaces at $\mu$ that minimize the sum of squared Riemannian distances to the data. Compared to tPCA, PGA searches through the same subspaces but uses a more natural loss function.

Both PGA and tPCA project on submanifolds that go through the Fréchet mean. When the data is not well-centered, this may be sub-optimal, and Geodesic PCA (GPCA) was proposed to alleviate this issue [17, 18]. GPCA first finds a geodesic $\gamma$ that best fits the data, then finds other orthogonal geodesics that go through some common point $b$ on $\gamma$ . In other words, GPCA removes the constraint of PGA that $b$ is the Fréchet mean. Extensions of GPCA have been proposed such as probabilistic methods [36] and Horizontal Component Analysis [29].

Pennec [26] proposes a more symmetric approach. Instead of using the exponential map at a base point, it parameterizes $K$ -dimensional subspaces as the barycenter loci of $K+1$ points. Nested sequences of subspaces (flags) can be formed by simply adding more points. In hyperbolic geometry, this construction coincides with the one based on geodesic hulls that we use in Section 3, except that it applies to points inside $\mathbb{H}^{d}$ instead of ideal points and thus needs more parameters to parameterize a flag (see Remark A.8).

By considering a more general type of submanifolds, Hauberg [14] gives another way to avoid the sensitive dependence on Fréchet mean. However, its bigger search space also makes the method computationally expensive, especially when the target dimension $K$ is bigger than $1$ .

In contrast with all methods so far, HoroPCA relies on horospherical projections instead of geodesic projections. This yields a generalization of PCA that depends only on the directions and not specific locations of the subspaces.

Dimension reduction in hyperbolic geometry

We now review some dimension reduction methods proposed specifically for hyperbolic geometry. Cvetkovski & Crovella [9] and Sala et al. [27] are examples of hyperbolic multidimensional scaling methods, which seek configurations of points in lower-dimensional hyperbolic spaces whose pairwise distances best approximate a given dissimilarity matrix. Unlike HoroPCA, they do not learn a projection that can be applied to unseen data.

Tran & Vu [33] constructs a map to lower-dimensional hyperbolic spaces whose preimages of compact sets are compact. Unlike most methods, it is data-agnostic and does not optimize any objective.

Benjamini & Makarychev [2] adapt the Euclidean Johnson–Lindenstrauss transform to hyperbolic geometry and obtain a distortion bound when the dataset size is not too big compared to the target dimension. They do not seek an analog of Euclidean directions or projections, but nevertheless implicitly use a projection method based on pushing points along horocycles, which shares many properties with our horospherical projections. In fact, the latter converges to the former as the ideal points get closer to each other.

From hyperbolic to Euclidean

Liu et al. [23] use “distances to centroids” to compute Euclidean representations of hyperbolic data. The Busemann functions we use bear resemblances to these centroid-based functions but are better analogs of coordinates along given directions, which is a central concept in PCA, and have better regularity properties [3]. Recent works have also used Busemann functions for hyperbolic prototype learning [19, 34]. These works do not define projections to lower-dimensional hyperbolic spaces. In contrast, HoroPCA naturally returns both hyperbolic representations (via horospherical projections) and Euclidean representations (via Busemann coordinates). This allows leveraging techniques in both settings.

7 Conclusion

We proposed HoroPCA, a method to generalize PCA to hyperbolic spaces. In contrast with previous PCA generalizations, HoroPCA preserves the core location-independence PCA property. Empirically, HoroPCA significantly outperforms previous methods on the reduction of hyperbolic data. Future extensions of this work include deriving a closed-form solution, analyzing the stability properties of HoroPCA, or using the concepts introduced in this work to derive efficient nearest neighbor search algorithms or neural network operations.

Acknowledgements

We gratefully acknowledge the support of NIH under No. U54EB020405 (Mobilize), NSF under Nos. CCF1763315 (Beyond Sparsity), CCF1563078 (Volume to Velocity), and 1937301 (RTML); ONR under No. N000141712266 (Unifying Weak Supervision); the Moore Foundation, NXP, Xilinx, LETI-CEA, Intel, IBM, Microsoft, NEC, Toshiba, TSMC, ARM, Hitachi, BASF, Accenture, Ericsson, Qualcomm, Analog Devices, the Okawa Foundation, American Family Insurance, Google Cloud, Swiss Re, Total, the HAI-AWS Cloud Credits for Research program, the Stanford Data Science Initiative (SDSI), and members of the Stanford DAWN project: Facebook, Google, and VMWare. The Mobilize Center is a Biomedical Technology Resource Center, funded by the NIH National Institute of Biomedical Imaging and Bioengineering through Grant P41EB027060. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views, policies, or endorsements, either expressed or implied, of NIH, ONR, or the U.S. Government.

References

Balazevic et al. [2019] Balazevic, I., Allen, C., and Hospedales, T. Multi-relational poincaré graph embeddings. In Advances in Neural Information Processing Systems, pp. 4463–4473, 2019.
Benjamini & Makarychev [2009] Benjamini, I. and Makarychev, Y. Dimension reduction for hyperbolic space. Proceedings of the American Mathematical Society, pp. 695–698, 2009.
Busemann [1955] Busemann, H. The geometry of geodesics. Academic Press Inc., New York, N. Y., 1955.
Chamberlain et al. [2017] Chamberlain, B., Clough, J., and Deisenroth, M. Neural embeddings of graphs in hyperbolic space. In CoRR. MLG Workshop 2017, 2017.
Chami et al. [2019] Chami, I., Ying, Z., Ré, C., and Leskovec, J. Hyperbolic graph convolutional neural networks. In Advances in neural information processing systems, pp. 4868–4879, 2019.
Chami et al. [2020a] Chami, I., Gu, A., Chatziafratis, V., and Ré, C. From trees to continuous embeddings and back: Hyperbolic hierarchical clustering. Advances in Neural Information Processing Systems, 33, 2020a.
Chami et al. [2020b] Chami, I., Wolf, A., Juan, D.-C., Sala, F., Ravi, S., and Ré, C. Low-dimensional hyperbolic knowledge graph embeddings. In Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, pp. 6901–6914, 2020b.
Cho et al. [2019] Cho, H., DeMeo, B., Peng, J., and Berger, B. Large-margin classification in hyperbolic space. In The 22nd International Conference on Artificial Intelligence and Statistics, pp. 1832–1840. PMLR, 2019.
Cvetkovski & Crovella [2011] Cvetkovski, A. and Crovella, M. Multidimensional scaling in the poincaré disk. arXiv preprint arXiv:1105.5332, 2011.
Fletcher et al. [2004] Fletcher, P. T., Conglin Lu, Pizer, S. M., and Sarang Joshi. Principal geodesic analysis for the study of nonlinear statistics of shape. IEEE Transactions on Medical Imaging, 23(8):995–1005, 2004.
Fréchet [1948] Fréchet, M. Les éléments aléatoires de nature quelconque dans un espace distancié. In Annales de l’institut Henri Poincaré, volume 10, pp. 215–310, 1948.
Ganea et al. [2018] Ganea, O.-E., Bécigneul, G., and Hofmann, T. Hyperbolic neural networks. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pp. 5350–5360, 2018.
Gu et al. [2018] Gu, A., Sala, F., Gunel, B., and Ré, C. Learning mixed-curvature representations in product spaces. In International Conference on Learning Representations, 2018.
Hauberg [2016] Hauberg, S. Principal curves on riemannian manifolds. IEEE Transactions on Pattern Analysis and Machine Intelligence, 38(9):1915–1921, 2016. doi: 10.1109/TPAMI.2015.2496166.
Hinton & Salakhutdinov [2006] Hinton, G. E. and Salakhutdinov, R. R. Reducing the dimensionality of data with neural networks. science, 313(5786):504–507, 2006.
Huckemann & Eltzner [2020] Huckemann, S. and Eltzner, B. Statistical methods generalizing principal component analysis to non-euclidean spaces. In Handbook of Variational Methods for Nonlinear Geometric Data, pp. 317–338. Springer, 2020.
Huckemann & Ziezold [2006] Huckemann, S. and Ziezold, H. Principal component analysis for riemannian manifolds, with an application to triangular shape spaces. Advances in Applied Probability, 38(2):299–319, 2006.
Huckemann et al. [2010] Huckemann, S., Hotz, T., and Munk, A. Intrinsic shape analysis: Geodesic pca for riemannian manifolds modulo isometric lie group actions. Statistica Sinica, pp. 1–58, 2010.
Keller-Ressel [2020] Keller-Ressel, M. A theory of hyperbolic prototype learning. arXiv preprint arXiv:2010.07744, 2020.
Kendall [1990] Kendall, W. S. Probability, convexity, and harmonic maps with small image i: uniqueness and fine existence. Proceedings of the London Mathematical Society, 3(2):371–406, 1990.
Krauthgamer & Lee [2006] Krauthgamer, R. and Lee, J. R. Algorithms on negatively curved spaces. In 2006 47th Annual IEEE Symposium on Foundations of Computer Science (FOCS’06), pp. 119–132. IEEE, 2006.
Lee [2013] Lee, J. M. Smooth manifolds. In Introduction to Smooth Manifolds, pp. 1–31. Springer, 2013.
Liu et al. [2019] Liu, Q., Nickel, M., and Kiela, D. Hyperbolic graph neural networks. In Advances in Neural Information Processing Systems, pp. 8230–8241, 2019.
Monath et al. [2019] Monath, N., Zaheer, M., Silva, D., McCallum, A., and Ahmed, A. Gradient-based hierarchical clustering using continuous representations of trees in hyperbolic space. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 714–722, 2019.
Nickel & Kiela [2017] Nickel, M. and Kiela, D. Poincaré embeddings for learning hierarchical representations. In Advances in neural information processing systems, pp. 6338–6347, 2017.
Pennec [2018] Pennec, X. Barycentric subspace analysis on manifolds. Annals of Statistics, 46(6A):2711–2746, 2018.
Sala et al. [2018] Sala, F., De Sa, C., Gu, A., and Ré, C. Representation tradeoffs for hyperbolic embeddings. In International Conference on Machine Learning, pp. 4460–4469, 2018.
Sarkar [2011] Sarkar, R. Low distortion delaunay embedding of trees in hyperbolic plane. In International Symposium on Graph Drawing, pp. 355–366. Springer, 2011.
Sommer [2013] Sommer, S. Horizontal dimensionality reduction and iterated frame bundle development. In International Conference on Geometric Science of Information, pp. 76–83. Springer, 2013.
Tay et al. [2018] Tay, Y., Tuan, L. A., and Hui, S. C. Hyperbolic representation learning for fast and efficient neural question answering. In Proceedings of the Eleventh ACM International Conference on Web Search and Data Mining, pp. 583–591, 2018.
Thurston [1978] Thurston, W. P. Geometry and topology of three-manifolds, 1978. Lecture notes. Available at http://library.msri.org/books/gt3m/.
Tifrea et al. [2018] Tifrea, A., Becigneul, G., and Ganea, O.-E. Poincare glove: Hyperbolic word embeddings. In International Conference on Learning Representations, 2018.
Tran & Vu [2008] Tran, D. A. and Vu, K. Dimensionality reduction in hyperbolic data spaces: Bounding reconstructed-information loss. In Seventh IEEE/ACIS International Conference on Computer and Information Science (icis 2008), pp. 133–139. IEEE, 2008.
Wang [2021] Wang, M.-X. Laplacian eigenspaces, horocycles and neuron models on hyperbolic spaces, 2021. URL https://openreview.net/forum?id=ZglaBL5inu.
Wu & Charikar [2020] Wu, X. and Charikar, M. Nearest neighbor search for hyperbolic embeddings. arXiv preprint arXiv:2009.00836, 2020.
Zhang & Fletcher [2013] Zhang, M. and Fletcher, T. Probabilistic principal geodesic analysis. In Advances in Neural Information Processing Systems, pp. 1178–1186, 2013.

Appendix A Horospherical Projection: Proofs and Discussions

We show that the horospherical projection $\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}$ in Section 3.2.2 is well-defined, shares many nice properties with Euclidean orthogonal projections, and has a closed-form expression.

More specifically, this section is organized as follows. In Section A.1, we show that horospherical projections are well-defined (Theorem A.4). In Section A.2, we show that they are base-point independent (A.7), non-expanding (A.12), and distance-preserving along a family of $K$ -dimensional submanifolds (A.10). Finally, in Section A.4, we explain how they can be computed efficiently in the hyperboloid model.

For an illustration in the case of $K=2$ ideal points in $\mathbb{H}^{3}$ , see Fig. 5. This figure might help explain the intuitions behind the theorems in Section A.1 and Section A.2.

A.1 Well-definedness

Recall that given a base point $b\in\mathbb{H}^{d}$ and $K>1$ ideal points $\{p_{1},\dots,p_{K}\}$ , we would like to define $\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}$ by

x\to M\cap S(p_{1},x)\cap S(p_{2},x)\cap\dots\cap S(p_{K},x),

where $M=\operatorname{GH}(b,p_{1},\dots,p_{K})$ is the target submanifold and $S(p_{j},x)$ is the horosphere centered at $p_{j}$ and passing through $x$ . For this definition to make sense, the intersection in the right hand side must contain exactly one point for each $x\in\mathbb{H}^{d}$ . Unfortunately, this is not the case: In fact, the intersection generally consists of two points. Nevertheless, we will show that there is a consistent way to choose one from these two points, making the function $\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}$ well-defined. This is the result of Theorem A.4.

First, to understand the above intersection, we give a more concrete description of $\cap_{j}S(p_{j},x)$ .

Lemma A.1.

Let $P=\operatorname{GH}(p_{1},p_{2},\dots,p_{K})$ . Then for every $x\in\mathbb{H}^{d}$ , the intersection of horospheres

S(x)=S(p_{1},x)\cap S(p_{2},x)\cap\dots\cap S(p_{k},x)

is precisely the orbit of $x$ under the group $G$ of rotations around $P$ .

Proof.

First, note that every rotation around $P$ preserves the horospheres $S(p_{j},x)$ - just like how every rotation around an axis preserves every sphere whose center is on that axis. It follows that $S(x)$ is preserved by $G$ . In particular, the orbit of $x$ under $G$ is contained in $S(x)$ .

It remains to show that $S(x)$ contains no other points. To this end, consider any $y\neq x$ in $S(x)$ . The perpendicular bisector $B$ of $x$ and $y$ is a totally geodesic hyperplane of $\mathbb{H}^{d}$ that contains every $p_{j}$ (because each $p_{j}$ is intuitively the center of a sphere that goes through $x$ and $y$ ). Thus, by the definition of geodesic hull, $B\supset P$ . In particular, the reflection through $B$ sends $x$ to $y$ and fixes every point in $P$ .

Now take any geodesic hyperplane $A$ that contains both $P$ and $y$ , so that the reflection through $A$ fixes $y$ and every point in $P$ . Then the composition of the reflections through $B$ and $A$ is a rotation that sends $x$ to $y$ and fixes every point in $P$ . In other words, it is a rotation around $P$ that sends $x$ to $y$ . Therefore, $y$ belongs to the orbit of $x$ under $G$ . ∎

Corollary A.2.

If $x\in P$ then $S(x)=\{x\}$ . Otherwise, let $\pi^{\operatorname{G}}_{P}(x)$ be the geodesic projection of $x$ onto $P$ , and $Q(x)$ be the geodesic submanifold that orthogonally complements $P$ at $\pi^{\operatorname{G}}_{P}(x)$ . Then $S(x)\subset Q(x)$ and is precisely the (hyper)sphere in $Q(x)$ that is centered at $\pi^{\operatorname{G}}_{P}(x)$ and passing through $x$ .

Proof.

This follows from Lemma A.1. If $x\in P$ then every rotation around $P$ fixes $x$ , so the orbit of $x$ is just itself.

Now consider the case $x\not\in P$ . All rotations around $P$ must preserve $\pi^{\operatorname{G}}_{P}(x)$ and the orthogonal complement $Q(x)$ of $P$ at $\pi^{\operatorname{G}}_{P}(x)$ . Furthermore, when restricted to the space $Q(x)$ , these rotations are precisely the rotations in $Q(x)$ around the point $\pi^{\operatorname{G}}_{P}(x)$ . Thus, for every $y\in Q(x)$ , the orbit of $y$ under $G$ is a sphere in $Q(x)$ centered at $\pi^{\operatorname{G}}_{P}(x)$ . In particular, $S(x)$ , which is the orbit of $x$ , is the sphere in $Q(x)$ that is centered at $\pi^{\operatorname{G}}_{P}(x)$ and passing through $x$ . ∎

A.2 gives the following characterization of the intersection

M\cap S(p_{1},x)\cap S(p_{2},x)\cap\dots\cap S(p_{K},x)=M\cap S(x):

Note that $P=\operatorname{GH}(p_{1},\dots,p_{K})$ is a geodesic submanifold of $M=\operatorname{GH}(b,p_{1},\dots,p_{K})$ and that $\dim P=\dim M-1$ . Thus, through every point $y\in P$ , there is a unique geodesic $\alpha$ on $M$ that goes through $y$ and is perpendicular to $P$ .

Corollary A.3.

If $x\in P$ then $M\cap S(x)=\{x\}$ . Otherwise, let $\alpha$ be the geodesic on $M$ that goes through $\pi^{\operatorname{G}}_{P}(x)$ and is perpendicular to $P$ . Then $M\cap S(x)$ consists of the two points on $\alpha$ whose distance to $\pi^{\operatorname{G}}_{P}(x)$ equals $d_{\mathbb{H}}(x,\pi^{\operatorname{G}}_{P}(x))$ .

Proof.

The case $x\in P$ is clear since $x\in M$ and $S(x)=\{x\}$ by A.2.

For the other case, let $Q(x)$ be the orthogonal complement of $P$ at $\pi^{\operatorname{G}}_{P}(x)$ . Then by A.2, $S(x)$ is precisely the sphere in $Q(x)$ that is centered at $c(x)$ and passing through $x$ .

Now note that $M\cap Q(x)=\alpha$ . Since $S(x)\subset Q(x)$ , this gives $M\cap S(x)=M\cap Q(x)\cap S(x)=\alpha\cap S(x)$ . We know that every sphere intersects every line through the center at two points. ∎

Therefore, to define $\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}(x)$ , we just need to choose one of the two poins in $M\cap S(x)$ in a consistent way (so that the map is differentiable). To this end, note that $P$ cuts $M$ into two half-spaces, and exactly one of them contains the base point $b$ . (Recall that $b\in M$ and $b\not\in P$ by the “independence” condition). We denote this half by $P_{b}$ . Then, while $S(x)$ contains two points in $M$ , it only contains one point in $P_{b}$ :

Theorem A.4.

Let $\alpha$ be the geodesic on $M$ that goes through $\pi^{\operatorname{G}}_{P}(x)$ and is perpendicular to $P$ . Let $\alpha^{+}$ be the half of $\alpha$ contained in $P_{b}$ . Then $\alpha^{+}$ intersects the sphere $S(x)$ at a unique point $x^{\prime}$ , which is also the unique intersection point between $P_{b}$ and $S(x)$ . Thus, we can define $\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}$ by

\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}(x)=\alpha^{+}\cap S(x)=P_{b}\cap S(x).

Equivalently, $\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}(x)$ is the point in $M\cap S(x)$ that is strictly closer to $b$ .

Proof.

By A.2, $S(x)$ is a sphere centered at $\pi^{\operatorname{G}}_{P}(x)$ . Then, since $\alpha^{+}$ is a geodesic ray starting at $\pi^{\operatorname{G}}_{P}(x)$ , it must intersect $S(x)$ at a unique point $x^{\prime}$ .

Next, we have

P_{b}\cap S(x)=P_{b}\cap M\cap S(x)=P_{b}\cap\alpha\cap S(x)=\alpha^{+}\cap S(x),

where the first equality holds because $P_{b}\subset M$ , the second because $M\cap S(x)=\alpha\cap S(x)$ by the proof of A.3, and the third because $P_{b}\cap\alpha=\alpha^{+}$ . Thus, $P_{b}\cap S(x)$ is precisely $x^{\prime}$ .

Finally, let $x^{\prime\prime}$ be the other point of $M\cap S(x)=\alpha\cap S(x)$ . Then $P$ is the perpendicular bisector of $x^{\prime}$ and $x^{\prime\prime}$ in $M$ . Thus, every point on the same side of $P$ in $M$ as $x^{\prime}$ (but not on the boundary $P$ ) is strictly closer to $x^{\prime}$ than to $x^{\prime\prime}$ . By definition, $b$ is one of such point. Thus, $\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}(x)$ is the point in $M\cap S(x)$ that is closer to $b$ . ∎

A.2 Geometric properties

From Lemma A.1 and Theorem A.4, we obtain another interpretation of $\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}$ : It maps $\mathbb{H}^{d}$ to $P_{b}\subset M$ by rotating every point $x\in\mathbb{H}^{d}$ around $P$ until it hits $P_{b}$ . In other words, we have

Theorem A.5 (The “open book” interpretation).

For any $x\not\in P$ , let $M_{x}=\operatorname{GH}(P\cup\{x\})$ . Then $P$ cuts $M_{x}$ into two half-spaces; we denote the half that contains $x$ by $P_{x}$ . Then, when restricted to $P_{x}$ , the map

\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}:P_{x}\to P_{b}

is simply a rotation around $P$ . ∎

Following this, we call $P$ the spine of the horosphere projection. The identity $\mathbb{H}^{d}=\cup_{x\not\in P}P_{x}$ can be thought of as an open book decomposition of $\mathbb{H}^{d}$ into pages $P_{x}$ that are bounded by the spine $P$ . The horosphere projection $\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}$ then simply acts by collapsing every page onto a specified page $P_{b}$ .

Here are some consequences of this interpretation:

Corollary A.6.

$\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}$ only depends on the spine $P$ and not specifically on $p_{1},\dots p_{K}$ . Thus, when we are not interested in the specific ideal points, we simply write $\pi^{\operatorname{H}}_{b,P}$ .

Proof.

As noted above, $\pi^{\operatorname{H}}_{b,p_{1},\dots,p_{K}}(x)$ can be obtained by rotating $x$ around $P$ until it hits $P_{b}$ . This operation does not used the exact positions of $p_{j}$ at all. ∎

Corollary A.7.

The choice of $b$ does not affect the geometry of the projection $\pi^{\operatorname{H}}_{b,P}$ . More precisely, for any two base points $b,b^{\prime}\not\in P$ , the horosphere projections

\pi^{\operatorname{H}}_{b,P}:\mathbb{H}^{d}\to P_{b}\ \ \text{and}\ \ \pi^{\operatorname{H}}_{b^{\prime},P}:\mathbb{H}^{d}\to P_{b^{\prime}}

only differ by a rotation $P_{b}\to P_{b^{\prime}}$ around $P$ .

Proof.

By Theorem A.5, when restricted to any page $P_{x}$ , the maps $\pi^{\operatorname{H}}_{b,P}:\mathbb{H}^{d}\to P_{b}$ and $\pi^{\operatorname{H}}_{b^{\prime},P}:\mathbb{H}^{d}\to P_{b^{\prime}}$ are just rotations around $P$ . The difference between any two rotations around $P$ is another rotation around $P$ . ∎

In particular, A.7 implies Theorem 4.1: See 4.1

As discussed in Section 4, this theorem helps reduces parameters and simplifies the computation of $\pi^{\operatorname{H}}_{b,P}$ .

Remark A.8.

A.6 and A.7 together imply that the HoroPCA algorithm (5) only depends on the geodesic hulls $\operatorname{GH}(p_{1}),\operatorname{GH}(p_{1},p_{2}),\dots\operatorname{GH}(p_{1},\dots,p_{K})$ of the ideal points and not the specific ideal points themselves. It follows that, theoretically, the search space of (5) has dimension $dK-\frac{1}{2}K(K+1)$ – the same as the dimension of the space of flags in Euclidean spaces.

In our implementation, for simplicity we parametrize the $K$ ideal points independently, which results in a suboptimal search space dimension $(d-1)K$ . Nevertheless, this is still slightly more efficient than the parametrizations used in PGA and BSA, which require have $(d+1)K$ -dimensions.

The following corollaries say that horospherical projections share a nice property with Euclidean orthogonal projections: When projecting to a $K$ -dimensional submanifold, they preserve the distances along $K$ dimensions and collapse the distances along the other $d-K$ orthogonal dimensions:

Corollary A.9.

$\pi^{\operatorname{H}}_{b,P}$ is invariant under rotations around $P$ . In other words, if a rotation around $P$ takes $x$ to $y$ then $\pi^{\operatorname{H}}_{b,P}(x)=\pi^{\operatorname{H}}_{b,P}(y)$ .

Consequently, every $x\not\in P$ belongs to a $(d-K)$ -dimensional submanifold that is collapsed to a point by $\pi^{\operatorname{H}}_{b,P}$ .

Proof.

The open book interpretation tells us that $\pi^{\operatorname{H}}_{b,P}(x)$ and $\pi^{\operatorname{H}}_{b,P}(y)$ are precisely the intersections of $P_{b}$ with $S(x)$ and $S(y)$ , respectively. If $y$ belongs to the rotation orbit $S(x)$ of $x$ then the rotation orbit $S(y)$ of $y$ is the same as $S(x)$ . Thus $\pi^{\operatorname{H}}_{b,P}(x)=\pi^{\operatorname{H}}_{b,P}(y)$ .

Hence, for every $x\in\mathbb{H}^{d}$ , $S(x)$ is collapsed to a point by $\pi^{\operatorname{H}}_{b,P}$ . To see that $\dim S(x)=d-K$ when $x\not\in P$ , recall that by A.2, if $Q(x)$ is the orthogonal complement of $P$ at $\pi^{\operatorname{G}}_{P}(x)$ then $S(x)$ is a hypersphere inside $Q(x)$ . Since the ideal points $p_{j}$ are assumed to be “affinely independent,” we have $\dim P=K-1$ , so $\dim Q(x)=d-(K-1)$ , and $\dim S(x)=\dim Q(x)-1=d-K$ . ∎

Corollary A.10.

For every $x\in\mathbb{H}^{d}$ , there exists a $K$ -dimensional totally geodesic submanifold (with boundary) that contains $x$ and is mapped isometrically to $P_{b}$ by $\pi^{\operatorname{H}}_{b,P}$ . If $x\not\in P$ then such a manifold is unique.

Proof.

If $x\not\in P$ then the submanifold $P_{x}$ in Theorem A.5 is a geodesic submanifold that contains $x$ and is mapped isometrically to $P_{b}$ by $\pi^{\operatorname{H}}_{b,P}$ . As in the proof of A.9, we have $\dim P=K-1$ and $\dim P_{x}=\dim P+1=K$ .

Since A.9 implies that the other $d-K$ dimensions are collapsed by $\pi^{\operatorname{H}}_{b,P}$ , no other distances from $x$ can be preserved. Thus, $P_{x}$ is the unique submanifold with the desired properties.

If $x\in P$ then $x\in P_{y}$ for every $y\not\in P$ . Note that this means every distance from $x$ is preserved by $\pi^{\operatorname{H}}_{b,P}$ . ∎

The following corollaries say that, like Euclidean orthogonal projections, horospherical projections never increase distances. Thus, minimizing distortion is roughly equivalent to maximizing projected distances. This is another motivation for Eq. 4.

Corollary A.11.

For every $x\in\mathbb{H}^{d}$ and every tangent vector $\vec{v}$ at $x$ ,

\|\pi^{\operatorname{H}}_{b,P}(\vec{v})\|_{\mathbb{H}}\leq\|\vec{v}\|_{\mathbb{H}}.

Proof.

This follows from A.10 and A.9:

If $x\in P$ then the proof of A.10 implies that $\pi^{\operatorname{H}}_{b,P}$ preserves every distance from $x$ . Thus, the desired inequality is actually an equality.

If $x\not\in P$ then by $\vec{v}$ has an orthogonal decomposition $\vec{v}=\vec{u}+\vec{u}^{\perp}$ , where $\vec{u}$ and $\vec{u}^{\perp}$ are tangent and perpendicular to $P_{x}$ , respectively. By A.10 and A.9, $\pi^{\operatorname{H}}_{b,P}$ preserves the length of $\vec{u}$ while collapsing $\vec{u}^{\perp}$ to $0$ . It follows that $\|\pi^{\operatorname{H}}_{b,P}(\vec{v})\|_{\mathbb{H}}=\|\vec{u}\|_{\mathbb{H}}\leq\|\vec{v}\|_{\mathbb{H}}$ . ∎

Corollary A.12.

$\pi^{\operatorname{H}}_{b,P}$ is non-expanding. In other words, for every $x,y\in\mathbb{H}^{d}$ ,

d_{\mathbb{H}}(\pi^{\operatorname{H}}_{b,P}(x),\pi^{\operatorname{H}}_{b,P}(y))\leq d_{\mathbb{H}}(x,y).

Proof.

We first show that for any path $\gamma(t)$ in $\mathbb{H}^{d}$ ,

\operatorname{length}(\pi^{\operatorname{H}}_{b,P}(\gamma))\leq\operatorname{length}(\gamma).

Indeed, note that the velocity vector of $\pi^{\operatorname{H}}_{b,P}(\gamma)$ is precisely $\pi^{\operatorname{H}}_{b,P}(\dot{\gamma})$ . Then, by A.11,

	$\displaystyle\operatorname{length}(\pi^{\operatorname{H}}_{b,P}(\gamma))$	$\displaystyle=\int\\|\pi^{\operatorname{H}}_{b,P}(\dot{\gamma}(t))\\|_{\mathbb{H}}dt$
		$\displaystyle\leq\int\\|\dot{\gamma}(t)\\|_{\mathbb{H}}dt=\operatorname{length}(\gamma).$

Now for any $x,y\in\mathbb{H}^{d}$ , let $\gamma(t)$ be the geodesic segment from $x$ to $y$ . Then $\pi^{\operatorname{H}}_{b,P}(\gamma)$ is a path connecting $\pi^{\operatorname{H}}_{b,P}(x)$ and $\pi^{\operatorname{H}}_{b,P}(y)$ . Thus, the projected distance is at most the length of $\pi^{\operatorname{H}}_{b,P}(\gamma)$ , which by the above argument is at most $\operatorname{length(\gamma)}=d_{\mathbb{H}}(x,y).$ ∎

A.3 Detour: A Review of the Hyperboloid Model

Theorem A.4 suggests that $\pi^{\operatorname{H}}_{b,P}(x)$ can be computed in three steps:

1.

Find the geodesic projection $\pi^{\operatorname{G}}_{P}(x)$ of $x$ onto $P$ .
2.

Find a geodesic ray $\alpha^{+}$ on $P_{b}$ that starts at $\pi^{\operatorname{G}}_{P}(x)$ and is orthogonal to $P$ .
3.

Return the unique point on $\alpha^{+}$ that is of distance $d_{\mathbb{H}}(x,\pi^{\operatorname{G}}_{P}(x))$ from $\pi^{\operatorname{G}}_{P}(x)$ .

It turns out that these subroutines are easier to implement in the hyperboloid model instead of the Poincaré model of hyperbolic spaces. Thus, we first briefly review the basic definitions and properties of this model. The readers who are already familiar with the hyperboloid model can skip to Section A.4, where we describe the full algorithm.

For a more detailed treatment, see [31].

Remark A.13.

The above steps are slightly different from the ones mentioned in Section 3.2.2. However, by Theorem A.4, these two descriptions are equivalent. We decided to use the “closer-to- $b$ ” description in Section 3.2.2 because it is slightly more self-contained, but the actual implementation will be based on the above three steps.

Minkowski spaces

We first describe Minkowski spaces, which is the embedding space where the hyperboloid model sits in as hypersurfaces.

The $(d+1)$ -dimensional Minkowski space $\mathbb{R}^{1,d}$ is like the flat Euclidean space $\mathbb{R}^{1+d}$ except with a dot product that has a negative sign in the first coordinate. More precisely, $\mathbb{R}^{1,d}$ is the vector space $\mathbb{R}^{1+d}$ equipped with the indefinite, non-degenerate bilinear form

B\left((t,x_{1},\dots,x_{d}),(u,y_{1},\dots,y_{d})\right)=-tu+\sum_{i=1}^{d}x_{i}y_{i},

which serves as the “dot product.”

Like with Euclidean spaces, the quantity $B(\vec{v},\vec{v})$ is called the (Minkowski) squared norm of the vector $\vec{v}$ . Two vectors $\vec{u},\vec{v}\in\mathbb{R}^{1,d}$ are called orthogonal if $B(\vec{u},\vec{v})=0$ . The orthogonal complement of a linear subspace $V\subset\mathbb{R}^{1,d}$ is the set $V^{\perp}=\{\vec{w}\in\mathbb{R}^{1,d}:B(\vec{w},\vec{v})=0\text{ for every }\vec{v}\in V\}$ , which still has dimension $\dim\mathbb{R}^{1,d}-\dim V$ . However, unlike in Euclidean spaces, vectors in $\mathbb{R}^{1,d}$ can have negative or zero squared norms, and linear subspaces can intersect their orthogonal complements.

Types of vectors in $\mathbb{R}^{1,d}$

Vectors with negative, zero, and positive (Minkowski) squared norms are called time-like, light-like, and space-like, respectively. Time-like and light-like vectors together form a solid double cone in $\mathbb{R}^{1,d}$ .

If $\vec{v}$ is a time-like or light-like vector, we call it future-pointing if its first coordinate is positive, otherwise we call it past-pointing. It follows from Cauchy-Schwarz inequality that.

Proposition A.14.

If $\vec{u},\vec{v}$ are future-pointing time-like or light-like vectors then $B(\vec{u},\vec{v})<0$ .

A linear subspace $V$ of $\mathbb{R}^{1,d}$ is called space-like if every non-zero vector in $V$ is space-like. In that case, the bilinear form $B(\cdot,\cdot)$ restricts to a positive-definite bilinear form on $V$ , thus making it isometric to an Euclidean vector space. It follows from A.14 that

Proposition A.15.

The orthogonal complement of a time-like vector is a space-like linear subspace.

The hyperboloid model of hyperbolic spaces

We are now ready to introduce the hyperboloid model $\mathbb{H}^{d}$ . It sits inside $\mathbb{R}^{1,d}$ in a similar way to how the unit sphere $\mathbb{S}^{d}$ sits inside the Euclidean space $\mathbb{R}^{1+d}$ .

Definition A.16.

The hyperboloid model of $d$ -dimensional hyperbolic spaces is the set $\mathbb{H}^{d}$ of future-pointing vectors in $\mathbb{R}^{1,d}$ with Minkowsi squared norm $-1$ .

Remark A.17.

For the rest of Appendix A, we will use $\mathbb{H}^{d}$ to denote this hyperboloid model as a subset of $\mathbb{R}^{1,d}$ , and not the abstract hyperbolic space or its other models (e.g. Poincaré ball). This should not lead to any ambiguities because we will not work with any other model.

The following properties of $\mathbb{H}^{d}$ further illustrate the analogy with spheres in Euclidean spaces.

Proposition A.18.

For every $\vec{x}\in\mathbb{H}^{d}$ , the tangent space $T_{\vec{x}}\mathbb{H}^{d}$ of $\mathbb{H}^{d}$ at $\vec{x}$ is (parallel to) the orthogonal complement of $\vec{x}$ in $\mathbb{R}^{1,d}$ .

Proposition A.19.

When restricted to each tangent space of $\mathbb{H}^{d}$ , the bilinear form $B(\cdot,\cdot)$ is positive definite. This defines a Riemannian metric on $\mathbb{H}^{d}$ which has constant curvature $-1$ . The stereographic projection

(t,x_{1},\dots,x_{d})\to\left(\frac{x_{1}}{1+t},\dots,\frac{x_{d}}{1+t}\right)

is an isometry between $\mathbb{H}^{d}$ and the Poincaré model. Its inverse map is given by

(y_{1},\dots,y_{d})\to\frac{(1+\sum_{i}y_{i}^{2},2y_{1},\dots,2y_{d})}{1-\sum_{i}y_{i}^{2}}.

Proposition A.20.

Every $k$ -dimensional geodesic submanifold of $\mathbb{H}^{d}$ is the intersection of $\mathbb{H}^{d}$ with a $(k+1)$ -dimensional linear subspace of $\mathbb{R}^{1,d}$ .

In the hyperboloid model, the ideal points of hyperbolic spaces are represented by light-like directions (instead of individual vectors):

Proposition A.21.

Each ideal point of $\mathbb{H}^{d}$ is represented by a $1$ -dimensional linear subspace spanned by some $\vec{v}\in\mathbb{R}^{1,d}$ with $B(\vec{v},\vec{v})=0$ . The map

(t,x_{1},\dots,x_{d})\to\left(\frac{x_{1}}{t},\dots,\frac{x_{d}}{t}\right)

gives a correspondence between a light-like vector $\vec{v}\in\mathbb{R}^{1,d}$ and an ideal point $p\in\mathbb{S}_{\infty}^{d-1}$ in the Poincaré model that is represented by $\operatorname{span}(\vec{v})$ . This correspondence is compatible with the stereographic projection in A.19. Its inverse map is given by

(y_{1},\dots,y_{d})\to(1,y_{1},\dots,y_{d}).

We conclude this section by noting that geodesic hulls in the hyperboloid model are closely related to linear spans in Minkowski space:

Proposition A.22.

Let $S$ be a set of vectors that are either in $\mathbb{H}^{d}$ or represent ideal directions of $\mathbb{H}^{d}$ . Then the geodesic hull of $S$ in $\mathbb{H}^{d}$ is the intersection of $\mathbb{H}^{d}$ with the linear span of $S$ .

In particular, since the spine $P$ in our setting (Theorem A.4) is the geodesic hull of some ideal points, it is cut out by the linear span of the corresponding ideal directions.

A.4 Computation of Horospherical Projections

Algorithm 1 Horospherical Projection

Input: point

x

, ideal points

\{p_{1},\dots,p_{K}\}

, base point

b

z\leftarrow\pi^{\operatorname{Mink}}_{\operatorname{span}(p_{1},\dots,p_{K})}(x)

{orthogonal projection in ambient space}

w\leftarrow z/\|z\|

{rescale to a vector on the hyperboloid}

u\leftarrow b-w

{subtraction in ambient space}

u\leftarrow u-\pi^{\operatorname{Mink}}_{\operatorname{span}(p_{1},\dots,p_{K})}(u)

{orthogonal projection in ambient space}

u\leftarrow u/\|u\|

{make

u

the unit tangent vector orthogonal to the spine

P

}

y\leftarrow\exp_{w}(d_{\mathbb{H}}(x,w)\cdot u)

{exponential map in

\mathbb{H}^{d}

}

return

y

Now we describe how the three steps mentioned at the beginning of Section A.3 can be implemented in the hyperboloid model. The results of this section are summarized in Algorithm 1. To transfer back and forth between the hyperboloid and Poincaré models, we use the formulas in A.19 and A.21.

Step 1: Computing geodesic projections

We first describe how to compute the geodesic (or closest-point) projection from $\mathbb{H}^{d}$ to a geodesic submanifold $P$ in $\mathbb{H}^{d}$ . This process is very similar to how geodesic projections work in Euclidean spheres: We first perform an orthogonal projection onto the linear subspace $V$ that cuts out $P$ , then rescale the result to get a vector on $\mathbb{H}^{d}$ .

Generally, orthogonal projections in Minkowski spaces are very similar to those in Euclidean spaces. However, since vectors in $\mathbb{R}^{1,d}$ can have norm zero, the orthogonal projection $\pi^{\operatorname{Mink}}_{V}$ is not well-defined for every subspace $V$ . Thus, a little extra argument is needed.

Proposition A.23 (Orthogonal projections onto time-containing linear subspaces).

Let $V$ be a linear subspace of $\mathbb{R}^{1,d}$ that contains some time-like vectors. Then

1.

$V\cap V^{\perp}=\{0\}$ . Consequently, we can define a linear orthogonal projection $\pi^{\operatorname{Mink}}_{V}:\mathbb{R}^{1,d}\to V$ as follows:

Since $\mathbb{R}^{1,d}=V\oplus V^{\perp}$ , every vector $\vec{x}\in\mathbb{R}^{1,d}$ can be uniquely written as $\vec{x}=\vec{z}+\vec{n}$ for some $\vec{z}\in V$ and $\vec{n}\in V^{\perp}$ . Then, we let $\pi^{\operatorname{Mink}}_{V}(\vec{x})\coloneqq\vec{z}.$
2.

Let $A$ be a matrix whose column vectors form a linear basis of $V$ . Let $B$ be the $(1+d)\times(1+d)$ symmetric matrix associated to the bilinear form $B$ , i.e. the diagonal matrix with entries $(-1,1,1,1,\dots,1)$ . Then $A^{\top}BA$ is non-singular, and the linear projection $\pi^{\operatorname{Mink}}_{V}$ is given by $\vec{x}\to A(A^{\top}BA)^{-1}A^{\top}B\vec{x}$
3.

If $\vec{x}$ is a future-pointing time-like vector then so is $\pi^{\operatorname{Mink}}_{V}(\vec{x})$ .

Proof.

1.

Let $\vec{v}$ be a time-like vector in $V$ . Then $\vec{v}^{\perp}$ is space-like by A.15. Since $V\cap V^{\perp}\subset V^{\perp}\subset\vec{v}^{\perp}$ , the intersection $V\cap V^{\perp}$ must be space-like. On the other hand, for every $\vec{w}\in V\cap V^{\perp}$ , we have $B(\vec{w},\vec{w})=0$ , so $\vec{w}$ must be light-like. It follows that $V\cap V^{\perp}=\{0\}$ . The $\mathbb{R}^{1,d}=V\oplus V^{\perp}$ part of the claim is a standard linear algebra fact.
2.

This follows from the same argument that deduces the formula for Euclidean orthogonal projections.
3.

To avoid cumbersome notations, let $\vec{z}=\pi^{\operatorname{Mink}}_{V}(\vec{x})$ . Then since $\vec{x}-\vec{z}$ is orthogonal to $V$ and in particular to $\vec{z}$ , we have the Pythagorean formula

$B(\vec{z},\vec{z})+B(\vec{x}-\vec{z},\vec{x}-\vec{z})=B(\vec{x},\vec{x}).$

On the other hand, since $\vec{x}-\vec{z}$ is orthogonal to $V$ and in particular orthogonal to some time-like vectors in $V$ , by A.15, it must be space-like. Thus, if $B(\vec{x},\vec{x})<0$ then the above equation implies $B(\vec{z},\vec{z})<0$ , which means $\vec{z}$ is time-lile.

Now either $\vec{z}$ or $-\vec{z}$ is future-pointing. In the latter case, since $\vec{x}$ is future-pointing, A.14 implies $B(-\vec{z},\vec{x})<0$ . On the other hand, we have $B(\vec{z},\vec{x})=B(\vec{z},\vec{z})<0$ , contradicting the above inequality. Thus, $\vec{z}$ is future-pointing. ∎

Proposition A.24 (Geodesic projections in $\mathbb{H}^{d}$ ).

Let $P$ be a geodesic submanifold of $\mathbb{H}^{d}$ . Recall that by A.20, $P=\mathbb{H}^{d}\cap V$ for some linear subspace $V$ of $\mathbb{R}^{1,d}$ . Then for every $\vec{x}\in\mathbb{H}^{d}$ , the geodesic projection $\pi^{\operatorname{G}}_{P}(\vec{x})$ of $\vec{x}$ onto $P$ in $\mathbb{H}^{d}$ is given by

\pi^{\operatorname{G}}_{P}(\vec{x})=\frac{\vec{z}}{\sqrt{-B(\vec{z},\vec{z})}},

where $\vec{z}=\pi^{\operatorname{Mink}}_{V}(\vec{x})$ is the linear orthogonal projection of $\vec{x}$ onto $V$ .

Proof.

Again, to avoid cumbersome notations, let

\vec{w}=\frac{\vec{z}}{\sqrt{-B(\vec{z},\vec{z})}}.

We will show that $\vec{w}=\pi^{\operatorname{G}}_{P}(\vec{x})$ . First, note that since $\vec{z}$ is a future-pointing time-like vector by A.23, $\vec{w}\in\mathbb{H}^{d}$ . Since $\vec{w}\in V$ , it follows that $\vec{w}\in P$ .

Let $W$ be the linear span of $\vec{w}$ and $\vec{x}$ , so that $W\cap\mathbb{H}^{d}$ is the geodesic $\gamma$ in $\mathbb{H}^{d}$ that connects $\vec{w}$ and $\vec{x}$ . Note that $\vec{w}$ and $\vec{x}-\vec{z}$ form a basis of $W$ . They are both orthogonal to the tangent space $T_{\vec{w}}P$ of $P$ at $\vec{w}$ because:

•

By A.18, $\vec{w}$ is orthogonal to every tangent vector of $\mathbb{H}^{d}$ at $\vec{w}$ .
•

$T_{\vec{w}}P$ is contained in $V$ , which is orthogonal to $\vec{x}-\vec{z}$ .

Thus, $W$ is orthogonal to $T_{\vec{w}}P$ . It follows that the geodesic $\gamma$ is orthogonal to $P$ , which means $\vec{w}=\pi^{\operatorname{G}}_{P}(\vec{x})$ . ∎

Step 2: Finding orthogonal geodesic ray

Recall that if $P$ is a geodesic submanifold of $\mathbb{H}^{d}$ and $\vec{b}\in\mathbb{H}^{d}$ is a point not in $P$ , then the geodesic hull $M$ of $P\cup\{\vec{b}\}$ in $\mathbb{H}^{d}$ is a geodesic submanifold of $\mathbb{H}^{d}$ with dimension $\dim P+1$ . Thus, $P$ cuts $M$ into two halves, and we denote the half that contains $\vec{b}$ by $P_{b}$ .

Given a point $\vec{w}\in P$ , there exists a unique geodesic ray $\alpha^{+}$ that starts at $\vec{w}$ , stays on $P_{b}$ , and is orthogonal to $P$ . In this section, we describe how to compute $\alpha^{+}$ .

Recall that by A.20, $P=V\cap\mathbb{H}^{d}$ for some linear subspace $V$ .

Proposition A.25.

The vector

\vec{u}=(\vec{b}-\vec{w})-\pi^{\operatorname{Mink}}_{V}(\vec{b}-\vec{w})

is tangent to $M$ and orthogonal to $P$ at $\vec{w}$ . Furthermore, it points toward the side of $P_{b}$ .

Proof.

By construction, $\vec{u}$ is orthogonal to $V$ . Note that $V$ contains both $\vec{w}$ and the tangent space $T_{\vec{w}}P$ of $P$ at $\vec{w}$ . Thus, $\vec{u}$ is orthogonal to both of them. Together with A.18, this implies that $\vec{u}$ is a tangent vector of $\mathbb{H}^{d}$ at $\vec{w}$ that is orthogonal to $P$ .

Next, note that $M$ is the intersection of $\mathbb{H}^{d}$ with the linear span of $V\cup\{\vec{b}\}$ , and since $\vec{b}$ and $\vec{w}$ both belong to this linear span, so does $\vec{u}$ . Thus, since $\vec{u}$ is a tangent vector of $\mathbb{H}^{d}$ at $\vec{w}$ , it must be tangent to $M$ .

By construction, $\vec{u}$ points from $\vec{w}$ toward the side of $\vec{b}$ instead of away from it. More rigorously, note that by essentially the same argument as above,

\vec{a}\coloneqq(\vec{b}-\vec{w})-\pi^{\operatorname{Mink}}_{\vec{w}}(\vec{b}-\vec{w})

is a tangent vector at $\vec{w}$ of the geodesic on $\mathbb{H}^{d}$ that goes from $\vec{w}$ to $\vec{b}$ . Also, note that

\vec{a}-\vec{u}=\pi^{\operatorname{Mink}}_{V}(\vec{b}-\vec{w})-\pi^{\operatorname{Mink}}_{\vec{w}}(\vec{b}-\vec{w}).

Since $V=T_{\vec{w}}P\oplus\operatorname{span}(\vec{w})$ is an orthogonal decomposition, this means

\vec{a}-\vec{u}=\pi^{\operatorname{Mink}}_{T_{\vec{w}}P}(\vec{b}-\vec{w}),

which in particular implies $\vec{a}-\vec{u}\in T_{\vec{w}}P$ . It follows that $\vec{u}$ is the projection of $\vec{a}$ onto the orthogonal complement of $T_{\vec{w}}P$ in $T_{\vec{w}}\mathbb{H}^{d}$ .

In an Euclidean vector space, the dot product of any vector with its projection onto any direction is positive unless the projection is zero. Since $T_{\vec{w}}\mathbb{H}^{d}$ is a space-like subspace, this statement applies to $\vec{a}$ and $\vec{u}$ . Note that $\vec{u}$ cannot be zero because otherwise $\vec{b}$ would belong to $V$ and hence $P$ . Thus, we conclude that $\vec{a}\cdot\vec{u}>0$ , which means that $\vec{u}$ points toward the side of $P_{b}$ . ∎

Thus, $\vec{u}$ is the tangent vector at $\vec{w}$ of the desired ray $\alpha^{+}$ . To compute points on this ray, we can use the exponential map at $\vec{w}$ .

Step 3: The exponential map

Finally, given a distance $d$ and a tangent vector $\vec{u}$ at $\vec{w}$ of a geodesic ray $\alpha^{+}$ in $\mathbb{H}^{d}$ , we need to compute the point on $\alpha^{+}$ that is of hyperbolic distance $d$ from $\vec{w}$ . This is based on the following lemma:

Lemma A.26.

Suppose that $\vec{u}$ is a unit tangent vector of $\mathbb{H}^{d}$ at $\vec{w}$ . Then

\alpha(t)=(\cosh t)\vec{w}+(\sinh t)\vec{u}

is a unit-speed geodesic with $\alpha(0)=\vec{w}$ and $\dot{\alpha}(0)=\vec{u}$ .

Proof.

Note that $\vec{w}$ and $\vec{u}$ are orthogonal by A.18. The facts that $\alpha(0)=\vec{w}$ , $\dot{\alpha}(0)=\vec{u}$ , and that $\alpha(t)$ is a unit-speed curve on $\mathbb{H}^{d}$ follow from simple, direct computations. To see that $\alpha(t)$ is a geodesic, note that it is the intersection of $\mathbb{H}^{d}$ with the linear subspace $\operatorname{span}(\vec{w},\vec{u})$ and apply A.20. ∎

It follows that

Proposition A.27.

If $\vec{u}$ is the tangent vector at the starting point $\vec{w}$ of a geodesic ray $\alpha^{+}$ in $\mathbb{H}^{d}$ then

(\cosh d)\vec{w}+(\sinh d)\frac{\vec{u}}{\sqrt{B(\vec{u},\vec{u})}}

is the point on $\alpha^{+}$ that is of distance $d$ from $\vec{w}$ . ∎

From the results in this section, we conclude that Algorithm 1 computes the horospherical projection $\pi^{\operatorname{H}}_{b,P}(x)$ .

Appendix B Geodesic Projection: Distortion Analysis

We show that in hyperbolic geometry, the geodesic projection $\pi^{\operatorname{G}}_{M}$ (also called closest-point projection) to a geodesic submanifold $M$ always strictly decreases distances unless the inputs are already in $M$ ; furthermore, all paths of distance at least $r$ from $M$ get shorter by at least $\cosh r$ times under the projection. Note that most ideas and results in this section already exist in the literature. For a comprehensive treatment, see [31].

This section is organized as follows. First, we use basic hyperbolic trigonometry to prove a bound on projected distances onto geodesics in dimension two (B.3). Then we generalize this bound to higher dimensions (Theorem B.4) and obtain an infinitesimal version of it (B.5). The main theorem (Theorem B.6) then follows immediately.

Hyperbolic trigonometry in two dimensions

We first recall the laws of sine and cosine for triangle in hyperbolic planes:

Proposition B.1.

Consider any triangle in $\mathbb{H}^{2}$ with edge lengths $A,B,C$ . Let $\alpha,\beta,\gamma$ be the angles at the vertices opposite to edges $A,B,C$ . Then

\frac{\sinh A}{\sin\alpha}=\frac{\sinh B}{\sin\beta}=\frac{\sinh C}{\sin\gamma},

\cos\gamma=\frac{\cosh A\cosh B-\cosh C}{\sinh A\sinh B}.

Proof.

See [31], section 2.6. ∎

Proposition B.2.

Consider any quadrilateral in $\mathbb{H}^{2}$ with vertices $x,y,z,w$ (in that order). Suppose that the angles at $z$ and $w$ are $90$ degrees. Then

	$\displaystyle\cosh(xy)=$	$\displaystyle\cosh(zw)\cosh(xw)\cosh(yz)$
		$\displaystyle-\sinh(xw)\sinh(yz).$

Here $\cosh(xy)$ is a shorthand for $\cosh(d_{\mathbb{H}}(x,y))$ .

Proof.

Applying the laws of sine to triangle $wxz$ gives

\frac{\sinh(xw)}{\sin\angle(zx,zw)}=\frac{\sinh(xz)}{\sin\angle(wx,wz)}=\sinh(xz),

\sin\angle(zx,zw)=\frac{\sinh(xw)}{\sinh(xz)}.

On the other hand, applying the law of cosine to triangle $xyz$ gives

\cos\angle(zx,zy)=\frac{\cosh(zx)\cosh(zy)-\cosh(xy)}{\sinh(zx)\sinh(zy)}.

Now note that $\sin\angle(zx,zw)=\cos\angle(zx,zy)$ . Thus,

\frac{\sinh(xw)}{\sinh(xz)}=\frac{\cosh(zx)\cosh(zy)-\cosh(xy)}{\sinh(zx)\sinh(zy)},

which can be simplified to

\cosh(zx)\cosh(zy)-\cosh(xy)=\sinh(xw)\sinh(zy).

Finally, note that the law of cosine for triangle $wxz$ gives

\frac{\cosh(wx)\cosh(wz)-\cosh(xz)}{\sinh(wx)\sinh(wz)}=\cos\angle(wx,wz)=0,

so $\cosh(xz)=\cosh(wx)\cosh(wz)$ . Plugging this into the above equation gives

		$\displaystyle\cosh(wx)\cosh(wz)\cosh(zy)-\cosh(xy)$
	$\displaystyle=$	$\displaystyle\sinh(xw)\sinh(zy),$

which can be arranged to get the desired identity. ∎

Now we are ready to prove a bound on the projected distances onto a geodesic in $\mathbb{H}^{2}$ :

Proposition B.3.

Let $L$ be a geodesic in $\mathbb{H}^{2}$ . Consider any $x,y\in\mathbb{H}^{2}$ that are on the same side of $L$ and are of distances at least $r$ from $L$ . Let $d=d_{\mathbb{H}}(x,y)$ and $d^{\prime}=d_{\mathbb{H}}(\pi^{\operatorname{G}}_{L}(x),\pi^{\operatorname{G}}_{L}(y))$ denote the original and projected distances of $x$ and $y$ . Then,

\sinh(d^{\prime}/2)\leq\frac{\sinh(d/2)}{\cosh r}.

Proof.

Let $r_{x}=d_{\mathbb{H}}(x,\pi^{\operatorname{G}}_{L}(x))$ and $r_{y}=d_{\mathbb{H}}(y,\pi^{\operatorname{G}}_{L}(y))$ . Then by B.2,

\cosh d=\cosh d^{\prime}\cosh r_{x}\cosh r_{y}-\sinh r_{x}\sinh r_{y}.

Since $\cosh a\cosh b-\sinh a\sinh b=\cosh(a-b)$ , this gives

\cosh d=(\cosh d^{\prime}-1)\cosh r_{x}\cosh r_{y}+\cosh(r_{x}-r_{y}).

Subtracting $1$ from both sides and using $\cosh(a)-1=2\sinh^{2}(a/2)$ gives

		$\displaystyle 2\sinh^{2}(d/2)$
	$\displaystyle=$	$\displaystyle 2\sinh^{2}(d^{\prime}/2)\cosh r_{x}\cosh r_{y}+2\sinh^{2}((r_{x}-r_{y})/2).$

Thus,

\sinh^{2}(d/2)\geq\sinh^{2}(d^{\prime}/2)\cosh r_{x}\cosh r_{y}.

Since $r_{x},r_{y}\geq r$ , this implies

\sinh^{2}(d/2)\geq\sinh^{2}(d^{\prime}/2)\cosh^{2}r,

which is equivalent to the desired inequality. ∎

General dimensions

Now we show that B.3 generalizes to higher dimensions.

Theorem B.4.

Let $M\subset\mathbb{H}^{d}$ be a totally geodesic submanifold. Consider any $x,y\in\mathbb{H}^{d}$ that are of distances at least $r$ from $M$ . Let $d=d_{\mathbb{H}}(x,y)$ and $d^{\prime}=d_{\mathbb{H}}(\pi^{\operatorname{G}}_{M}(x),\pi^{\operatorname{G}}_{M}(y))$ denote the original and projected distances of $x$ and $y$ . Then,

\sinh(d^{\prime}/2)\leq\frac{\sinh(d/2)}{\cosh r}.

Proof.

This is clearly true if $\pi^{\operatorname{G}}_{M}(x)=\pi^{\operatorname{G}}_{M}(y)$ or if both $x$ and $y$ belong to $M$ . Thus, without loss of generality, we can assume that $\pi^{\operatorname{G}}_{M}(x)\neq\pi^{\operatorname{G}}_{M}(y)$ and $\pi^{\operatorname{G}}_{M}(x)\neq x$ .

Let $L\subset M$ be the geodesic joining $\pi^{\operatorname{G}}_{M}(x)$ and $\pi^{\operatorname{G}}_{M}(y)$ , so that $d_{\mathbb{H}}(x,M)=d_{\mathbb{H}}(x,L)$ and $d_{\mathbb{H}}(y,M)=d_{\mathbb{H}}(y,L)$ . Let $L_{x}$ be the half-plane of $\operatorname{GH}(L\cup\{x\})$ that is bounded by $L$ and contains $x$ . Note that if we rotate $y$ around $L$ until it hits $L_{x}$ at a point $\bar{y}$ then

(i)

$\pi^{\operatorname{G}}_{L}(\bar{y})=\pi^{\operatorname{G}}_{L}(y)$ ,
(ii)

$d_{\mathbb{H}}(\bar{y},L)=d_{\mathbb{H}}(y,L)$ ,
(iii)

$\pi^{\operatorname{H}}_{x,L}(y)=\bar{y}$ by the “open book interpretation” of horospherical projections (see the beginning of Section A.2).

Thus, using (i), (ii) and applying B.3 to $L,x,\bar{y}$ gives,

\sinh(d^{\prime}/2)\leq\frac{\sinh(d_{\mathbb{H}}(x,\bar{y})/2)}{\cosh r}.

By (iii) and A.12,

d_{\mathbb{H}}(x,\bar{y})\leq d_{\mathbb{H}}(x,y).

Combining these two inequalities gives the desired inequality. ∎

Since $\sinh(t)\approx t$ as $t\to 0$ , applying Theorem B.4 as $y\to x$ gives

Corollary B.5.

Under the geodesic projection $\pi^{\operatorname{G}}_{M}(\cdot)$ , every tangent vector $\vec{v}$ at $x$ gets at least $\cosh(d_{\mathbb{H}}(x,M))$ times shorter. ∎

This implies the main theorem of this section:

Theorem B.6.

Let $\gamma(t)$ be any smooth path in $\mathbb{H}^{d}$ that stays a distance at least $r$ away from $M$ . Then

\operatorname{length}(\pi^{\operatorname{G}}_{M}(\gamma))\leq\frac{\operatorname{length}(\gamma)}{\cosh r}.

In particular, $\operatorname{length}(\pi^{\operatorname{G}}_{M}(\gamma))<\operatorname{length}(\gamma)$ unless $\gamma$ is entirely in $M$ .

Proof.

Note that the velocity vector of $\pi^{\operatorname{G}}_{M}(\gamma)$ is precisely $\pi^{\operatorname{G}}_{M}(\dot{\gamma})$ . Then, by B.5,

	$\displaystyle\operatorname{length}(\pi^{\operatorname{G}}_{M}(\gamma))$	$\displaystyle=\int\\|\pi^{\operatorname{G}}_{M}(\dot{\gamma}(t))\\|_{\mathbb{H}}dt$
		$\displaystyle\leq\int\frac{\\|\dot{\gamma}(t)\\|_{\mathbb{H}}}{\cosh r}dt=\frac{\operatorname{length}(\gamma)}{\cosh r}.$

Since $\cosh r\geq 1$ , this implies

\operatorname{length}(\pi^{\operatorname{G}}_{M}(\gamma))\leq\operatorname{length}(\gamma).

This can only be an equality if the $\cosh r$ factor is constantly $1$ , which means $\gamma$ stays in $M$ . ∎

When $\gamma$ is a geodesic, this gives a proof of 3.5: See 3.5

Remark B.7.

The above proposition does not imply that if $x$ and $y$ are of distances at least $r$ from $M$ then $d_{\mathbb{H}}(x,y)$ gets at least $\cosh r$ times shorter under the geodesic projection. This is because the geodesic segment between $x$ and $y$ can get closer to $M$ in the middle. This is inevitable when $x$ and $y$ are extremely far apart: For example, if $d_{\mathbb{H}}(x,M)=d_{\mathbb{H}}(y,M)=r$ , then the triangle inequality implies

d_{\mathbb{H}}(\pi^{\operatorname{G}}_{M}(x),\pi^{\operatorname{G}}_{M}(y))\geq d_{\mathbb{H}}(x,y)-2r.

Hence, if $d_{\mathbb{H}}(x,y)\gg r$ then the shrink factor cannot be much smaller than $1$ . Generally, the rule of thumb is that geodesic projections shrink distances by a big factor if the input points are not very far from each other but quite far from the target submanifold.

	Sarkar			PyTorch
dimension	$2$	$10$	$50$	$2$	$10$	$50$
Balanced Tree	0.088	0.136	0.152	0.097	0.022	0.022
Phylo Tree	0.639	0.639	0.638	0.292	0.286	0.284
Diseases	0.247	0.212	0.210	0.126	0.054	0.055
CS Ph.D.	0.388	0.384	0.380	0.189	0.140	0.139

Table 4: Average distortion for embeddings learned with different methods. Sarkar refers to combinatorial embeddings learned with Sarkar’s construction [28, 27] while PyTorch refers to embeddings learned using the PyTorch code from [13].

Appendix C Additional Experimental Details

We first provide more details about our experimental setup, including a description of how we implemented different methods and then provide additional experimental results.

	Balanced Tree		Phylo Tree		Diseases		CS Ph.D.
	distortion ( $\downarrow$ )	variance ( $\uparrow$ )	distortion ( $\downarrow$ )	variance ( $\uparrow$ )	distortion ( $\downarrow$ )	variance ( $\uparrow$ )	distortion ( $\downarrow$ )	variance ( $\uparrow$ )
PCA	0.91	1.09	0.99	0.28	0.95	0.75	0.96	1.25
tPCA	0.83	4.01	0.27	62.54	0.64	35.45	0.55	86.10
PGA	0.67 $\pm$ 0.013	11.15 $\pm$ 0.941	0.08 $\pm$ 0.001	1110.30 $\pm$ 0.577	0.55 $\pm$ 0.006	56.82 $\pm$ 0.493	0.25 $\pm$ 0.002	270.28 $\pm$ 0.487
BSA	0.58 $\pm$ 0.018	16.74 $\pm$ 1.278	0.08 $\pm$ 0.000	1110.05 $\pm$ 0.130	0.54 $\pm$ 0.008	57.73 $\pm$ 0.486	0.25 $\pm$ 0.001	270.98 $\pm$ 0.363
hMDS	0.12	58.62	0.08	498.05	0.32	568.22	0.18	1241.34
HoroPCA	0.06 $\pm$ 0.010	62.72 $\pm$ 0.924	0.01 $\pm$ 0.000	1190.27 $\pm$ 2.974	0.09 $\pm$ 0.010	161.32 $\pm$ 1.185	0.05 $\pm$ 0.004	994.69 $\pm$ 321.359

Table 5: Dimensionality reduction results on 10-dimensional combinatorial embeddings (Sarkar’s construction) reduced to two dimensions. Results are averaged over 5 runs for non-deterministic methods. Best in bold and second best underlined.

C.1 Datasets

Embedding

The datasets we use are structured as graphs with nodes and edges. We map these graphs to the hyperbolic space using different embedding methods. The datasets from [27] can be found in the open-source implementation⁷⁷7https://github.com/HazyResearch/hyperbolics and we compute hyperbolic embeddings for different dimensions using the PyTorch code from [13]. We also consider embeddings computed with Sarkar’s combinatorial construction [28] and report the average distortion measure with respect to the original graph distances in Table 4.

For classification experiments, we follow the experimental protocol in [8] and use the embeddings provided in the open-source implementation.⁸⁸8https://github.com/hhcho/hyplinear

Centering

For simplicity, we center the input embeddings so that they have a Fréchet mean of zero. This makes computations of projections more efficient. To do so, we compute the Fréchet mean of input hyperbolic points using gradient-descent, and then apply an isometric hyperbolic reflection that maps the mean to the origin.⁹⁹9This reflection can be computed in closed-form using circle inversions, see for instance Section 2 in [27].

C.2 Implementation Details in the Poincaré Model

Section A.4 (in particular, A.24 and Algorithm 1) describes how geodesic projections and horospherical projections can be computed efficiently in the hyperboloid model. However, because the Poincaré ball model is useful for visualizations and is popular in the ML literature, here we also give high-level descriptions of simple alternative methods to compute these projections directly in the Poincaré model.

This subsection is organized as follows. We first describe an implementation of geodesic projections in the Poincaré model (Section C.2.1). We then describe our implementation of all baseline methods, which rely on geodesic projections (Section C.2.2). Finally, we describe how HoroPCA could also be implemented in the Poincaré model (Section C.2.3).

C.2.1 Computing Geodesic Projections

Recall that geodesic projections map points to a target submanifold such that the projection of a point is the point on the submanifold that is closest to it. In the Poincaré model, when the target submanifold goes through the origin, geodesic projections can be computed efficiently as follows.

Consider any geodesic submanifold $P$ that contains the origin in the Poincaré Ball: this must be a linear space since any geodesic that goes through the origin in this model is a straight-line. The Euclidean reflection with respect to $P$ is a hyperbolic isometry (i.e. preserves distances). Given any point $x$ , we can compute its Euclidean (and hyperbolic) reflection $r_{P}(x)$ with respect to $P$ . Then, the (hyperbolic geometry) midpoint between $x$ and $r_{P}(x)$ belongs to $P$ and is the geodesic projection of $x$ onto $P$ . There is a closed-form expression for this midpoint, which can be derived using Mobïus operations.

C.2.2 Implementation of Baselines

We now detail our baseline implementation.

PCA and tPCA

We used the Singular Value Decomposition (SVD) PyTorch implementation to implement both PCA and tPCA. Before using the SVD algorithm, tPCA maps points to the tangent space at the Frćhet mean using the logarithmic map. Having centered the data, this mapping can be done using the logarithmic map at the origin which is simply:

\log_{\mathbf{o}}(x)=\operatorname{tanh^{-1}}(||x||)\frac{x}{||x||}.

PGA and BSA

Both PGA and BSA rely on geodesic projections. When the data is centered, the target submanifolds in BSA and PGA are simply linear spaces going through the origin. We can therefore compute the geodesic projections in these methods using the computations described above. To test their sensitivity to the base point, we also implemented two baselines that perturb the base point, by adding Gaussian noise.

hMDS

To implement hMDS, we simply implemented the formulas in the original paper (see Algorithm 2 in [27]).

hAE

For hyperbolic autoencoders, we used two fully-connected layers following the approach of [12]. One layer was used for encoding into the reduced number of dimensions, and one layer was used for decoding into the original number of dimensions. We trained these networks by minimizing the reconstruction error, and used the intermediate hidden representations as low-dimensional representations.

C.2.3 HoroPCA Implementation

The definition of horospherical projections onto $K>1$ directions involves taking the intersection of $K$ horospheres $S(p_{1},x)\cap\dots\cap S(p_{K},x)$ (Section 3.2.2). Although Section A.4 shows that horospherical projections can be computed efficiently in the hyperboloid model without actually computing these intersections, it is also possible to implement HoroPCA by directly computing these intersections. This can be done in a simple iterative fashion.

For example, in the Poincaré ball, horospheres can be represented as Euclidean spheres. Note that the intersection of two $d$ -dimensional Euclidean hyperspheres $S(p_{0},r_{0})$ and $S(p_{1},r_{1})$ is a $d-1$ -dimensional hypersphere which can be represented by a center, a radius, and $d-1$ -dimensional subspace. The radius is uniquely determined by the radii of the original spheres $r_{0},r_{1}$ and the distance between their centers $\|p_{0}-p_{1}\|$ and analyzing the 2-dimensional case, for which there is a simple closed formula. The center can similarly be found as a weighted combination of the centers $p_{0},p_{1}$ based on the formula for the 2D case. Lastly, this $d-1$ -dimensional subspace can be represented by noting it is the orthogonal space to the vector $p_{1}-p_{0}$ . Thus the intersection of these $d$ -dimensional hyperspheres can be easily computed, and this process can be iterated to find the intersection of $K$ spheres $S(p_{i},r_{i})$ .

C.3 Additional Experimental Results

C.3.1 Distortion Analysis

We first analyze the average distortion incurred by horospherical and geodesic projections on a toy example with synthetically-generated data. We generate 1000 points the Poincaré disk by sampling tangent vectors from a multivariate Gaussian, and mapping them to the disk using the exponential map at the origin. We then sample 100 random ideal points (directions) in the disk, and consider the corresponding straight-line geodesics pointing towards these directions. We project the datapoints onto these geodesics and measure the average distortion from before and after projection and visualize the results in a histogram (Fig. 6). As we can see, Horospherical projections achieve much lower distortion than geodesic projections on average, suggesting that these projections may better preserve information such as distances in higher-dimensional datasets.

C.3.2 Dimensionality Reduction Results

Sarkar embeddings

Table 2 showed dimensionality reduction results for embeddings learned with optimization. Here, we consider the same reduction experiment on combinatorial embeddings and report the results in Table 5. The results confirm the trends observed in Table 2: HoroPCA outperforms baseline methods, with significant improvements on distance preservation.

More dimension/component configurations

We consider the reduction of 50-dimensional PyTorch embeddings of the Diseases and CS Ph.D. datasets. We plot average distortion and explained Fréchet variance for different number of components in Fig. 7. HoroPCA significantly outperforms all previous generalizations of PCA. HoroPCA also outperforms hMDS, which is a competitive baseline, but not a PCA method as discussed before.

HoroPCA: Hyperbolic Dimensionality Reduction via Horospherical Projections

Abstract

1 Introduction

2 Background

2.1 The Poincaré Model of Hyperbolic Space

Geodesics

Geodesic submanifolds

2.2 Directions in Hyperbolic space

Ideal points

Busemann coordinates

Horospheres

2.3 Geodesic Projections

2.4 Manifold Statistics

3 Generalizing PCA to the Hyperbolic Space

3.1 Hyperbolic Flags

Definition 3.1.

Remark 3.2.

Remark 3.3.

3.2 Projections via Horospheres

3.2.1 Projecting onto K=1K=1 Directions

Proposition 3.4.

Proof.

3.2.2 Projecting onto K>1K>1 Directions

Proposition 3.5.

Computation

3.3 Intrinsic Variance Objective

4 HoroPCA

Euclidean PCA

HoroPCA

Base point independence

Theorem 4.1.

5 Experiments

5.1 Experimental Setup

Baselines

Datasets

Evaluation metrics

5.2 Dimensionality Reduction

5.3 Hyperbolic Data Whitening

5.4 Visualizations

6 Related Work

PCA methods in Riemannian manifolds

Dimension reduction in hyperbolic geometry

From hyperbolic to Euclidean

7 Conclusion

Acknowledgements

References

Appendix A Horospherical Projection: Proofs and Discussions

A.1 Well-definedness

Lemma A.1.

Proof.

Corollary A.2.

Proof.

Corollary A.3.

Proof.

Theorem A.4.

Proof.

A.2 Geometric properties

Theorem A.5 (The “open book” interpretation).

Corollary A.6.

Proof.

Corollary A.7.

Proof.

Remark A.8.

Corollary A.9.

Proof.

Corollary A.10.

Proof.

Corollary A.11.

Proof.

Corollary A.12.

Proof.

A.3 Detour: A Review of the Hyperboloid Model

Remark A.13.

Minkowski spaces

Types of vectors in ℝ1,d\mathbb{R}^{1,d}

Proposition A.14.

Proposition A.15.

The hyperboloid model of hyperbolic spaces

Definition A.16.

Remark A.17.

HoroPCA: Hyperbolic Dimensionality Reduction via
Horospherical Projections

3.2.1 Projecting onto $K=1$ Directions

3.2.2 Projecting onto $K>1$ Directions

Types of vectors in $\mathbb{R}^{1,d}$

Proposition A.24 (Geodesic projections in $\mathbb{H}^{d}$ ).