Year: 2020

Abstract:
We use hyperbolic Poisson kernel to construct the horocycle neuron model on hyperbolic spaces, which is a spectral generalization of the classical neuron model. We prove a universal approximation theorem for horocycle neurons. As a corollary, we obtain a state-of-the-art result on the expressivity of f a,p, which is used in the hyperbolic multiple linear regression. Our experiments get state-of-the-art results on the Poincare-embedding subtree classification task and the classification accuracy of the two-dimensional visualization of images.

heading: 1 INTRODUCTION
text: Conventional deep network techniques attempt to use architecture based on compositions of simple functions to learn representations of Euclidean data (LeCun et al., 2015). They have achieved remarkable successes in a wide range of applications (Hinton et al., 2012; He et al., 2016). Geometric deep learning, a niche field that has caught the attention of many authors, attempts to generalize conventional learning techniques to non-Euclidean spaces (Bronstein et al., 2017; Monti et al., 2017).
There has been growing interest in using hyperbolic spaces in machine learning tasks because they are well-suited for tree-like data representation (Ontrup & Ritter, 2005; Alanis-Lobato et al., 2016; Nickel & Kiela, 2017; Chamberlain et al., 2018; Nickel & Kiela, 2018; Sala et al., 2018; Ganea et al., 2018b; Tifrea et al., 2019; Chami et al., 2019; Liu et al., 2019; Balazevic et al., 2019; Yu & Sa, 2019; Gulcehre et al., 2019; Law et al., 2019). Many authors have introduced hyperbolic analogs of classical learning tools (Ganea et al., 2018a; Cho et al., 2019; Nagano et al., 2019; Grattarola et al., 2019; Mathieu et al., 2019; Ovinnikov, 2020; Khrulkov et al., 2020; Shimizu et al., 2020).
Spectral methods are successful in machine learning, from nonlinear dimensionality reduction (Belkin & Partha, 2002) to clustering (Shi & Malik, 2000; Ng et al., 2002) to hashing (Weiss et al., 2009) to graph CNNs (Bruna et al., 2014) to spherical CNNs (Cohen et al., 2018) and to inference networks (Pfau et al., 2019). Spectral methods have been applied to learning tasks on spheres (Cohen et al., 2018) and graphs (Bruna et al., 2014), but not yet on hyperbolic spaces. This paper studies a spectral generalization of the FC (affine) layer on hyperbolic spaces.
Before presenting the spectral generalization of the affine layer, we introduce some notations. Let (·, ·)E be the inner product, | · | the Euclidean norm, and ρ an activation function. The Poincaré ball model of the hyperbolic space Hn(n≥2) is a manifold {x∈Rn : |x|<1} equipped with a Riemannian metric ds2Hn= ∑n i=1 4(1−|x|
2)−2dx2i . The boundary of Hn under its canonical embedding in Rn is the unit sphere Sn−1. The classical neuron y=ρ((x,w)E+b) is of input x∈Rn, output y∈R, with trainable parameters w∈Rn, b∈R. An affine layer Rn → Rm is a concatenation of m neurons. An alternative representation of the neuron x 7→ρ((x,w)E+b) is given by 1
x∈Rn 7→ ρ(λ(x, ω)E+b), ω∈Sn−1, λ, b∈R. (1) This neuron is constant over any hyperplane that is perpendicular to a fixed direction ω. In Hn, a horocycle is a n−1 dimensional sphere (one point deleted) that is tangential to Sn−1. Horocycles are hyperbolic counterparts of hyperplanes (Bonola, 2012). Horocyclic waves 〈x, ω〉H := 12 log 1−|x|2 |x−ω|2 are constant over any horocycle that is tangential to Sn−1 at ω. Therefore,
x∈Hn 7→ ρ(λ〈x, ω〉H+b), ω∈Sn−1, λ, b∈R (2) 1if w 6= (0, . . . , 0), one can take ω = w/|w|, λ = |w|; else, one can take λ = 0 and any ω ∈ Sn−1.
generalizes the classical neuron model (1), and a concatenation of finitely many (2) generalizes the FC (affine) layer. We call (2) a horocycle neuron. Figure 1 (middle) is an example on H2.
The neuron models in (1, 2) are related to spectral theory because (·, ω)E (respectively 〈·, ω〉H ) are building blocks of the Euclidean (respectively hyperbolic) Laplacian eigenspace. Moreover, many L2 spaces have a basis given by Laplacian eigenfunctions (Einsiedler & Ward, 2017). On one side, all Euclidean (respectively hyperbolic) eigenfunctions are some kind of “superposition” of (·, ω)E (respectively 〈·, ω〉H ). On the other side, neural networks based on (1) (respectively (2)) represent functions that are another kind of “superposition” of (·, ω)E (respectively 〈·, ω〉H ). They heuristically explain why the universal approximation property is likely to hold for networks constructed by (1) and (2). By using the Hahn Banach theorem, an injectivity theorem of Helgason, and integral formula, we prove that finite sums of horocycle neurons (2) are universal approximators (Theorem 2).
Let p ∈ Hn, Tp(Hn) be the tangent space of Hn at p, a ∈ Tp(Hn), ⊕ be the Möbius addition (Ungar, 2008). We remind the reader that the following functions
f1a,p(x) = 2|a| 1− |p|2 sinh −1 ( 2(−p⊕ x, a)E (1− | − p⊕ x|2)|a| ) (3)
are building blocks of many hyperbolic learning tools (Ganea et al., 2018a; Mathieu et al., 2019; Shimizu et al., 2020). Figure 1 illustrates examples of different neuron models (1, 2, 3) on H2.
In Lemma 1, we shall present a close relationship between (2) and (3). Using this relationship and Theorem 2, we obtain a novel result on the expressivity of f1a,p (Corollary 1).
This article contributes to hyperbolic learning. We first apply spectral methods, such as the horocycle, to hyperbolic deep learning. We prove results on the expressivity of horocycle neurons (2) and f1a,p (3). With horocycle neurons, we obtain state-of-the-art results on the Poincaré-embedding subtree classification task and the classification accuracy of the 2-D visualization of images in in the experiment.

heading: 2 RELATED WORK
text: Universal approximation There is a vast literature on universal approximation (Cybenko, 1989; Hornik et al., 1989; Funahashi, 1989; Leshno et al., 1993). Cybenko (1989)’s existential approach uses the Hahn Banach theorem and Fourier transform of Radon measures. To prove Theorem 2, we also use the Hahn Banach theorem, and additionally an integral formula (7) and an injectivity Theorem 1 of Helgason. Generalizing integral formulas and injectivity theorems is easier than generalizing Fourier transform of Radon measures on most non-Euclidean spaces. (Carroll & Dickinson, 1989) uses the inverse Radon transform to prove universal approximation theorems. This method relates to ours, as injectivity theorems are akin to inverse Radon transforms. However, using the injectivity theorem is an existential approach while using the inverse Radon transform is a constructive one.
Spectral methods Spectral methods in Bronstein et al. (2017); Bruna et al. (2014); Cohen et al. (2018) use a basis of L2(X) given by eigenfunctions, whereX is a finite graph or the sphere. Because L2(Hn) has no eigenfunctions as a basis, our approach is different from theirs.
Hyperbolic deep learning One part of hyperbolic learning concerns embedding data into the hyperbolic space (Nickel & Kiela, 2017; Sala et al., 2018). Another part concerns learning architectures with hyperbolic data as the input (Ganea et al. (2018a); Cho et al. (2019)). Ganea et al. (2018a) proposes two ways to generalize the affine layer on hyperbolic spaces: one by replacing the linear and bias part of an affine map with (25, 26) of their paper; another one by using a concatenation of f1a,p in
their hyperbolic multiple linear regression (MLR). The latter seems more relevant to ours. A level set of f1a,p is a hypercycle that has the same distance to a chosen geodesic hypersurface, while a level set of a horocycle neuron is a horocycle that has the same “spectral” distance to an ideal point at infinity. Based on functions similar to f1a,p, Mathieu et al. (2019); Shimizu et al. (2020) build the gyroplane layer and Poincaré FC layer. Ganea et al. (2018a); Cho et al. (2019) take geodesics as decision hyperplanes, while we (initially) take horocycles. We shall construct the horocycle multiple linear regression (MLR), where decision hypersurfaces are geodesics. Geodesics decision hyperplanes (Ganea et al., 2018a; Cho et al., 2019) and geodesic decision hypersurfaces here arise from different methods. Khrulkov et al. (2020) investigates hyperbolic image embedding, where prototypes (or models) of each class are center-based. We study a different one, and we shall call our prototypes end-based.

heading: 3 HYPERBOLIC SPACES
text: This section reviews facts from hyperbolic geometry that are used in the proof of Theorem 2. For the reader who is not interested in the proof, (4) is enough for the implementation.
Hyperbolic metric We use the Poincaré model. The hyperbolic space Hn(n≥2) is the manifold {x∈Rn : |x|<1} equipped with a Riemannian metric ds2 = ∑n i=1 4(1−|x|2)−2dx2i . Let o be the origin of Hn. The distance function dHn satisfies dHn(o, x)=2 arctanh |x|.
Geodesics, horocycles and corresponding points Geodesics in Hn are precisely circular arcs that are orthogonal to Sn−1. Horocycles in Hn are precisely (n−1)-dimensional spheres that are tangential to Sn−1 (Helgason, 1970). Horocycles are hyperbolic analogs of hyperplanes. Figure 2 illustrates geodesics and horocycles on H2. Hyperbolic Poisson kernel The Poisson kernel for Hn is P (x, ω)= (
1−|x|2 |x−ω|2
)n−1 , where
x∈Hn, ω∈Sn−1 (Helgason (1970)[p.108]). The function 〈·, ω〉H defined by
〈x, ω〉H = 1 2(n− 1) logP (·, ω) = 1 2 log
1− |x|2 |x− ω|2 (4)
is constant over any horocycle that is tangential to Sn−1 at ω (Figure 1 (middle), (6)).
Riemannian volume The Riemannian volume induced by the metric ds2 on Hn is dVol = 2n(1− |x|2)−ndx1 . . . dxn. (5)
Horocycles Let Ξ be the set of horocycles of Hn, and let Ξω be the set of all horocycles that are tangential to Sn−1 at ω. Given λ∈R, we let ξλ,ω be the unique horocycle that connects ω and tanh (λ/2) · ω. We have Ξω = ∪λ∈R{ξλ,ω} and Ξ = ∪ω∈Sn−1Ξω. The length of any geodesic (that ends at ω) line segment cut by ξλ1,ω and ξλ2,ω equals |λ1 − λ2| (A.2). Therefore |λ1 − λ2| is a natural distance function defined on Ξω, and the map λ→ ξλ,ω is an isometry between R and Ξω. This isometry is closely related to 〈·, ω〉H (A.3): for any x ∈ ξλ,ω , 〈x, ω〉H = λ/2. (6) The annoying /2 in (6) is a tradeoff that the metric here is different from that in Helgason (2000).
Integral formula For fixed ω ∈ Sn−1, Hn=∪λ∈Rξλ,ω. Let dVolξλ,ω be the measure induced by ds2 on ξλ,ω . Let L be a family of geodesics that end at ω, δ > 0, and U=L ∩ (∪λ≤α≤λ+δξα,ω). For l ∈ L, dH(l ∩ ξλ,ω, l ∩ ξλ+δ,ω)=δ (A.2), hence dVol(U) = δ · dVolξλ,ω (U ∩ ξλ,ω) and therefore∫
Hn f(x)dVol(x) = ∫ R (∫ ξλ,ω f(z)dVolξλ,ω (z) ) dλ. (7)
The above proof (for Hn) is essentially the same as that in (Helgason, 2000)[p.37] (for H2). To further convince the reader that (7) holds for all n, we give another simple proof in A.4.
Injectivity theorem With respect to the canonical measure on Ξ, Helgason (1970)[p.13] proved Theorem 1 (Helgason). If f ∈ L1(Hn) and ∫ ξ f(z)dVolξ(z) = 0 for a.e ξ ∈ Ξ, then f = 0 a.e..
Theorem 1 demonstrates that if the integral of f ∈ L1(Hn) over almost every horocycle is zero then f is also zero. This theorem and the integral formula (7) are essential for the proof of Theorem 2.

heading: 4 LEARNING ARCHITECTURES AND EIGENFUNCTIONS OF THE LAPLACIAN
text: In this section, we discuss a heuristic connection between the representation properties of eigenfunctions and classical neurons, and then we define some horocycle-related learning tools.

heading: 4.1 EIGENSPACES AND NEURON MODELS
text: On a Riemannian manifold X , the Laplace-Beltrami LX is the divergence of the gradient, and it has a well-known representation property (Einsiedler & Ward, 2017): if X is a compact Riemannian manifold or bounded domain in Rn, then L2(X) has a basis given by eigenfunctions. This statement is false if X is Rn or Hn (Hislop, 1994).
Eigenspaces of on Rn and Hn Our work is motivated by the theory of eigenspaces, in which Euclidean (respectively hyperbolic) eigenfunctions are obtained from (x, ω)E (respectively 〈x, ω〉H ) by some kind of superposition. For example, all smooth eigenfunctions of LRn are precisely the functions (M. Hashizume & Okamoto, 1972)[p.543]
f(x) = ∫ Sn−1 eλ(x,ω)EdT (ω), (8)
and eigenfunctions of LHn are precisely the functions (Helgason, 1970)[Theorem 1.7, p.139]
f(x) = ∫ Sn−1 eλ〈x,ω〉HdT (ω), (9)
where T in (8) and (9) are some technical linear forms of suitable functional spaces on Sn−1.
Neuron models By (8) and (1), Euclidean eigenfunctions (respectively classical neurons) are superpositions of (·, ω)E and exp (respectively ρ), with homogeneity and additivity. By (9) and (2), hyperbolic eigenfunctions (respectively horocycle neurons) are superpositions of 〈·, ω〉H and exp (respectively ρ). The representation property of eigenfunctions on compact manifolds and bounded domains suggests that the universal approximation property is likely to hold for networks constructed by (·, ω)E or 〈·, ω〉H . However, this heuristic is not proof (A.5).

heading: 4.2 HOROCYCLE BASED LEARNING ARCHITECTURES
text: Horocycle neuron In the implementation of the horocycle neuron (2), we take 1 2
log (
1−|x|2 |x−ω|2+ + ) for 〈x, ω〉H , where is a small constant to ensure numerical stability. For updating ω, we use the sphere optimization algorithm (Absil et al., 2008; Bonnabel, 2013) (A.6).
Horocycle feature and horocycle decision hypersurface Given a non-origin point x ∈ Hn, for y ∈ Hn we define hx(y) = 〈y, x/|x|〉H and call it the horocycle feature attached to x. This feature is useful in the Poincaré embedding subtree classification task (see the experiment and Figure 3[left]). The horocycle is the hyperbolic analog of the Euclidean hyperplane, and therefore it could be a possible choice of decision hypersurface, which may arise from a level set of a horocycle feature.
End-based clustering and end prototype Natural clustering is a topic in representation learning (Bengio et al., 2013), and the common prototype-based clusters are center-based (Tan et al., 2005). We propose a type of clustering that embeds high-dimensional data in Hn and places prototypes in Sn−1. Figure 3[right] is an example for n = 2. For ω ∈ Sn−1 and any b ∈ R, the function x ∈ Hn 7→ − log ( 1−|x|2 |x−ω|2 ) + b measures the relative distance of Hn from ω in Gromov’s bordification theory (Bridson & Haefliger (2009)[II.8], A.18). Moreover, we define Dist : Hn ×Sn−1 ×R→ R by
Dist(x, ω, b) = − log ( 1− |x|2
|x− ω|2
) + b = −2〈x, ω〉H + b. (10)
It is a relative distance function, and this is why Dist may assume negative values and why there is a bias term b in (10). Consider classes Cls = {C1, C2, . . . , CM} and labeled training examples {(X1, Y 1), . . . , (XN , Y N )}, where Xi ∈ RD are D-dimensional input features and Y i ∈ {1, 2, . . . ,M}. Each example Xi belongs to the class CY i . In light of (10), our goal is to find a neural network NNθ : RD → Hn that is parameterized by θ, prototypes ω1, . . . , ωM ∈ Sn−1, and real numbers b1, . . . , bM ∈ R such that
# { 1≤i≤N : Y i = arg min
1≤j≤M
( Dist(NNθ(X i), ωj , bj) )}
N (11)
is maximized. We call {NNθ(Xj) : 1 ≤ j ≤ N} the end-based clustering and ωi end prototypes (in hyperbolic geometry, the end is an equivalence class of parallel lines in Figure 2[left]). In experiments, we take NNθ = Exp ◦ NN′θ, where NN ′ θ : R
D → Rn is a standard neural network parameterized by θ and Exp : Rn → Hn is the exponential map of the hyperbolic space.
Horocycle layer, horocycle multiple linear regression (MLR) and geodesic decision hypersurfaces We call a concatenation of (2) a horocycle layer, and we shall carefully describe a prototypical learning framework for end-based clusterings. Using the same notions as in the previous paragraph, the classification task has M classes, and NNθ = Exp ◦NN′θ : RD → Hn is a deep network. For prototypes ω1, . . . , ωM ∈ Sn−1, real numbers b1, . . . , bM ∈ R, and any exampleX , our feedforward for prediction will be
x = NNθ(X), (Feature descriptor) SCj(X) = −Dist(x, ωj , bj), (Scores; Similarity)
X ∈ Cargmax 1≤j≤M (SCj(X)). (Classifier)
The goal is to maximize the accuracy (11), and then we need a loss function for the backpropagation. Following the convention of prototypical networks (Snell et al., 2017; Yang et al., 2018), we choose an increasing function ρ (in our experiments, ρ(x) = x or ρ = tanh. 2) and let the distribution over classes for an input X (with label Y ) be
pθ(Y = Cj |X) ∝ e−ρ(Dist(NNθ(X),mj ,bj)) = e−ρ(−SCj(X)). 2One often takes ρ(x) = x2 in metric learning, which is improper here because Dist(x) could be negative.
Therefore, given a batch of training examples, the loss function is L = − ∑
(Xj ,Y j)∈Batch log pθ(Y = CY j |Xj) #Batch . (12)
The training proceeds by minimizing L, and we call this framework a horocycle MLR. The set of parameters of the framework is {θ} ∪ {ω1, . . . , ωM} ∪ {b1, . . . , bM}. It is worth mentioning that decision boundaries of the horocycle MLR are geodesics, which follows from
SCi(X)=SCj(X)⇐⇒ log ( 1−|x|2
|x−ωi|2
) −bi = log ( 1−|x|2
|x−ωj |2
) −bj ⇐⇒
|x−ωi| |x−ωj | = e bj−bi 2
and the theorem of Apollonian circles (A.7).
Poisson neuron and Poisson multiple linear regression (MLR) Although 〈x, ω〉H (4) is wellmotivated by the theory of eigenspaces (9) and fits naturally into metric learning (see 10 or also Corollary 1), it is only defined on Hn. Some readers might not be convinced that the neuron has to be defined on hyperbolic spaces. Therefore, we try to remove the log in (4) and define the Poisson neuron model by Pρw,λ,b(x) = ρ ( λ |w| 2−|x|2 |x−w|2 + b ) for w ∈ Rn, λ, b ∈ R, which is well-defined on Rn\{w}. Notice that if |x| < |w| then |w| 2−|x|2 |x−w|2 = e 2〈x/|w|,w/|w|〉H . In A.8, Figure 7 illustrates an example of a Poisson neuron on R2. In the implementation, we take |w| 2−|x|2
|x−w|2+ for |w|2−|x|2 |x−w|2 , where
is a small constant for numerical stability. We call a concatenation of Poisson neurons a Poisson layer, and we use it with a deep neural network NNθ : RD → Rn to construct the Poisson MLR, which is similar to the horocycle MLR. Let w1, . . . , wM ∈ Rn and b1, . . . , bM ∈ R, the feedforward for prediction of our framework is
x = NNθ(X), SCj(X) = BatchNorm(P ρ wj ,−1,bj (x)), X ∈ Cargmax 1≤j≤M (SCj(X)). (13)
We let the pθ(Y = Cj |X) ∝ eSCj(X) and take (12) as the loss. This framework is called a Poisson MLR. We use the usual optimization algorithms to update parameters in the Poisson neuron. The BatchNorm(Ioffe & Szegedy, 2015) seems crucial for (13) in the experiment. Figure 4 illustrates that high-confidence prediction regions (deep red areas) of the Poisson MLR are compact sets, in contrast to classical classifiers Hein et al. (2019)[Theorem 3.1]. We shall use this figure to explain an experiment in Section 6.4.

heading: 5 REPRESENTATIONAL POWER
text: In this section, ρ is a continuous sigmoidal function (Cybenko, 1989), ReLU(Nair & Hinton, 2010), ELU(Clevert et al., 2016), or Softplus(Dugas et al., 2001). We remind the reader that ρ is sigmoidal if lim t→∞ ρ(t) = 1 and lim t→−∞ ρ(t) = 0. The following theorem justifies the representational power of horocycle neurons. Theorem 2. Let K be a compact set in Hn, and 1≤p<∞. Then finite sums of the form
F (x) = N∑ i=1 αiρ(λi〈x, ωi〉H+bi), ωi∈Sn−1, αi, λi, bi∈R (14)
are dense in Lp(K,µ), where µ is either dVol (5) or the induced Euclidean volume.
We provide a sketch of the proof here and go through the details in A.9. It suffices to prove the theorem for a sigmoidal function ρ and µ = dVol , as other cases follow from this one. Assume that these finite sums are not dense in Lp(K, dVol). By the Hahn-Banach theorem, there exists some nonzero h∈Lq(K, dVol), where q=p/(p− 1) if p>1 and q=∞ if p=1, such that ∫ K F (x)h(x)dVol(x) = 0 for all finite sums of the form (14). Extend h to be a function H that is defined on Hn by assigning H(x)=h(x) if x∈K and H(x)=0 if x∈Hn\K. Using the property of sigmoidal functions, the bounded convergence theorem, and the integral formula (7), we prove that the integration of H on almost every horocycle is zero. By the injectivity Theorem 1, H is almost everywhere zero, which contradicts our assumption and completes the proof.
In A.10, we shall prove the same result for Poisson neurons. In A.11, we prove the following lemma, which demonstrates a close relationship between horocycle neurons and the widely used f1a,p (3).
Lemma 1. Let K be a compact set in Hn, ω ∈ Sn−1, and > 0. There are c, d ∈ R, p ∈ Hn, and a ∈ Tp(Hn) such that the function D(x) = cf1a,p(x) + d− 〈x, ω〉H satisfies ||D||Lp(K,dVol) < .
This lemma suggests that 〈·, ω〉H is a boundary point of some “compactification” of the space of f1a,p. The above lemma together with Theorem 2 implies
Corollary 1. Let K be a compact set in Hn and 1≤p<∞. Finite sums of the form
F (x) = N∑ i=1 αiρ(cif 1 ai,pi(x) + di), pi ∈ H n, ai ∈ Tpi(H n), αi, ci, di ∈ R,
are dense in Lp(K,µ), where µ = dVol or µ is the induced Euclidean volume.
This result provides novel insights into the hyperbolic neural network (Ganea et al., 2018a), gyroplane layer (Mathieu et al., 2019), and Poincaré FC layer (Shimizu et al., 2020). Although level sets of f1a,p are hypercycles, our proof of Lemma 1 relies on the theory of horocycles. It would be interesting to have more natural approaches to treat the expressivity of f1a,p.

heading: 6 EXPERIMENTS
text: In this section, we first play with the MNIST toy. Next, we apply a horocycle feature to the Poincaré embedding subtree classification task. After that, we construct 2-D clusterings of image datasets by using the horocycle MLR. Finally, we provide evidence for further possible applications of the Poisson MLR. We use the framework or some functions of Tensorflow, Keras, and scikit-learn (Abadi et al., 2015; Chollet et al., 2015; Pedregosa et al., 2011).

heading: 6.1 MNIST
text: The MNIST (LeCun et al., 1998) task is popular for testing hyperbolic learning tools (Ontrup & Ritter, 2005; Nagano et al., 2019; Mathieu et al., 2019; Grattarola et al., 2019; Ovinnikov, 2020; Khrulkov et al., 2020). We train two different classifiers. A.12, A.14, and code contain details. The first one is a single horocycle layer followed by the softmax classifier. The average test error rate after 600 epochs is 1.96%, and Theorem 2 provides the rationale for this experiment (A.13). The second one is a Poisson MLR. It is the best hyperbolic geometry related MNIST classifier (Table 1). In this table, Ontrup & Ritter (2005) uses the hyperbolic SOM, Grattarola et al. (2019) uses the adversarial autoencoder, and Khrulkov et al. (2020) uses the hyperbolic MLR. Our experiment performs well on MNIST suggests that horocycle and Poisson neurons are computationally efficient and easily coordinate with classical learning tools (such as the convolutional layer and the softmax).

heading: 6.2 POINCARÉ EMBEDDING SUBTREE CLASSIFICATION
text: Given a Poincaré embedding (Nickel & Kiela, 2017) PE : {WordNet noun} → HD of 82114 nouns and given a node x ∈ {WordNet noun}, the task is to classify all other nodes as being part of the subtree rooted at x (Ganea et al., 2018a). Our model is logistic regression, where the horocycle feature p ∈ {WordNet noun} 7→ hPE(x)(PE(p)/s) (s is a hyperparameter lying in [1, 1.5]) is the only predictor, and the dependent variable is whether p is in the subtree rooted at x. The decision hypersurface of this model is a horocycle, as illustrated in Figure 3 (left). In the experiment, we pre-train three different Poincaré embeddings3 in each of H2,H3,H5,H10. For each x ∈ {animal, group, location, mammal, worker} and D ∈ {2, 3, 5, 10}, we randomly select one of three pre-trained Poincaré embedding PE : {WordNet noun} → HD and then test the model. Table 2 reports the F1 classification scores and two standard deviations of 100 trials for each {x,D}. Different Poincaré embeddings account for the most variance of the performance. Our model is different from the existing ones. Firstly, we take the horocycle as the decision hypersurface, while others take the geodesic. Secondly, we train a logistic regression on top of the horocycle feature attached to PE(x), which is efficiently calculated, while others train the hyperbolic MLR with different parametrizations. On the number of parameters, we have three (independent of D), Ganea et al. (2018a) has 2D, and Shimizu et al. (2020) has D + 1. The number of parameters explains why our model is prominent in low dimensions.

heading: 6.3 END-BASED CLUSTERING FOR 2D DIMENSION REDUCTION
text: In this experiment, we use the horocycle MLR (Section 4.2) to construct end-based clusterings NNθ : R
D → H2 for MNIST, Fashion-MNIST(Xiao et al., 2017), and CIFAR-10(Krizhevsky, 2012). We take NNθ = Exp ◦ NN′θ, where Exp is the exponential map of H2 and NN ′ θ : R
D → R2 is a network with four convolutional blocks for MNIST/Fashion-MNIST or a ResNet-32 structure for CIFAR-10. A.16 and code contain details.
3https://github.com/dalab/hyperbolic_cones
Figure 5 illustrates end-based clusterings for MNIST, Fashion-MNIST, and CIFAR-10, with performance reported in the caption. Our accuracy for Fashion-MNIST is 8% higher than all numbers presented in McInnes et al. (2020). Moreover, Table 3 compares the numbers of Yang et al. (2018); Ghosh & Kirby (2020), and ours for MNIST, and our methods are similar. We all use convolutional networks as the (Feature descriptor) and prototype-based functions as the loss. However, Yang et al. (2018); Ghosh & Kirby (2020) use the center-based prototype loss, while we use the end-based (12). Yang et al. (2018)[Figure 1] points out that the traditional CNN is good at linearly separating feature representations, but the learned features are of large intra-class variations. The horocycle MLR leads to the inter-class separability in the same way (angle accounts for label difference) a traditional CNN does. At the same time, it also obtains intra-class compactness (Figure 5).

heading: 6.4 POISSON MLR
text: Using a Poisson MLR whose feature descriptor is a ResNet-32 structure, we obtain a classifier with a test error rate of 6.46% on CIFAR-10. It is on par with other methods with similar network structures (Yang et al., 2018). Moreover, we apply Poisson MLR to the classification task of flowers (Tensorflow), which is a typical example of overfitting. Replacing the MLR part of the Keras model (Tensorflow) with a Poisson MLR, the new Poisson model shows better generalization performance (Figure 6). A.17 and code contain the details. This subsection provides evidence for further applications of horocycles.

heading: 7 CONCLUSION
text: Based on the spectral theory of hyperbolic spaces, we introduce several horocycle-related learning tools. They find applications in the hyperbolic neural networks, the Poincaré embedding subtree classification task, and the visualization and classification of image datasets. We give an existential proof of a universal approximation theorem for shallow networks constructed by horocycle neurons or f1a,p. Hopefully, it will trigger further research on the expressivity problems, such as constructive approaches, quantitative results, and benefit of depth (Mhaskar & Poggio, 2016), on horocycle neurons, f1a,p, and similar functions on more general manifolds.

heading: A APPENDIX

heading: A.1 NOTATIONS AND SYMBOLS

heading: Default Notations
text: Notation Description Related formula
R The set of real numbers
Rn n dimensional Euclidean space x ∈ Rn, x = (x1, . . . , xn) (·, ·)E Euclidean inner product x ∈ Rn, y ∈ Rn, (x, y)E = ∑n i=1 xiyi 〈·, ·〉H Hyperbolic analogue of (·, ·)E x ∈ Hn, y ∈ Sn−1, 〈x, ω〉H = 12 log 1−|x|2 |x−ω|2
| · | Euclidean norm x ∈ Rn, |x| = √ (x, x)E
Hn n dimensional hyperbolic space as a set, Hn = {x ∈ Rn : |x| < 1} Tp(X) Tangent space of X at p
T (X) Tangent space of X T (X) = ∪p∈XTp(X) ds2Hn The canonical metric on Hn with curva-
ture -1 ds2Hn=
∑n i=1 4(1−|x| 2)−2dx2i
dVol Riemannian volume on Hn dVol = 2n(1− |x|2)−ndx1 . . . dxn Lp(K, dVol) Lp space Lp(K, dVol) = { f | ∫ K |f |pdVol <∞ } || · ||Lp(K,dVol) Lp norm f measurable on K, ||f ||Lp(K,dVol) = (∫ K |f |pdVol ) 1 p Sn−1 n− 1 dimensional sphere as a set, Sn−1 = {x ∈ Rn : |x| = 1}
P (·, ·) Hyperbolic Poisson kernel x ∈ Hn, ω ∈ Sn−1, P (x, ω) = (
1−|x|2 |x−ω|2 )n−1 f1a,p Model in the hyperbolic MLR f1a,p(x) = 2|a| 1−|p|2 sinh −1 ( 2(−p⊕x,a)E (1−|−p⊕x|2)|a|
) dHn The hyperbolic distance function
Ξ The space of horocycles
Ξω The set of horocycles that are tangential to Sn−1 at ω
LX Laplace-Beltrami operator on X
hx The horocycle feature function hx(y) = 〈y, x/|x|〉H ξλ,ω The unique horocycle connecting ω and tanhλ/2 · ω. MLR Multiple linear regression
dim dimension
IK the indicator function of K
Dist Relative distance function Dist(x, ω, b) = −2〈x, ω〉H + b Cls Set of classes Cls = {C1, C2, . . . , CM} NNθ A network parameterized by θ
NN′θ A network parameterized by θ
Exp Exponential map of the hyperbolic space
(X1, Y 1) Labeled sample
SCj Score function
pθ(Y = Cj |X) Prediction probability L Loss function
Pρw,λ,b Poisson neuron P ρ w,λ,b(x) = ρ
( λ |w|
2−|x|2 |x−w|2 + b ) PE Poincaré embedding

heading: Conventional symbols
text: Symbol In most cases it refers
n,m, i integers
x, y, w points in Rn or Hn, or real numbers
o the origin of Rn or Hn
b, c, d, α, δ real numbers
λ real or complex number
t real number, represent the timestamp in optimization
ω point in Sn−1
ρ an activation function
f, g functions
K a compact set
X a manifold
p a point in Hn or on a manifold
a an element in Tp(Hn)
ξ a horocycle
µ a measure
L a family of geodesics lines
l a geodesics line
U a set in Hn
F, h,H functions
M number of classes
D dimension

heading: A.2 PROOF OF THE ISOMETRY
text: Given ω∈Sn−1 and λ∈R, we let ξλ,ω the unique horocycle that connects ω and tanh (λ/2) · ω. The length of any geodesic (that ends at ω) line segment cut by ξλ1,ω and ξλ2,ω equals |λ1 − λ2|. This fact is obvious in the half-space model.
There is a Riemannian isometry F : {z ∈ Rn : |z| < 1} → {(x1, · · · , xn) : x1 > 0} (the latter is with the metric ds2 = dx 2 1+···+dx 2 n
x21 ) such that F (ω) = ∞ and F (o) = (1, 0, . . . , 0). Using
dHn(o, tanh(λi/2)ω) = |λi|, d{(x1,··· ,xn):x1>0}((1, 0, . . . , 0), (e±λi , 0, . . . , 0)) = |λi|, F (ω) =∞ and F (o) = (1, 0, . . . , 0), we have F (tanh(λi/2)ω) = (eλi , 0, . . . , 0). Therefore, F maps ξλi,ω to {(x1, x2, . . . , xn) : x1 = eλi}. Any geodesic (that ends at ω) line segment cut by ξλ1,ω and ξλ2,ω is mapped by F to {(t, α2, . . . , αn) : (t− eλ1)(t− eλ2) < 0} for some fixed αj . It is easy to check the length of this segment with respect to dx 2 1+···+dx 2 n
x21 (as αi are constants, the metric reduces to dx21/x 2 1
on this segment) is |λ1 − λ2|.

heading: A.3 PROOF OF (6)
text: Because x is on ξλ which is a sphere with center 1+tanhλ/2 2 ω and radius 1−tanhλ/2 2 , we have∣∣∣x− 1+tanhλ/22 ω∣∣∣2 = ∣∣∣ 1−tanhλ/22 ∣∣∣2, which leads to |x|2−(1+tanhλ/2)(x, ω)E+tanhλ/2|ω|2 =
0, and then 1+tanhλ/22 |x− ω| 2 = 1−tanhλ/22 (|ω 2| − |x|2), and finally 〈x, ω〉H = 12 log |ω|2−|x|2 |x−ω|2 = 1 2 log 1+tanhλ/2 1−tanhλ/2 = λ/2.

heading: A.4 ANOTHER PROOF OF THE INTEGRAL FORMULA (7)
text: We use Hn for the upper half space model {(x1, · · · , xn) : x1 > 0} with the Riemannian volume dx1···dxnxn1 . Let ω = (∞, 0, . . . , 0) and o be (1, 0, . . . , 0) as in (A.2), then ξλ,ω = {(x1, x2, . . . , xn) : x1 = eλ}. The induced Riemannian metric on ξλ,ω (respectively volume dVolξλ,ω ) is dx22+···+dx 2 n e2λ (respectively dx2···dxn e(n−1)λ ). For any integral function f on Hn, using change of variable x1 = eλ∫ Hn f(x1, . . . , xn) dx1 · · · dxn xn1 = ∫ λ ∫ (x2,...,xn)∈Rn−1 f(eλ, x2, . . . , xn) dx2 · · · dxn enλ eλdλ
= ∫ λ ∫ (x2,...,xn)∈Rn−1 f(eλ, x2, . . . , xn) dx2 · · · dxn e(n−1)λ dλ
= ∫ λ ∫ ξλ,ω f(z)dVolξλ,ω (z)dλ.
The above identity is equivalent to the integral formula ∫ Hn
f(x)dVol(x) =∫ R (∫ ξλ,ω f(z)dVolξλ,ω (z) ) dλ. presented in (7), according to the Riemannian isometry in (A.2).

heading: A.5 THE HEURISTIC IS NOT A PROOF
text: The spectral theory does not directly lead to universal approximation theorems because of the following: 1, superpositions in (1, 2) and (8, 9) are different (similarly, although another kind of superposition in Hilbert’s 13th problem (Hilbert, 1935; Arnold, 2009) was a driving force for universal approximation theorems (Nielsen, 1987), the former is hardly relevant for networks (Girosi & Poggio, 1989)); 2, desired representation properties of hyperbolic eigenfunctions are unknown, partially because Hn is non-compact; 3, results in spectral theory favor Hilbert spaces, while universal approximation theorems embrace more than L2 space.

heading: A.6 OPTIMIZATION
text: The parameters update for the horocycle unit (2) involves the optimization problem on the sphere (for ω) and the hyperbolic space (for x). We use a standard algorithm of sphere optimization (Absil et al., 2008) to update ω, and in the supplement we present an optimization approach based on the geodesic polar-coordinates to update x.
In the implementation of a horocycle layer, the forward propagation is trivial, while the backpropagation involves optimization on the sphere and hyperbolic space. In the following, η is the learning rate, αt is the value of α (α may be η, s, z, ω, . . .) at the t-th step, TpX is the tangent fiber at p, ∇ is the gradient, and∇H is the hyperbolic gradient. It suffices to consider the layer s=〈z, ω〉.
Optimization on the sphere The parameter update of ω in s=〈z, ω〉 involves the optimization on the sphere. The projection of ∂Lθ∂s ∇s(ωt) = ∂Lθ ∂s zt−ωt |zt−ωt|2 ∈ TωtR n onto TωtS n−1 is given by Absil et al. (2008)[p.48]
vt = ∂Lθ ∂s zt − ωt |zt − ωt|2 − ∂Lθ ∂s ( zt − ωt |zt − ωt|2 , ωt ) ωt = ∂Lθ ∂s zt − (zt, ωt)ωt |zt − ωt|2 .
Two well-known update algorithms of wt Absil et al. (2008)[p.76] are:
ωt+1 = cos (ηt|vt|)ωt − sin (ηt|vt|)|vt|−1vt; ωt+1 = (ωt − ηtvt)/|ωt − ηtvt|.

heading: A.7 A PROOF OF APOLLONIUS THEOREM
text: Theorem 3 (Apollonius). Given distinct ω1, ω2 ∈ Sn−1 and a positive number λ, the locus {x : |x− ω1| = λ|x− ω2|} is a sphere orthogonal to Sn−1.
Proof. If λ is one then it is trivial. We assume now λ is not one. By |x− ω1| = λ|x− ω2|, we can have ∣∣∣∣x− ω1 − λω21− λ
∣∣∣∣2 = |ω1 − λω2|2|1− λ|2 − 1. The locus is a sphere with center O = ω1−λω21−λ and radius R = √ |ω1−λω2|2 |1−λ|2 − 1. The theorem of Apollonius (in all dimension) claims that this sphere is orthogonal to Sn−1. To prove this, it suffices to prove |oO|2 = 1 +R2 (recall o is the origin of Hn), which follows from∣∣∣∣ω1 − λω21− λ ∣∣∣∣2 = √ |ω1 − λω2|2 |1− λ|2 − 1 2 + 1.
A.8 INVERSION
On Rn ∪ {∞}, given the sphere {x : |x− w0| = r}, the corresponding inversion is given by
Iv(x) = w0 + r2(x− w0) |x− w0|2 .
For x ∈ Rn ∪ {∞}, Iv(x) is called the inverse of x with respect to {x : |x− w0| = r}.

heading: A.9 PROOF OF THEOREM 2
text: Theorem 2 Let K be a compact set in Hn, and 1≤p<∞. Then finite sums of the form
F (x) = N∑ i=1 αiρ(λi〈x, ωi〉H+bi), ωi∈Sn−1, αi, λi, bi∈R
are dense in Lp(K,µ), where µ is either dVol (5) or the induced Euclidean volume.
Proof. We first treat the case ρ is sigmoidal and µ = dVol . Assume that these finite sums are not dense in Lp(K, dVol). By the Hahn-Banach theorem, there exists some nonzero h∈Lq(K, dVol), where q=p/(p − 1) if p>1 and q=∞ if p=1, such that ∫ K F (x)h(x)dVol(x) = 0 for all fi-
nite sums of the form (14). As K is a compact set, by Hölder’s inequality, ∫ K |h(x)| dVol ≤
( ∫ K dVol)1/p||h||Lq(K,dVol), which leads to h∈L1(K, dVol). Extend h to be a function H
that is defined on Hn by assigning H(x)=h(x) if x∈K and H(x)=0 if x∈Hn\K. Then H∈L1(Hn, dVol)∩Lq(Hn, dVol) and∫
Hn F (x)H(x)dVol(x) = 0 (15)
for all finite sums of the form (14). For any ω∈Sn−1 and λ, b∈R, we set Fω,λ,b(x) = ρ(λ(〈x, ω〉H−b)). These functions are uniformly bounded, as |Fω,λ,b(x)|≤1. Moreover,
lim λ→∞ Fω,λ,b(x) = { 1 if 〈x, ω〉H>b, 0 if 〈x, ω〉H<b.
(16)
According to (15), for all ω, λ, b, we have ∫ Hn
Fω,λ,b(x)H(x)dVol(x) = 0. Functions {Fω,λ,b}λ∈R converge pointwise as λ→∞, and they are uniformly bounded by |H|∈L1(Hn, dVol). By the bounded convergence theorem, for all ω∈Sn−1, b∈R, we have∫
{x:〈x,ω〉H>b} H(x)dVol(x) = 0. (17)
By the integral formula (7) (with notations defined there), (6) and (17), for all b∈R,∫ ∞ 2b (∫ ξt,ω H(z)dVolξt,ω (z) ) dt = 0. (18)
Taking the derivative of ∫∞ 2b (∫ ξt,ω H(z)dVolξt,ω (z) ) dt with respect to b, we deduce from (18) that∫
ξ2b,ω H(z)dVolξ2b,ω (z) = 0 for a.e. b∈R. In other words, the integration of H on a.e. ξ ∈ Ξω is zero. This fact is valid for all ω∈Sn−1. Therefore, the integration of H on a.e. ξ ∈ Ξ is zero. By the injectivity Theorem 1, H = 0 a.e., which contradicts our assumption. Therefore, finite sums of the form (14) are dense in Lp(K, dVol). The case ρ is ReLU, ELU or Softplus and µ = dVol follows from the above case and the fact that x 7→ ρ(x+ 1)− ρ(x) is sigmoidal. The case µ is the Euclidean volume follows from previous cases and the fact that the Euclidean volume on compact K is bounded from above by λdVol for some constant λ.

heading: A.10 UNIVERSAL APPROXIMATION THEOREM FOR POISSON NEURONS.
text: In this section, ρ is a continuous sigmoidal function (Cybenko, 1989), ReLU(Nair & Hinton, 2010), ELU(Clevert et al., 2016), or Softplus(Dugas et al., 2001). We also recall the Poisson neuron:
Pρw,λ,b(x) = ρ
( λ |w|2 − |x|2
|x− w|2 + b
) , w ∈ Rn, λ, b ∈ R.
Theorem 4. Let K be a compact set in Hn, and 1≤p<∞. Then finite sums of the form
F (x) = N∑ i=1 αiP ρ ωi,λi,bi (x), ωi∈Sn−1, αi, λi, bi∈R (19)
are dense in Lp(K,µ), where µ is either dVol (5) or the induced Euclidean volume.
Proof. We first treat the case ρ is sigmoidal and µ = dVol . Assume that these finite sums are not dense in Lp(K, dVol). By the Hahn-Banach theorem, there exists some nonzero h∈Lq(K, dVol), where q=p/(p − 1) if p>1 and q=∞ if p=1, such that ∫ K F (x)h(x)dVol(x) = 0 for all fi-
nite sums of the form (19). As K is a compact set, by Hölder’s inequality, ∫ K |h(x)| dVol ≤
( ∫ K dVol)1/p||h||Lq(K,dVol), which leads to h∈L1(K, dVol). Extend h to be a function H that is defined on Hn by assigning H(x)=h(x) if x∈K and H(x)=0 if x∈Hn\K. Then H∈L1(Hn, dVol)∩Lq(Hn, dVol) and∫
Hn F (x)H(x)dVol(x) = 0 (20)
for all finite sums of the form (19). For any ω∈Sn−1, λ ∈ R, and b > 0, we set
Fω,λ,b(x) = P ρ ω,λ,−λb(x) = ρ
( λ ( 1− |x|2 |x− ω|2 − b )) .
These functions are uniformly bounded, as |Fω,λ,b(x)|≤1. Moreover,
lim λ→∞ Fω,λ,b(x) =  1 if 1−|x|2 |x−ω|2>b,
0 if 1−|x| 2 |x−ω|2<b. (21)
According to (20), for all ω, λ, b, we have ∫ Hn
Fω,λ,b(x)H(x)dVol(x) = 0. Functions {Fω,λ,b}λ∈R converge pointwise as λ→∞, and they are uniformly bounded by |H|∈L1(Hn, dVol). By the bounded convergence theorem, for all ω∈Sn−1, b∈R, we have∫
{x:〈x,ω〉H>(log b)/2} H(x)dVol(x) = ∫ { x:
1−|x|2 |x−ω|2
>b }H(x)dVol(x) = 0. (22)
By the integral formula (7) (with notations defined there), (6) and (22), for all b∈R,∫ ∞ log b (∫ ξt,ω H(z)dVolξt,ω (z) ) dt = 0. (23) Taking the derivative of ∫∞ log b (∫ ξt,ω H(z)dVolξt,ω (z) ) dt with respect to b, we deduce from (23) that∫
ξlog b,ω H(z)dVolξlog b,ω (z) = 0 for a.e. b>0. In other words, the integration of H on a.e. ξ ∈ Ξω is zero. This fact is valid for all ω∈Sn−1. Therefore, the integration of H on a.e. ξ ∈ Ξ is zero. By the injectivity Theorem 1, H = 0 a.e., which contradicts our assumption. Therefore, finite sums of the form (19) are dense in Lp(K, dVol). The case ρ is ReLU, ELU or Softplus and µ = dVol follows from the above case and the fact that x 7→ ρ(x+ 1)− ρ(x) is sigmoidal. The case µ is the Euclidean volume follows from previous cases and the fact that the Euclidean volume on compact K is bounded from above by λdVol for some constant λ.
We refere the reader to the difference of (16) and (21), (17) and (22), and (18) and (23). However, basically the proofs are the same. The points are the integral formula (7), the injectivity Theorem 1 and the fact that level sets of horocycle/Poisson neurons are horocycles. Moreover, as a corollary of Theorem 4, we have
Corollary 2. Let K be a compact set in Rn, and 1≤p<∞. Then finite sums of the form
F (x) = N∑ i=1 αiP ρ wi,λi,bi (x), wi∈Rn, αi, λi, bi∈R
are dense in Lp(K,µ), where µ is the Euclidean volume.
Proof. Because K is compact, there exists a positive number R such that K ⊂ {x ∈ Rn : |x| < R}. By the above theorem, finite sums of the form
F (x) = N∑ i=1 αiP ρ wi,λi,bi (x), wi∈Sn−1, αi, λi, bi∈R
are dense in Lp(K/R, µ). Then the corollary follows from
Pρw,λ,b(x) = P ρ w/R,λ,b(x/R).

heading: A.11 PROOF OF THE LEMMA 1
text: Recall
f1a,p(x) = 2|a|
1− |p|2 sinh−1
( 2(−p⊕ x, a)E
(1− | − p⊕ x|2)|a|
) . (24)
The proof of Lemma 1 follows from the following direct computation.
Proof. Let t ∈ (0, 1). Take pt = tω and at = −ω, then we have
−pt ⊕ x = −t(1− 2t(ω, x)E + |x|2)ω + (1− t2)x
1− 2t(ω, x)E + t2|x|2 .
Let Ft(x) = 2(−pt⊕x,at)E
(1−|−pt⊕x|2)|at| , then
Ft(x) = 2(−pt ⊕ x, at)E
(1− | − pt ⊕ x|2)|at| =
2 t(1−2t(ω,x)E+|x| 2)−(1−t2)(x,ω)E
1−2t(ω,x)E+t2|x|2
1− |−t(1−2t(ω,x)E+|x| 2)ω+(1−t2)x|2
(1−2t(ω,x)E+t2|x|2)2
= 2t(1− 2t(ω, x)E + t2|x|2)(1− 2t(ω, x)E + |x|2)− 2(1− t2)(1− 2t(ω, x)E + t2|x|2)(x, ω)E (1− 2t(ω, x)E + t2|x|2)2 − | − t(1− 2t(ω, x)E + |x|2)ω + (1− t2)x|2 = At(x)/Bt(x),
where At, Bt are defined as the corresponding numerator and denominator. We have
At(x)|t=1 = 2|x− ω|4
Bt(x)|t=1 = 0 ∂Bt(x)/∂t|t=1 = 2|x− ω|2(|x|2 − 1).
Let Gt(x) = sinh−1(Ft(x)) + log 1−t1+t , then
Gt(x) = log
( At(x)
Bt(x) +
√ 1 + A2t (x)
B2t (x)
) + log
1− t 1 + t = log ( (1− t)At (1 + t)Bt + √ (1− t)2 (1 + t)2 + (1− t)2A2t (x) (1 + t)2B2t (x) ) .
By L’Hôpital’s rule,
lim t<1,t→1 (1− t)At(x) (1 + t)Bt(x) = −At(x) + (1− t)A′t(x) Bt(x) + (1 + t)B′t(x) ∣∣∣ t=1 = |x− ω|2 2− 2|x|2 .
Therefore,
lim t<1,t→1 Gt(x) = log
( |x− ω|2
1− |x|2
) .
For t < 1, we take pt = tω, at = −ω, ct = t 2−1 4 , dt = 1 2 log 1+t 1−t , then for all x ∈ K,
lim t<1,t→1 ctf 1 at,pt(x) + dt = limt<1,t→1 −1 2 Gt(x) = 1 2 log
( 1− |x|2
|x− ω|2
) = 〈x, ω〉H .
If there exists c1, c2 such that |ctf1at, pt(x) + dt|(= |Gt(x)|/2) ≤ c2 for all t ∈ (c1, 1), x ∈ K, then by the dominated convergence theorem, there exists t such that ||ctf1at,pt(x) + dt − 〈x, ω〉H ||Lp(K,m) < , which proves the lemma. Note that
(1− t)At(x) (1 + t)Bt(x) = 2|x− ω|4(1− t) +
∑4 j=1 Uj(x, ω)(1− t)j+1
−2|x− ω|2(|x|2 − 1)(1− t)(1 + t) + ∑4 l=2 Ll(x, ω)(1− t)l(1 + t)
= 2|x− ω|4 +
∑4 j=1 Uj(x, ω)(1− t)j
2|x− ω|2(1− |x|2)(1 + t) + ∑4 l=2 Ll(x, ω)(1− t)l−1(1 + t) ,
where Uj and Ll are continuous functions defined on K × {ω}. There exist positive numbere c3, c4 and c1 ∈ (0, 1) such that for all x ∈ K and t ∈ (c1, 1),
c3 ≤ 2|x− ω|4 ≤ c4, c3 ≤ 2|x− ω|2(1− |x|2)(1 + t) ≤ c4,
c3 2 ≥ | 4∑ j=1 Uj(x, ω)(1− t)j |,
c3 2 ≥ | 4∑ l=2 Ll(x, ω)(1− t)l−1(1 + t)|.
Therefore, for x ∈ K and t ∈ (c1, 1), we have
c3 2c4 + c3 ≤ (1− t)At(x) (1 + t)Bt(x) ≤ 2c4 + c3 c3 .
This implies that for t ∈ (c1, 1), Gt|K and therefore |ctf1at,pt + dt||K are uniformly bounded, which finishes the proof of the lemma.

heading: A.12 THE FIRST MNIST CLASSIFIER IN 6.1
text: At the preprocessing stage, we compute the projection of the 28× 28 input pattern on the 40 principal components and then scale them so that the scaled 40-dimensional PCA features are within the unit ball. In our network,
1. Input layer: scaled 40-dimensional PCA features; 2. First layer: 40 inputs/1000 outputs horocycle layer (tanh activation); 3. Last layer: 1000 inputs/10 outputs affine layer; 4. Loss: cross entroy loss.
Take learning rate = 1, learning rate decay = 0.999, and batch size = 128, and run it three times. The average test error rates after 600 epochs is 1.96%.
PCA follows LeCun et al. (1998)(C.3), where 40 PCA is used for the quadratic network. Quadratic network has a similar structure to ours, because our neuron are contructed by quotient of quadratic functions followed by log.

heading: A.13 HOROCYCLE LAYER FOLLOWED BY MLR CAN APPROXIMATE THE CLASSFICATION FUNCTION
text: Suppose the MNIST classification function M is defined on ∪9j=0Kj ⊂ H40, where Ki are relatively compact and M|Kj = j. By Theorem 2, for 0≤j≤9, there exist Fj(x) = ∑Nj i=1 αj,iρ(λj,i〈x, ωj,i〉H+bj,i) such that Fj approximates IKj , where I is the indicator function. Therefore, a network with the first (horocycle) layer given by ρ(λj,i〈x, ωj,i〉H+bj,i)(0≤j≤9, 1≤i≤Nj) followed by a classical MLR with parameters given by αj,i(0≤j≤9, 1≤i≤Nj) (with arg max for prediction) approximatesM.

heading: A.14 THE SECOND MNIST CLASSIFIER IN 6.1
text: At the preprocessing stage, we do data augmentation by letting each image 1 step toward each of its 4 corners, so that our traning set has 300000 examples. In our network,
1. Input layer: (28,28, 1); 2. First block: 32-filters 3× 3 convolution, ReLU, 2× 2 max-pooling, BatchNorm; 3. Second block: 64-filters 3× 3 convolution, ReLU, BatchNorm; 4. Thrid block: 64-filters 3× 3 convolution,ReLU,2× 2 max-pooling, BatchNorm;
5. Fourth block: 128-filters 3× 3 convolution, ReLU, 2× 2 max-pooling, BatchNorm; 6. Fifth block: FC 1000, ReLU, BatchNorm;
7. Last block: 1000 input/10 output Poisson layer, sigmoid, BatchNorm;
8. Loss: cross entroy loss.
In optimization, we take Adam(Kingma & Ba, 2015). The batch size is 128 in the first 5 epochs, and 1024 in the next 15 epochs. After 5 epochs, we set ωi in the Poisson layer to be non-trainable. We train our network five times, the average test error rate after 20 epochs is 0.35%.
The in |w| 2−|x|2
|x−w|2+ is an important hyperparameter for the numerical stability. We train this MNIST model with ∈ {10−1, 10−2, 10−4, 10−6, 10−8, 10−10, 10−20}. They all show robust performance.

heading: A.15 EXPERIMENT OF POINCARE TREE CLASSIFICATION TASK
text: Given a Poincaré embedding (Nickel & Kiela, 2017) PE : {WordNet noun} → HD of the 82114 WordNet noun nodes and given a node x, the task is to classify all other nodes as being part of the subtree rooted at x (Ganea et al., 2018a). Our model is a logistic regression, where the horocycle feature p ∈ {WordNet noun} 7→ hPE(x)(PE(p)/s) (s is a hyperparameter lying in [1, 1.5]) is the only predictor, and the dependent variable is whether p is in the subtree rooted at x. Let P be the set of all nodes in the Poincare embedding, and let p range from P .
1. Input: hPE(x)(PE(p)/s) (s is a hyperparameter.) 2. Only layer: 1 input/1 output affine layer. (two parameters: one for input, one for bias.)
3. Loss: Logistic. (with respect to 1 if p in the tree rooted at x; 0 else.)
In each training, x is one of {animal, group, location, mammal, worker}, dim is one of {2,3,5,10}, and Poincaré embeddings are from the animation_train.py of Ganea et al. (2018b) 4 (with tree=wordnet_full, model=poincare, dim=dim, seed randomly ∈ {7, 8, 9}). All nodes in the subtree rooted at x are divided into training nodes (80%) and test nodes (20%). The same splitting procedure applies for the rest nodes. We choose s that has the best training F1, and then record the corresponding test F1. For each x and dim, we do the training 100 times. The average test F1 classification scores are recorded in Table 2.
The horocycle feature performs well here because it is compatible with the Poincaré embedding algorithm. Let x be a node that is not at the origin. It seems that the Poincaré embedding algorithm tends to pull all nodes that are from the subtree rooted at x towards the direction of x|x| , therefore
y → 〈 y, x|x| 〉 H is a suitable feature for this task.

heading: A.16 END-BASED CLUSTERING IN H2
text: For MNIST, at the preprocessing stage, we do data augmentation by letting each image 1 step toward each of its 4 corners, so that our traning set has 300000 examples. Our network for H2 embedding of MNIST dataset is
1. Input layer: (28,28, 1);
2. First block: 32-filters 3× 3 convolution, ReLU, 2× 2 max-pooling, BatchNorm; 3. Second block: 64-filters 3× 3 convolution, ReLU, BatchNorm; 4. Thrid block: 64-filters 3× 3 convolution,ReLU,2× 2 max-pooling, BatchNorm; 5. Fourth block: 128-filters 3× 3 convolution, ReLU, 2× 2 max-pooling, BatchNorm; 6. Fifth block: FC 1000, ReLU, BatchNorm;
7. Sixth block: FC 2, ReLU, BatchNorm, Exp;
8. Last block: 2 input/10 output horocycle layer, sigmoid;
4https://github.com/dalab/hyperbolic_cones
9. Loss: cross entroy loss,
where Exp is the exponential map ToH2(= R2)→ H2. We apply the data augmentation as in A.14. In optimization, learning rate is 0.1, learning rate decay is 0.99, batch size is 128, epochs is 50.
Our network, data augmentation and optimization for H2 embedding of Fashion-MNIST dataset is completely the same as that for MNIST.
For MNIST and Fashion-MNIST we use sphere optimization. We would like to remark that there are interesting new features in sphere optimization. Because the S1 is compact, for any continuous function f , there exists x = argmaxS1f . The derivative of f at x vanish, so the usual optimization algorithm to find the minimum will fail in the general case. In our experiments, we solve this problem by adding the following tricks:
1. Observation: if the class Cα are all close to ω ∈ S1, and the end prototype ωα for the class Cα is around −ω, then ωα is a maximum point of the loss function and therefore can not be improved through normal SGD. We solve this problem by adopting an idea(supervised variation) of k-means clustering. In each early epochs, optimization consists of two parts. In the first part, the normal SGD applies. In the second part, we move end prototypes (ωi) to the average direction of the class (using training data).
2. Observation: if the class Cα and class Cβ are all close to ω ∈ S1, and the end prototype ωα, ωβ are also both around ω, then all points in class Cα and class Cβ , end prototypes ωα, ωβ will all be pulling to ω by the SGD, and finally the network can not distinguish class Cα and class Cβ . We solve this problem by adding a loss if two prototypes are close.
With these small tricks, our 2D end-based clustering algorithm is very stable for MNIST and FashionMNIST. We run it on MNIST 10 times, and they all get a test acc around 99% within 20 epochs.
Suppose the classification task has M classes and the prototype of the i-th class is ωi. We write down the additional loss function for the second observation as follows
i = RandomChoice({1, . . . ,M}) j = RandomChoice({1, . . . ,M} \ {i}) d = (ωi, ωj)E
LObservation2 = arctanh(10× ReLU(d− 0.9− )),
where is a small constant for numerical stability.
For CIFAR-10, our network for H2 embedding of CIFAR-10 dataset is
1. Input layer: (32,32, 3); 2. First block: ResNet-32/128 output; 3. Second block: FC 2, ReLU, BatchNorm, Exp; 4. Last block: 2 input/10 output horocycle layer; 5. Loss: cross entroy loss.
In the data augmentation, we apply horizontal/vertical shifts and horizontal flip. We use Adam. The batch size is 32 in the first 100 epochs, or 1024 in the next 50 epochs. The weights of the horocycle layer are fixed at the beginning of the training and are non-trainable, which follows an idea of Mettes et al. (2019).

heading: A.17 POISSON MLR
text: For CIFAR-10, we use a ResNet-32 structure as the feature descriptor, and we apply horizontal/vertical shifts and horizontal flip. In our network,
1. Input layer: (32,32, 3); 2. First block: ResNet-32/128 output; 3. Second block: FC 128, ReLU, BatchNorm;
4. Last block: 128 input/10 output Poisson layer, BatchNorm; 5. Loss: cross entroy loss.
We use Adam. The batch size is 32 in the first 80 epochs, or 1024 in the next 20 epochs. Test acc greater than 93.5%.
For the classification task of flowers (Tensorflow), The dataset of 3670 photos of flowers contains 5 classes: daisy, dandelion, roses, sunflowers and tulips. The keras model is
1. Input layer: (180,180, 3); 2. First block: 16-filters 3× 3 convolution, ReLU, 2× 2 max-pooling; 3. Second block: 32-filters 3× 3 convolution, ReLU, 2× 2 max-pooling; 4. Thrid block: 64-filters 3× 3 convolution,ReLU,2× 2 max-pooling; 5. Fourth block: FC 128, ReLU; 6. Last block: 128 input/10 output FC layer; 7. Loss: cross entroy loss.
Our Poisson model is
1. Input layer: (180,180, 3); 2. First block: 16-filters 3× 3 convolution, ReLU, 2× 2 max-pooling; 3. Second block: 32-filters 3× 3 convolution, ReLU, 2× 2 max-pooling; 4. Thrid block: 64-filters 3× 3 convolution,ReLU,2× 2 max-pooling; 5. Fourth block: FC 128, ReLU; 6. Last block: BatchNorm, 128 input/10 output Poisson layer, sigmoid, BatchNorm; 7. Loss: cross entroy loss.
We use 2936 photos for training and the rest 734 for testing. We train two models 5 times, and Figure 6 records the average of training and test accuracies in 20 epochs. Although the test accuracy of the Poisson model is bad in the beginning, it is fairly higher (around 4%) than the test accuracy of the keras model in the end.
A.18 Dist IS A RELATIVE DISTANCE FUNCTION
As Hn is a complete metric space, the absolute distance between any pair pf points (x, ω) ∈ Hn × Sn−1 is always +∞. This absolute distance is not useful, hence we look for a relative one. Moreover, if a function D reasonably measures the relative distance of Hn from ω ∈ Sn−1, then so does D + c for any constant c. This section explains why Dist (10) is a reasonable relative distance function. For (x, ω, b) ∈ Hn × Sn−1 ×R, we recall
Dist(x, ω, b) = − log ( 1− |x|2
|x− ω|2
) + b = −2〈x, ω〉H + b.

heading: A.18.1 A NAIVE VIEWPOINT
text: For fixed ω ∈ Sn−1 and b ∈ R, Dist(·, ω, b) is defined on Hn, and its level sets are the set of horocycles that are tangential to Sn−1 at ω: Ξω = ∪λ∈R{ξλ,ω}. For λ ∈ R, using (6), we have for all x ∈ ξλ,ω ,
Dist(x, ω, b) = −λ+ b.
The above identity implies the following: if x1, x2 are on the same horocycle ξλ,ω then Dist(x1, ω, b) = Dist(x2, ω, b) = −λ + b. Moreover, as λ moves to ∞ (equivalently speaking, the horocycle ξλ,ω moves to ω), then for x ∈ ξλ,ω the function value Dist(x, ω, b) goes to −∞; as λ moves to −∞ (equivalently speaking, the ξλ,ω moves away from ω), then for x ∈ ξλ,ω the function value Dist(x, ω, b) goes to∞. The above observation heuristically explains that Dist(·, ω, b) measures the relative distance of Hn from ω.

heading: A.18.2 THE BUSEMANN FUNCTION VIEWPOINT.
text: Fix ω ∈ Sn−1, and then let c : [0,∞) → Hn be the unique geodesic ray (with unit speed) that satisfies
c[0] = (0, . . . , 0), c(∞) = ω.
For the purpose of this section, we do not need the definition of Busemann functions (we refer the interested reader to Bridson & Haefliger (2009)[II.8]). Instead, it suffices to know the following result of the theory: let dHn be the hyperbolic distance function, and for any x ∈ Hn,
lim t→∞
(dHn(x, c(t))− dHn(c(0), c(t))) = −2〈x, ω〉H . (25)
We read the above (25) in the following way: for fixed t > 0
x ∈ Hn 7→ dHn(x, c(t))− dHn(c(0), c(t))
is a function that measures the relative distance of {x, c(0)} from c(t), and therefore the left hand side of (25)
x ∈ Hn 7→ lim t→∞ (dHn(x, c(t))− dHn(c(0), c(t)))
is a function that measures the relative distance of {x, c(0)} from limt→∞ c(t) = ω, and finally the right hand side of (25) x ∈ Hn 7→ −2〈x, ω〉H is a function that measures the relative distance of {x, c(0)} from ω. Moreover, if the geodesic ray c starts from a different y ∈ Hn, and then there will be a corresponding bias term added to (25). Therefore, for any b ∈ R and ω ∈ Sn−1,
x ∈ Hn 7→ Dist(x, ω, b) = −2〈x, ω〉H + b
is a function that measures the relative distance of Hn from ω.

Review summaries:
This paper develops a MLR based on hyperbolic geometry. The idea is based on well-known concept of horocycle and horospheres which are known to be hyperbolic counterpart of line and plane in Euclidean geometry (see Coxter). Then the authors show the universal approximation which kind of follows similarly from the Euclidean counterpart. In fact we can probably conject that this universal approximation holds for any manifolds with constant sectional curvature.

This paper proposes new neural models for hyperbolic space, which unlike previous hyperbolic NN works, relies on the notion of horocycle in the Poincare disk. This novel framework has connections to spectral learnig in hyperbolic space. Representation theorems alla Cybenko for layers constructed from these neurons are presented. Finally, various experiments on clustering and classifying datasets using these neurons to generate hyperbolic embeddings are presented.