Year: 2019

Abstract:
Interest has been rising lately towards methods representing data in non-Euclidean spaces, e.g. hyperbolic or spherical, that provide specific inductive biases useful for certain real-world data properties, e.g. scale-free, hierarchical or cyclical. However, the popular graph neural networks are currently limited in modeling data only via Euclidean geometry and associated vector space operations. Here, we bridge this gap by proposing mathematically grounded generalizations of graph convolutional networks (GCN) to (products of) constant curvature spaces. We do this by i) introducing a unified formalism that can interpolate smoothly between all geometries of constant curvature, ii) leveraging gyro-barycentric coordinates that generalize the classic Euclidean concept of the center of mass. Our class of models smoothly recover their Euclidean counterparts when the curvature goes to zero from either side. Empirically, we outperform Euclidean GCNs in the tasks of node classification and distortion minimization for symbolic data exhibiting non-Euclidean behavior, according to their discrete curvature.

heading: 1 INTRODUCTION
text: Graph Convolutional Networks. The success of convolutional networks and deep learning for image data has inspired generalizations for graphs for which sharing parameters is consistent with the graph geometry. Bruna et al. (2014); Henaff et al. (2015) are the pioneers of spectral graph convolutional neural networks in the graph Fourier space using localized spectral filters on graphs. However, in order to reduce the graph-dependency on the Laplacian eigenmodes, Defferrard et al. (2016) approximate the convolutional filters using Chebyshev polynomials leveraging a result of Hammond et al. (2011). The resulting method (discussed in appendix A) is computationally efficient and superior in terms of accuracy and complexity. Further, Kipf & Welling (2017) simplify this approach by considering first-order polynomials approximations obtaining high scalability. The proposed graph convolutional networks (GCN) is interpolating node embeddings via a symmetrically normalized adjacency matrix, while this weight sharing can be understood as an efficient diffusion-like regularizer. Recent works extend GCNs to achieve state of the art results for link prediction (Zhang & Chen, 2018), graph classification (Hamilton et al., 2017; Xu et al., 2018) and node classification (Klicpera et al., 2019; Veličković et al., 2018).
Euclidean geometry in ML. In machine learning (ML), data is most often represented in a Euclidean space for various reasons. First, some data is intrinsically Euclidean, such as positions in 3D space in classical mechanics. Second, intuition is easier in such spaces, as they possess an appealing vectorial structure allowing basic arithmetic and a rich theory of linear algebra. Finally, a lot of quantities of interest such as distances and inner-products are known in closed-form formulae and can be computed very efficiently on the existing hardware. These operations are the basic building blocks for most of today’s popular machine learning models. Thus, the powerful simplicity and efficiency of Euclidean geometry has led to numerous methods achieving state-of-the-art on tasks as diverse as machine translation (Bahdanau et al., 2014; Vaswani et al., 2017), speech recognition (Graves et al., 2013), image classification (He et al., 2016) or recommender systems (He et al., 2017).
Riemannian ML. In spite of this success, certain types of data (e.g. hierarchical, scale-free or spherical data) have been shown to be better represented by non-Euclidean geometries (Defferrard et al., 2019; Bronstein et al., 2017; Nickel & Kiela, 2017; Gu et al., 2019), leading in particular to the
rich theories of manifold learning (Roweis & Saul, 2000; Tenenbaum et al., 2000) and information geometry (Amari & Nagaoka, 2007). The mathematical framework in vigor to manipulate nonEuclidean geometries is known as Riemannian geometry (Spivak, 1979). Although its theory leads to many strong and elegant results, some of its basic quantities such as the distance function d(·, ·) are in general not available in closed-form, which can be prohibitive to many computational methods.
Representational Advantages of Geometries of Constant Curvature. An interesting trade-off between general Riemannian manifolds and the Euclidean space is given by manifolds of constant sectional curvature. They define together what are called hyperbolic (negative curvature), elliptic (positive curvature) and Euclidean (zero curvature) geometries. As discussed below and in appendix B, Euclidean spaces have limitations and suffer from large distortion when embedding certain types of data such as trees, e.g. fig. 1. In these cases, the hyperbolic and spherical spaces have representational advantages providing a better inductive bias for the respective data.
The hyperbolic space can be intuitively understood as a continuous tree: the volume of a ball grows exponentially with its radius, similarly as how the number of nodes in a binary tree grows exponentially with its depth. Its tree-likeness properties have long been studied mathematically (Gromov, 1987; Hamann, 2017; Ungar, 2008) and it was proven to better embed complex networks (Krioukov et al., 2010), scale-free graphs and hierarchical data compared to the Euclidean geometry (Cho et al., 2019; Sala et al., 2018; Ganea et al., 2018b; Gu et al., 2019; Nickel & Kiela, 2018; 2017; Tifrea et al., 2019). Several important tools or methods found their hyperbolic counterparts, such as variational autoencoders (Mathieu et al., 2019; Ovinnikov, 2019), attention mechanisms (Gulcehre et al., 2018), matrix multiplications, recurrent units and multinomial logistic regression (Ganea et al., 2018a).
Similarly, spherical geometry provides benefits for modeling spherical or cyclical data (Defferrard et al., 2019; Matousek, 2013; Davidson et al., 2018; Xu & Durrett,
2018; Gu et al., 2019; Grattarola et al., 2018; Wilson et al., 2014).
Computational Efficiency of Constant Curvature Spaces (CCS). CCS are some of the few Riemannian manifolds to possess closed-form formulae for geometric quantities of interest in computational methods, i.e. distance, geodesics, exponential map, parallel transport and their gradients. We also leverage here the closed expressions for weighted centroids.
“Linear Algebra” of CCS: Gyrovector Spaces. In order to study the geometry of constant negative curvature in analogy with the Euclidean geometry, Ungar (1999; 2005; 2008; 2016) proposed the elegant non-associative algebraic formalism of gyrovector spaces. Recently, Ganea et al. (2018a) have linked this framework to the Riemannian geometry of the space, also generalizing the building blocks for non-Euclidean deep learning models operating with hyperbolic data representations.
However, it remains unclear how to extend in a principled manner the connection between Riemannian geometry and gyrovector space operations for spaces of constant positive curvature (spherical). By leveraging Euler’s formula and complex analysis, we present to our knowledge the first unified gyro framework that smoothly interpolates between geometries of constant curvatures irrespective of their signs. This is possible when working with the Poincaré ball and stereographic spherical projection models of respectively hyperbolic and spherical spaces.
How should one adapt graph neural networks to non-flat geometries of constant curvature?
In this work, we propose constant curvature GCNs to model non-Euclidean data. Node embeddings lie in spaces of constant curvature or product of those instead of a Euclidean space, thus leveraging both the representational power of these geometries and the effectiveness of GCNs.
Concurrent to our work, Chami et al. (2019); Liu et al. (2019) propose hyperbolic graph neural networks using tangent space aggregation.

heading: 2 THE GEOMETRY OF CONSTANT CURVATURE SPACES
text: Riemannian Geometry. A manifoldM of dimension d is a generalization to higher dimensions of the notion of surface, and is a space that locally looks like Rd. At each point x ∈M,M can be associated a tangent space TxM, which is a vector space of dimension d that can be understood as a first order approximation ofM around x. A riemannian metric g is given by an inner-product gx(·, ·) at each tangent space TxM, gx varying smoothly with x. A given g defines the geometry ofM, because it can be used to define the distance between x and y as the infimum of the lengths of smooth paths γ : [0, 1]→M from x to y, where the length is defined as `(γ) := ∫ 1 0 √ gγ(t)(γ̇(t), γ̇(t))dt . Under certain assumptions, a given g also defines a curvature at each point. Unifying all curvatures κ. There exist several models of respectively constant positive and negative curvatures. For positive curvature, we choose the stereographic projection of the sphere, while for negative curvature we choose the Poincaré model which is the stereographic projection of the Lorentz model. As explained below, this choice allows us to generalize the gyrovector space framework and unify spaces of both positive and negative curvature κ into a single model which we call the κ-stereographic model. The κ-stereographic model. For a curvature κ ∈ R and a dimension d ≥ 2, it is defined as stdκ = {x ∈ Rd | −κ‖x‖22 < 1} equipped with its Riemannian metric gκx = 4(1+κ||x||2)2 I =: (λ κ x)
2I. Note in particular that when κ ≥ 0, stdκ is Rd, while when κ < 0 it is the open ball of radius 1/ √ −κ. Gyrovector spaces & Riemannian geometry. As discussed in section 1, the gyrovector space formalism is used to generalize vector spaces to the Poincaré model of hyperbolic geometry (Ungar, 2005; 2008). In addition, important quantities from Riemannian geometry can be rewritten in terms of the Möbius vector addition and scalar-vector multiplication (Ganea et al., 2018a). We here extend gyrovector spaces to the κ-stereographic model, i.e. allowing positive curvature.
For κ > 0 and any point x ∈ stdκ, we will denote by x̃ the unique point of the sphere of radius κ− 1 2 in Rd+1 whose stereographic projection is x. As detailed in appendix C.2.2, it is given by
x̃ := (λκxx, κ − 12 (λκx − 1)). (1)
For x,y ∈ stdκ, we define the κ-addition, in the κ-stereographic model by:
x⊕κ y = (1− 2κxTy − κ||y||2)x + (1 + κ||x||2)y
1− 2κxTy + κ2||x||2||y||2 ∈ stdκ. (2)
The κ-addition is defined in all the cases except for spherical geometry and x = y/(κ‖y‖2) as stated by the following theorem proved in Appendix C.2.1.
Theorem 1 (Definiteness of κ-addition). We have 1 − 2κxTy + κ2||x||2||y||2 = 0 if and only if κ > 0 and x = y/(κ‖y‖2).
For s ∈ R and x ∈ stdκ (and |s tan−1κ ‖x‖| < κ 1 2π/2 if κ > 0), the κ-scaling in the κ-stereographic model is given by: s⊗κ x = tanκ (s · tan−1κ ||x||) x
||x|| ∈ stdκ, (3)
where tanκ equals κ−1/2 tan if κ > 0 and (−κ)−1/2 tanh if κ < 0. This formalism yields simple closed-forms for various quantities including the distance function inherited from the Riemannian manifold (stdκ, g κ), the exp and log maps, and geodesics, as shown by the following theorem.
Theorem 2 (Extending gyrovector spaces to positive curvature). For x,y ∈ stdκ, x 6= y, v 6= 0, (and x 6= −y/(κ‖y‖2) if κ > 0), the distance function is given bya:
dκ(x,y) = 2|κ|−1/2 tan−1κ ‖ − x⊕κ y‖, (4)
the unit-speed geodesic from x to y is unique and given by
γx→y(t) = x⊕κ (t⊗κ (−x⊕κ y)) , (5)
and finally the exponential and logarithmic maps are described as: expκx(v) = x⊕κ ( tanκ ( |κ| 12 λ κ x||v||
2
) v
||v||
) ; logκx(y) =
2|κ|− 12 λκx tan−1κ ||−x⊕κy|| −x⊕κ y || − x⊕k y||
(6) aWe write −x⊕ y for (−x)⊕ y and not −(x⊕ y).

heading: Proof sketch:
text: The case κ ≤ 0 was already taken care of by (Ganea et al., 2018a). For κ > 0, we provide a detailed proof in Appendix C.2.2. The exponential map and unit-speed geodesics are obtained using the Egregium theorem and the known formulas in the standard spherical model. The distance then follows from the formula dκ(x,y) = ‖ logκx(y)‖x which holds in any Riemannian manifold.
Around κ = 0. One notably observes that choosing κ = 0 yields all corresponding Euclidean quantities, which guarantees a continuous interpolation between κ-stereographic models of different curvatures, via Euler’s formula tan(x) = −i tanh(ix) where i := √ −1. But is this interpolation differentiable with respect to κ? It is as shown by the following theorem, proved in Appendix C.2.3.
Theorem 3 (Smoothness of stdκ w.r.t. κ around 0). Let v 6= 0 and x,y ∈ Rd, such that x 6= y (and x 6= −y/(κ‖y‖2) if κ > 0). Quantities in Eqs. (4,5,6) are well-defined for |κ| < 1/min(‖x‖2, ‖y‖2), i.e. for κ small enough. Their first order derivatives at 0− and 0+ exist and are equal. Moreover, for the distance we have:
dκ(x,y) = 2‖x− y‖ − 2κ ( ‖x− y‖3/3 + (xTy)‖x− y‖2 ) +O(κ2). (7)
Note that for xTy ≥ 0, this tells us that an infinitesimal change of curvature from zero to small negative, i.e. towards 0−, while keeping x,y fixed, has the effect of increasing their distance.
As a consequence, we have a unified formalism that interpolates smoothly between all three geometries of constant curvature.
3 κ-GCNS
We start by introducing the methods upon which we build. We present our models for spaces of constant sectional curvature, in the κ-stereographic model. However, the generalization to cartesian products of such spaces (Gu et al., 2019) follows naturally from these tools.

heading: 3.1 GRAPH CONVOLUTIONAL NETWORKS
text: The problem of node classification on a graph has long been tackled with explicit regularization using the graph Laplacian (Weston et al., 2012). Namely, for a directed graph with adjacency matrix A, by adding the following term to the loss: ∑ i,jAij‖f(xi) − f(xj)‖2 = f(X)TLf(X), where
L = D−A is the (unnormalized) graph Laplacian, Dii := ∑ k Aik defines the (diagonal) degree matrix, f contains the trainable parameters of the model and X = (xji )ij the node features of the model. Such a regularization is expected to improve generalization if connected nodes in the graph tend to share labels; node i with feature vector xi is represented as f(xi) in a Euclidean space.
With the aim to obtain more scalable models, Defferrard et al. (2016); Kipf & Welling (2017) propose to make this regularization implicit by incorporating it into what they call graph convolutional networks (GCN), which they motivate as a first order approximation of spectral graph convolutions, yielding the following scalable layer architecture (detailed in appendix A):
H(t+1) = σ ( D̃− 1 2 ÃD̃− 1 2H(t)W(t) ) , (8)
where Ã = A + I has added self-connections, D̃ii = ∑ k Ãik defines its diagonal degree matrix, σ is a non-linearity such as sigmoid, tanh or ReLU = max(0, ·), and W(t) and H(t) are the parameter and activation matrices of layer t respectively, with H(0) = X the input feature matrix.
3.2 TOOLS FOR A κ-GCN
Learning a parametrized function fθ that respects hyperbolic geometry has been studied in (Ganea et al., 2018a): neural layers and hyperbolic softmax. We generalize their definitions into the κ-stereographic model, unifying operations in positive and negative curvature. We explain how curvature introduces a fundamental difference between left and right matrix multiplications, depicting the Möbius matrix multiplication of (Ganea et al., 2018a) as a right multiplication, independent for each embedding. We then introduce a left multiplication by extension of gyromidpoints which ties the embeddings, which is essential for graph neural networks.
3.3 κ-RIGHT-MATRIX-MULTIPLICATION
Let X ∈ Rn×d denote a matrix whose n rows are d-dimensional embeddings in stdκ, and let W ∈ Rd×e denote a weight matrix. Let us first understand what a right matrix multiplication is in Euclidean space: the Euclidean right multiplication can be written row-wise as (XW)i• = Xi•W. Hence each d-dimensional Euclidean embedding is modified independently by a right matrix multiplication. A natural adaptation of this operation to the κ-stereographic model yields the following definition.
Definition 1. Given a matrix X ∈ Rn×d holding κ-stereographic embeddings in its rows and weights W ∈ Rd×e, the κ-right-matrix-multiplication is defined row-wise as
(X⊗κ W)i• = expκ0 ((log κ 0 (X)W)i•) = tanκ ( ||(XW)i•|| ||Xi•|| tan−1κ (||X•i||) ) (XW)i• ||(XW)i•||
where expκ0 and log κ 0 denote the exponential and logarithmic map in the κ-stereographic model.
This definition is in perfect agreement with the hyperbolic scalar multiplication for κ < 0, which can also be written as r ⊗κ x = expκ0 (r log κ 0 (x)). This operation is known to have desirable properties such as associativity (Ganea et al., 2018a).
3.4 κ-LEFT-MATRIX-MULTIPLICATION AS A MIDPOINT EXTENSION
For graph neural networks we also need the notion of message passing among neighboring nodes, i.e. an operation that combines / aggregates the respective embeddings together. In Euclidean space such an operation is given by the left multiplication of the embeddings matrix with the (preprocessed) adjacency Â: H(l+1) = σ(ÂZ(l)) where Z(l) = H(l)W(l). Let us consider this left multiplication. For A ∈ Rn×n, the matrix product is given row-wise by:
(AX)i• = Ai1X1• + · · ·+AinXn•
This means that the new representation of node i is obtained by calculating the linear combination of all the other node embeddings, weighted by the i-th row of A. An adaptation to the κ-stereographic model hence requires a notion of weighted linear combination. We propose such an operation in stdκ by performing a κ-scaling of a gyromidpoint − whose definition is reminded below. Indeed, in Euclidean space, the weighted linear combination αx+βy can be re-written as (α+β)mE(x,y;α, β) with Euclidean midpoint mE(x,y;α, β) := αα+βx + β α+βy. This motivates generalizing the above operation to stdκ as follows.
Definition 2. Given a matrix X ∈ Rn×d holding κ-stereographic embeddings in its rows and weights A ∈ Rn×n, the κ-left-matrix-multiplication is defined row-wise as
(A κ X)i• := ( ∑ j Aij)⊗κ mκ(X1•, · · · ,Xn•;Ai1, · · · , Ain). (9)
The κ-scaling is motivated by the fact that dκ(0, r ⊗κ x) = |r|dκ(0,x) for all r ∈ R, x ∈ stdκ. We remind that the gyromidpoint is defined when κ ≤ 0 in the κ-stereographic model as (Ungar, 2010):
mκ(x1, · · · ,xn;α1, · · · , αn) = 1
2 ⊗κ ( n∑ i=1 αiλ κ xi∑n j=1 αj(λ κ xj − 1) xi ) , (10)
with λκx = 2/(1 + κ‖x‖2). Whenever κ > 0, we have to further require the following condition:∑ j αj(λ κ xj − 1) 6= 0. (11)
For two points, one can calculate that (λκx− 1) + (λκy − 1) = 0 is equivalent to κ‖x‖‖y‖ = 1, which holds in particular whenever x = −y/(κ‖y‖2). See fig. 3 for illustrations of gyromidpoints. Our operation κ satisfies interesting properties, proved in Appendix C.2.4:
Theorem 4 (Neuter element & κ-scalar-associativity). We have In κ X = X, and for r ∈ R,
r ⊗κ (A κ X) = (rA) κ X.
The matrix A. In most graph neural networks, the matrix A is intented to be a preprocessed adjacency matrix, i.e. renormalized by the diagonal degreee matrix Dii = ∑ k Aik. This normalization is often taken either (i) to the left: D−1A, (ii) symmetric: D− 1 2AD− 1 2 or (iii) to the right: AD−1. Note that the latter case makes the matrix right-stochastic1, which is a property that is preserved by matrix product and exponentiation. For this case, we prove the following result in Appendix C.2.5:
Theorem 5 (κ-left-multiplication by right-stochastic matrices is intrinsic). If A,B are rightstochastic, φ is a isometry of stdκ and X,Y are two matrices holding κ-stereographic embeddings:
∀i, dκ ((A κ φ(X))i•, (B κ φ(Y))i•) = dκ((A κ X)i•, (B κ Y)i•). (12)
The above result means that A can easily be preprocessed as to make its κ-left-multiplication intrinsic to the metric space (stdκ, dκ). At this point, one could wonder: does there exist other ways to take weighted centroids on a Riemannian manifold? We comment on two plausible alternatives.
1M is right-stochastic if for all i, ∑
j Mij = 1.
Fréchet/Karcher means. They are obtained as arg minx ∑ i αidκ(x,xi)
2; note that although they are also intrinsic, they usually require solving an optimization problem which can be prohibitively expensive, especially when one requires gradients to flow through the solution − moreover, for the space stdκ, it is known that the minimizer is unique if and only if κ ≥ 0. Tangential aggregations. They are defined by lifting the points in a chosen tangent space via the logarithmic map, performing a linear combination and then projecting back via the exponential map, and were in particular used in the recent works of Chami et al. (2019) and Liu et al. (2019). The below theorem describes that for the κ-stereographic model, this operation is also intrinsic, i.e. commutes with isometries. We prove it in Appendix C.2.6.
Theorem 6 (Tangential aggregation is intrinsic). Define the tangential aggregation of x1, . . . ,xn ∈ stdκ w.r.t. weights {αi}1≤i≤n, at point x ∈ stdκ (for xi 6= −x/(κ‖x‖2) if κ > 0) by:
tgκx(x1, ...,xn;α1, ..., αn) := exp κ x ( n∑ i=1 αi log κ x(xi) ) . (13)
For any isometry φ of stdκ, we have
tgφ(x)({φ(xi)}; {αi}) = φ(tgx({xi}; {αi})). (14)

heading: 3.5 LOGITS
text: Finally, we need the logit and softmax layer, a neccessity for any classification task. We here use the model of (Ganea et al., 2018a), which was obtained in a principled manner for the case of negative curvature. We leave for future work the adaptation of their analysis to positive curvature and use in our experiments the straightforwardly adaptated formula to positive curvature, which we detail appendix D.
3.6 κ-GCN
We are now ready to introduce our κ-stereographic GCN (Kipf & Welling, 2017), denoted by κ-GCN2. Assume we are given a graph with node level features G = (V,A,X) where X ∈ Rn×d with each row Xi• ∈ stdκ and adjacency A ∈ Rn×n. We first perform a preprocessing step by mapping the Euclidean features to stdκ via the projection X 7→ X/(2 √ |κ|||X||max), where ||X||max denotes the maximal Euclidean norm among all stereographic embeddings in X. For l ∈ {0, . . . , L − 2}, the (l + 1)-th layer of κ-GCN is given by:
H(l+1) = σ⊗κ ( Â κ ( H(l) ⊗κ W(l) )) , (15)
where H(0) = X, σ⊗κ(x) := expκ0(σ(log κ 0(x))) is the Möbius version (Ganea et al., 2018a) of a pointwise non-linearity σ and Â = D̃− 1 2 ÃD̃− 1 2 . The final layer is a κ-logit layer (appendix D):
H(L) = softmax ( Â logitκ ( H(L−1),W(L−1) )) , (16)
where W(L−1) contains the parameters ak and pk of the κ-logits layer. A very important property of κ-GCN is that its architecture recovers the Euclidean GCN when we let curvature go to zero:
κ-GCN κ−→0−−−→ GCN.

heading: 4 EXPERIMENTS
text: We evaluate the architectures introduced in the previous sections on the tasks of node classification and minimizing embedding distortion for several synthetic as well as real datasets. We detail the training setup and model architecture choices to appendix E.
2To be pronounced “kappa” GCN; the greek letter κ being commonly used to denote sectional curvature
Minimizing Distortion Our first goal is to evaluate the graph embeddings learned by our GCN models on the representation task of fitting the graph metric in the embedding space. We desire to minimize
the average distortion, i.e. defined similarly as in (Gu et al., 2019): 1n2 ∑ i,j (( d(xi,xj) dG(i,j) )2 − 1 )2 ,
where d(xi,xj) is the distance between the embeddings of nodes i and j, while dG(i, j) is their graph distance (shortest path length).
We create three synthetic datasets that best reflect the different geometries of interest: i) “Tree‘”: a balanced tree of depth 5 and branching factor 4 consisting of 1365 nodes and 1364 edges. ii) “Torus”: We sample points (nodes) from the (planar) torus, i.e. from the unit connected square; two nodes are connected by an edge iff their toroidal distance (the warped distance) is smaller than a fixed R = 0.01; this gives 1000 nodes and 30626 edges. iii) “Spherical Graph”: we sample points (nodes) from S2, connecting nodes iff their distance is smaller than 0.2, leading to 1000 nodes and 17640 edges.
For the GCN models, we use 1-hot initial node features. We use two GCN layers with dimensions 16 and 10. The non-Euclidean models do not use additional non-linearities between layers. All euclidean parameters are updated using the ADAM optimizer with learning rate 0.01. Curvatures are learned using (stochastic) gradient descent and learning rate of 0.0001. All models are trained for 10000 epochs and we report the minimal achieved distortion. The results shown in table 1 reveal the benefit of our models. One can notice that estimated curvatures correspond to our geometric knowledge about these specific datasets.

heading: 4.1 NODE CLASSIFICATION
text: We consider the popular node classification datasets Citeseer (Sen et al., 2008), Cora-ML (McCallum et al., 2000) and Pubmed (Namata et al., 2012). Node labels correspond to the particular subfield the published document is associated with. Dataset statistics and splitting details are deferred to the appendix E due to the lack of space.
Curvature Estimations of Datasets To understand how far are the real graphs of the above datasets from the Euclidean geometry, we first estimate the graph curvature of the four studied datasets using the deviation from the Parallelogram Law (Gu et al., 2019) as detailed in appendix F. Curvature histograms are shown in fig. 4. It can be noticed that the datasets are mostly non-Euclidean, thus offering a good motivation to apply our constant-curvature GCN architectures.
Training Details We trained the Euclidean models with the hyperparameters chosen as reported in (Klicpera et al., 2019). Namely, for GCN we use one hidden layer of size 64, dropout on the embeddings and the adjacency of rate 0.5 as well as L2-regularization for the weights of the first layer with λ = 0.02. Only for Cora-ML we had to adjust the regularization factor λ to 0.002 to ensure similar scores as achieved in (Klicpera et al., 2019).
All Non-Euclidean models use biased-L2 regularization with α = 10 and λ = 2e − 2. Euclidean models used L2 regularization with the same parameter λ. We used a combination of dropout and dropconnect for the non-Euclidean models. All models have the same number of parameters. We use 2 GCN layers, hidden dimension 64. Product models split hidden dimension into [32, 32] and also input features equally. Non-Euclidean models do not use additional non-linearities. Euclidean parameters use a learning rate of 0.01 for all models using ADAM. The curvatures are learned using gradient descent with a learning rate of 0.01. We show the values of the learned curvatures in appendix E. We use early stopping: we first train for a maximum of 2000 epochs, then we check every 200 epochs for improvement in the validation cross entropy loss; if that is not observed, we stop.
Node classification results. These are shown in table 2. It can be seen that our models are competitive with the two Euclidean GCN considered (with or without non-linearities), showcasing the benefit of our proposed architecture.

heading: 5 CONCLUSION
text: In this paper, we introduced a natural extension of graph convolutional networks to the stereographic models of both positive and negative curvatures in a unified manner. We show how this choice of models permits to smoothly interpolate between positive and negative curvature, allowing the curvature of the model to be trained independent of an initial sign choice. We hope that our models will open new exciting directions into non-Euclidean graph neural networks.

heading: A GCN - A BRIEF SURVEY

heading: A.1 CONVOLUTIONAL NEURAL NETWORKS ON GRAPHS
text: One of the pioneering works on neural networks in non-Euclidean domains was done by (Defferrard et al., 2016). Their idea was to extend convolutional neural networks for graphs using tools from graph signal processing.
Given a graphG = (V,A), where A is the adjacency matrix and V is a set of nodes, we define a signal on the nodes of a graph to be a vector x ∈ Rn where xi is the value of the signal at node i. Consider the diagonalization of the symmetrized graph Laplacian L̃ = UΛUT , where Λ = diag(λ1, . . . , λn). The eigenbasis U allows to define the graph Fourier transform x̂ = UTx ∈ Rn. In order to define a convolution for graphs, we shift from the vertex domain to the Fourier domain:
x ?G y = U (( UTx ) ( UTy )) Note that x̂ = UTx and ŷ = UTy are the graph Fourier representations and we use the element-wise product since convolutions become products in the Fourier domain. The left multiplication with U maps the Fourier representation back to a vertex representation. As a consequence, a signal x filtered by gθ becomes y = Ugθ(Λ)UTx where gθ = diag(θ) with θ ∈ Rn constitutes a filter with all parameters free to vary. In order to avoid the resulting complexity O(n), (Defferrard et al., 2016) replace the non-parametric filter by a polynomial filter:
gθ(Λ) = K−1∑ k=0 θkΛ k
where θ ∈ RK resulting in a complexityO(K). Filtering a signal is unfortunately still expensive since y = Ugθ(Λ)U
Tx requires the multiplication with the Fourier basis U, thus resulting in complexity O(n2). As a consequence, (Defferrard et al., 2016) circumvent this problem by choosing the Chebyshev polynomials Tk as a polynomial basis, gθ(Λ) = ∑K k=0 θkTk(Λ̃) where Λ̃ = 2Λ λmax
− I. As a consequence, the filter operation becomes y = ∑K k=0 θkTk(L̂)x where L̂ = 2L λmax
− I. This led to a K-localized filter since it depended on the K-th power of the Laplacian. The recursive nature of these polynomials allows for an efficient filtering of complexity O(K|E|), thus leading to an computationally appealing definition of convolution for graphs. The model can also be built in an analogous way to CNNs, by stacking multiple convolutional layers, each layer followed by a non-linearity.

heading: A.2 GRAPH CONVOLUTIONAL NETWORKS
text: (Kipf & Welling, 2017) extended the work of (Defferrard et al., 2016) and inspired many follow-up architectures (Chen et al., 2018; Hamilton et al., 2017; Abu-El-Haija et al., 2018; Wu et al., 2019). The core idea of (Kipf & Welling, 2017) is to limit each filter to 1-hop neighbours by setting K = 1, leading to a convolution that is linear in the Laplacian L̂:
gθ ? x = θ0x + θ1L̂x
They further assume λmax ≈ 2, resulting in the expression
gθ ? x = θ0x− θ1D− 1 2AD− 1 2x
To additionally alleviate overfitting, (Kipf & Welling, 2017) constrain the parameters as θ0 = −θ1 = θ, leading to the convolution formula
gθ ? x = θ(I + D − 12AD− 1 2 )x
Since I+D− 1 2AD− 1 2 has its eigenvalues in the range [0, 2], they further employ a reparametrization trick to stop their model from suffering from numerical instabilities:
gθ ? x = θD̃ − 12 ÃD̃− 1 2x where Ã = A + I and D̃ii = ∑n j=1 Ãij .
Rewriting the architecture for multiple features X ∈ Rn×d1 and parameters Θ ∈ Rd1×d2 instead of x ∈ Rn and θ ∈ R, gives
Z = D̃− 1 2 ÃD̃− 1 2XΘ ∈ Rn×d2
The final model consists of multiple stacks of convolutions, interleaved by a non-linearity σ: H(k+1) = σ ( D̃− 1 2 ÃD̃− 1 2H(k)Θ(k) ) where H(0) = X and Θ ∈ Rn×dk .
The final output H(K) ∈ Rn×dK represents the embedding of each node i as hi = Hi• ∈ RdK and can be used to perform node classification:
Ŷ = softmax ( D̃− 1 2 ÃD̃− 1 2H(K)W ) ∈ Rn×L
where W ∈ RdK×L, with L denoting the number of classes.
In order to illustrate how embeddings of neighbouring nodes interact, it is easier to view the architecture on the node level. Denote by N (i) the neighbours of node i. One can write the embedding of node i at layer k + 1 as follows:
h (k+1) i = σ Θ(l) ∑ j∈Ni∪{i} h (k) j√ |N (j)||N (i)|  Notice that there is no dependence of the weight matrices Θ(l) on the node i, in fact the same parameters are shared across all nodes. In order to obtain the new embedding h(k+1)i of node i, we average over all embeddings of the neighbouring nodes. This Message Passing mechanism gives rise to a very broad class of graph neural networks (Kipf & Welling, 2017; Veličković et al., 2018; Hamilton et al., 2017; Gilmer et al., 2017; Chen et al., 2018; Klicpera et al., 2019; Abu-El-Haija et al., 2018).
To be more precise, GCN falls into the more general category of models of the form
z (k+1) i = AGGREGATE (k)({h(k)j : j ∈ N (i)};W (k))
h (k+1) i = COMBINE (k)(h (k) i , z (k+1) i ;V (k))
Models of the above form are deemed Message Passing Graph Neural Networks and many choices for AGGREGATE and COMBINE have been suggested in the literature (Kipf & Welling, 2017; Hamilton et al., 2017; Chen et al., 2018).

heading: B GRAPH EMBEDDINGS IN NON-EUCLIDEAN GEOMETRIES
text: In this section we will motivate non-Euclidean embeddings of graphs and show why the underlying geometry of the embedding space can be very beneficial for its representation. We first introduce a measure of how well a graph is represented by some embedding f : V −→ X , i 7→ f(i): Definition 3. Given an embedding f : V −→ X , i 7→ f(i) of a graph G = (V,A) in some metric space X , we call f a D-embedding for D ≥ 1 if there exists r > 0 such that
r · dG(i, j) ≤ dX (f(i), f(j)) ≤ D · r · dG(i, j)
The infimum over all such D is called the distortion of f .
The r in the definition of distortion allows for scaling of all distances. Note further that a perfect embedding is achieved when D = 1.

heading: B.1 TREES AND HYPERBOLIC SPACE
text: Trees are graphs that do not allow for a cycle, in other words there is no node i ∈ V for which there exists a path starting from i and returning back to i without passing through any node twice. The
number of nodes increases exponentially with the depth of the tree. This is a property that prohibits Euclidean space from representing a tree accurately. What intuitively happens is that "we run out of space". Consider the trees depicted in fig. 5. Here the yellow nodes represent the roots of each tree. Notice how rapidly we struggle to find appropriate places for nodes in the embedding space because their number increases just too fast.
Moreover, graph distances get extremely distorted towards the leaves of the tree. Take for instance the green and the pink node. In graph distance they are very far apart as one has to travel up all the way to the root node and back to the border. In Euclidean space however, they are very closely embedded in a L2-sense, hence introducing a big error in the embedding.
This problem can be very nicely illustrated by the following theorem: Theorem 7. Consider the tree K1,3 (also called 3-star) consisting of a root node with three children. Then every embedding {x1, . . . , x4} with xi ∈ Rk achieves at least distortion 2√3 for any k ∈ N.
Proof. We will prove this statement by using a special case of the so called Poincaré-type inequalities (Deza & Laurent, 1996):
For any b1, . . . , bk ∈ R with ∑k i=1 bi = 0 and points x1, . . . , xk ∈ Rn it holds that
k∑ i,j=1 bibj ||xi − xj ||2 ≤ 0
Consider now an embedding of the tree x1, . . . , x4 where x1 represents the root node. Choosing b1 = −3 and bi = 1 for i 6= 1 leads to the inequality ||x2 − x3||2 + ||x2 − x4||2 + ||x3 − x4||2 ≤ 3||x1 − x2||2 + 3||x1 − x3||2 + 3||x1 − x4||2
The left-hand side of this inequality in terms of the graph distance is
dG(2, 3) 2 + dG(2, 4) 2 + dG(3, 4) 2 = 22 + 22 + 22 = 12
and the right-hand side is
3 · dG(1, 2)2 + 3 · dG(1, 3)2 + 3 · dG(1, 4)2 = 3 + 3 + 3 = 9 As a result, we always have that the distortion is lower-bounded by √
12 9 = 2√ 3
Euclidean space thus already fails to capture the geometric structure of a very simple tree. This problem can be remedied by replacing the underlying Euclidean space by hyperbolic space.
Consider again the distance function in the Poincaré model, for simplicity with c = 1: dP(x, y) = cosh −1 ( 1 + 2 ||x− y||2
(1− ||x||2)(1− ||y||2)
)
Assume that the tree is embedded in the same way as in fig. 5, just restricted to lie in the disk of radius 1√
c = 1. Notice that as soon as points move closer to the boundary (||x|| −→ 1), the fraction
explodes and the resulting distance goes to infinity. As a result, the further you move points to the border, the more their distance increases, exactly as nodes on different branches are more distant to each other the further down they are in the tree. We can express this advantage in geometry in terms of distortion:
Theorem 8. There exists an embedding x1, . . . , x4 ∈ P2 for K1,3 achieving distortion 1 + for > 0 arbitrary small.
Proof. Since the Poincaré distance is invariant under Möbius translations we can again assume that x1 = 0. Let us place the other nodes on a circle of radius r. Their distance to the root is now given as
dP(xi, 0) = cosh −1 ( 1 + 2 ||xi||2
1− ||xi||2
) = cosh−1 ( 1 + 2 r2
1− r2 ) By invariance of the distance under centered rotations we can assume w.l.o.g. x2 = (r, 0). We further embed
• x3 = ( r cos( 23π), r sin( 2 3π) ) = ( − r2 , √ 3 2 r ) • x4 = ( r cos( 43π), r sin( 4 3π) ) = ( − r2 , √ 3 2 r ) .
This procedure gives:
dP(x2, x3) = cosh −1 1 + 2 || ( 3r 2 , − √ 3 2 r ) ||2
(1− r2)2
 = cosh−1(1 + 2 3r2 (1− r2)2 ) If we let the points now move to the border of the disk we observe that
cosh−1 ( 1 + 2 3r 2
(1−r2)2 ) cosh−1 ( 1 + 2 r 2
1−r2 ) r−→1−−−→ 2 But this means in turn that we can achieve distortion 1 + for > 0 arbitrary small. QED.
The tree-likeliness of hyperbolic space has been investigated on a deeper mathematical level. (Sarkar, 2011) show that a similar statement as in theorem 8 holds for all weighted or unweighted trees. The interested reader is referred to (Hamann, 2017; Sarkar, 2011) for a more in-depth treatment of the subject.
Cycles are the subclasses of graphs that are not allowed in a tree. They consist of one path that reconnects the first and the last node: (v1, . . . , vn, v1). Again there is a very simple example of a cycle, hinting at the limits Euclidean space incurs when trying to preserve the geometry of these objects (Matousek, 2013).
Theorem 9. Consider the cycle G = (V,E) of length four. Then any embedding (x1, . . . , x4) where xi ∈ Rk achieves at least distortion √ 2.
Proof. Denote by x1, x2, x3, x4 the embeddings in Euclidean space where x1, x3 and x2, x4 are the pairs without an edge. Again using the Poincaré-type inequality with b1 = b3 = 1 and b2 = b4 = −1 leads to the short diagonal theorem (Matousek, 2013):
||x1 − x3||2 + ||x2 − x4||2 ≤ ||x1 − x2||2 + ||x2 − x3||2 + ||x3 − x4||2 + ||x4 − x1||2
The left hand side of this inequality in terms of the graph distance is dG(1, 3)2 + dG(2, 4)2 = 22 + 22 = 8 and the right hand side is 12 + 12 + 12 + 12 = 4. Therefore any embedding has to shorten one diagonal by at least a factor √ 2.
It turns out that in spherical space, this problem can be solved perfectly in one dimension for any cycle.
Theorem 10. Given a cycle G = (V,E) of length n, there exists an embedding {x1, . . . , xn} achieving distortion 1.
Proof. We model the one dimension spherical space as the circle S1. Placing the points at angles 2πin and using the arclength on the circle as the distance measure leads to an embedding of distortion 1 as all pairwise distances are perfectly preserved.
Notice that we could also use the exact same embedding in the two dimensional stereographic projection model with c = 1 and we would also obtain distortion 1. The difference to the Poincaré disk is that spherical space is finite and the border does not correspond to infinitely distant points. We therefore have no since we do not have to pass to a limit.

heading: C SPHERICAL SPACE AND ITS GYROSTRUCTURE
text: Contrarily to hyperbolic geometry, spherical geometry is not only in violation with the fifth postulate of Euclid but also with the first. Notice that, shortest paths are not unique as for antipodal (oppositely situated) points, we have infinitely many geodesics connecting the two. Hence the first axiom does not hold. Notice that the third postulate holds as we stated it but it is sometimes also phrased as: "A circle of any center and radius can be constructed". Due to the finiteness of space we cannot have arbitrary large circles and hence phrased that way, the third postulate would not hold. Finally, we replace the fifth postulate by:
• Given any straight line l and a point p not on l, there exists no shortest line g passing through p but never intersecting l.
The standard model of spherical geometry suffers from the fact that its underlying space depends directly on the curvature κ through a hard constraint −κ〈x,x〉 = 1 (similarly to the Lorentz model of hyperbolic geometry). Indeed, when κ→ 0, the domain diverges to a sphere of infinite radius which is not well defined. For hyperbolic geometry, we could circumvent the problem by moving to the Poincaré model, which is the stereographic projection of the Lorentz model, relaxing the hard constraint to an inequality. A similar solution is also possible for the spherical model.

heading: C.1 STEREOGRAPHIC PROJECTION MODEL OF THE SPHERE
text: In the following we construct a model in perfect duality to the construction of the Poincaré model. Fix the south pole z = (0,− 1√
κ ) of the sphere of curvature κ > 0, i.e. of radius R := κ− 1 2 . The
stereographic projection is the map:
Φ : SnR −→ Rn,x′ 7→ x = 1
1 + √ κx′n+1 x′1:n
with the inverse given by Φ−1 : Rn −→ SnR,x 7→ x′ = ( λκxx,
1√ κ
(λκx − 1) )
where we define λκx = 2 1+κ||x||2 . Again we take the image of the sphere SnR under the extended projection Φ((0, . . . , 0,− 1κ )) = 0, leading to the stereographic model of the sphere. The metric tensor transforms as:
gκij = (λ κ x) 2δij

heading: C.2 GYROVECTOR SPACE IN THE STEREOGRAPHIC MODEL

heading: C.2.1 PROOF OF THEOREM 1
text: Using Cauchy-Schwarz’s inequality, we have A := 1− 2κxTy+κ2||x||2||y||2 ≥ 1− 2|κ|‖x‖‖y‖+ κ2||x||2||y||2 = (1 − |κ|‖x‖‖y‖)2 ≥ 0. Since equality in the Cauchy-Schwarz inequality is only reached for colinear vectors, we have that A = 0 is equivalent to κ > 0 and x = y/(κ‖y‖2).

heading: C.2.2 PROOF OF THEOREM 2
text: Let us start by proving that for x ∈ Rn and v ∈ TxRn the exponential map is given by
expκx(v) = λκx
( cosκ (λ κ x||v||)− √ κxT v||v|| sinκ (λ κ x||v||) ) x + 1√ κ sinκ (λ κ x||v||) v||v||
1 + (λκx − 1) cosκ (λκx||v||)− √ κλκxx T v ||v|| sinκ (λ κ x||v||)
(17)
Indeed, take a unit speed geodesic γx,v(t) starting from x with direction v. Notice that the unit speed geodesic on the sphere starting from x′ ∈ Sn−1 is given by Γx′,v′(t) = x′ cosκ(t) + 1√κ sinκ(t)v
′. By the Egregium theorem, we know that Φ(γx,v(t)) is again a unit speed geodesic in the sphere
where Φ−1 : x 7→ x′ = ( λκxx, 1√ κ (λκx − 1) )
. Hence Φ(γx,v(t)) is of the form of Γ for some x′ and v′. We can determine those by
x′ = Φ−1(γ(0)) = Φ−1(x) = ( λκxx,
1√ κ
(λκx − 1) )
v′ = Γ̇(0) = ∂Φ−1(y)
∂y γ(0)γ̇(0)
Notice that∇xλκx = −κ(λκx)2x and we thus get
v′ = ( −2κ(λκx)2xTvx + λκxv − √ κ(λκx) 2xTv ) We can obtain γx,v again by inverting back by calculating γx,v(t) = Φ(Γx′,v′(t)), resulting in
γx,v(t) = (λκx cosκ(t)−
√ κ(λκx) 2xTv sinκ (t))x + 1√ κ λκx sinκ (t)v
1 + (λκx − 1) cosκ (t)− √ κ(λκx) 2xTv sinκ (t)
Denoting gκx(v,v) = ||v||2λκx we have that expκx(v) = γx, 1√ gκx (v,v) v
(√ gκx(v,v) ) which concludes
the proof of the above formula of the exponential map. One then notices that it can be re-written in terms of the κ-addition. The formula for the logarithmic map is easily checked by verifying that it is indeed the inverse of the exponential map. Finally, the distance formula is obtained via the well-known identity dκ(x,y) = ‖ logκx(y)‖x where ‖v‖x = √ gκx(v,v).
Note that as expected, expκx(v)→κ→0 x + v, converging to the Euclidean exponential map.

heading: C.2.3 PROOF OF THEOREM 3
text: We first compute a Taylor development of the κ-addition w.r.t κ around zero:
x⊕κ y = (1− 2κxTy − κ||y||2)x + (1 + κ||x||2)y
1− 2κxTy + κ2||x||2||y||2 (18)
= [(1− 2κxTy − κ||y||2)x + (1 + κ||x||2)y][1 + 2κxTy +O(κ2)] (19) = (1− 2κxTy − κ||y||2)x + (1 + κ||x||2)y + 2κxTy[x + y] +O(κ2) (20) = (1− κ||y||2)x + (1 + κ||x||2)y + 2κ(xTy)y +O(κ2) (21) = x + y + κ[‖x‖2y − ‖y‖2x + 2(xTy)y] +O(κ2). (22)
We then notice that using the Taylor of ‖ · ‖2, given by ‖x + v‖2 = ‖x‖2 + 〈x,v〉+O(‖v‖22) for v→ 0, we get
‖x⊕κ y‖ = ‖x + y‖+ κ〈‖x‖2y − ‖y‖2x + 2(xTy)y,x + y〉+O(κ2) (23) = ‖x + y‖+ κ(xTy)‖x + y‖2 +O(κ2). (24)
Finally Taylor developments of tanκ(|κ| 1 2u) and |κ|− 12 tan−1κ (u) w.r.t κ around 0 for fixed u yield
For κ→ 0+, tanκ(|κ| 1 2u) = κ− 1 2 tan(κ 1 2u) (25)
= κ− 1 2 (κ 1 2u+ κ 3 2u3/3 +O(κ 52 ) (26)
= u+ κu3/3 +O(κ2). (27)
For κ→ 0−, tanκ(|κ| 1 2u) = (−κ)− 12 tanh((−κ) 12u) (28)
= (−κ)− 12 ((−κ) 12u− (−κ) 32u3/3 +O(κ 52 ) (29) = u+ κu3/3 +O(κ2). (30)
The left and right derivatives match, hence even though κ 7→ |κ| 12 is not differentiable at κ = 0, the function κ 7→ tanκ(|κ| 1 2u) is. A similar analysis yields the same conclusion for κ 7→ |κ|− 12 tan−1κ (u) yielding For κ→ 0, |κ|− 12 tan−1κ (u) = u− κu3/3 +O(κ2). (31) Since a composition of differentiable functions is differentiable, we consequently obtain that ⊗κ, expκ, logκ and dκ are differentiable functions of κ, under the assumptions on x,y,v stated in Theorem 3. Finally, the Taylor development of dκ follows by composition of Taylor developments:
dκ(x,y) = 2‖κ‖− 1 2 tan−1κ (‖(−x)⊕κ y‖)
= 2(‖x− y‖+ κ((−x)Ty)‖x− y‖2)(1− (κ/3)(‖x− y‖+O(κ))2) +O(κ2) = 2(‖x− y‖+ κ((−x)Ty)‖x− y‖2)(1− (κ/3)‖x− y‖2) +O(κ2) = 2‖x− y‖ − 2κ ( (xTy)‖x− y‖2 + ‖x− y‖3/3 ) +O(κ2).

heading: C.2.4 PROOF OF THEOREM 4
text: If A = In then for all i we have ∑ j Aij = 1, hence
(In X)i• = 1
2 ⊗κ ∑ j δijλ κ xj∑ k δik(λ κ xk − 1) xj  (32) = 1 2 ⊗κ ( λκxi (λκxi − 1) xi ) (33) = 1
2 ⊗κ (2⊗κ xi) (34)
= xi (35) = (X)i•. (36)
For associativity, we first note that the gyromidpoint is unchanged by a scalar rescaling of A. The property then follows by scalar associativity of the κ-scaling.

heading: C.2.5 PROOF OF THEOREM 5
text: It is proved in (Ungar, 2005) that the gyromidpoint commutes with isometries. The exact same proof holds for positive curvature, with the same algebraic manipulations. Moreover, when the matrix A is right-stochastic, for each row, the sum over columns gives 1, hence our operation κ reduces to a gyromidpoint. As a consequence, our κ commutes with isometries in this case. Since isometries preserve distance, we have proved the theorem.

heading: C.2.6 PROOF OF THEOREM 6
text: We begin our proof by stating the left-cancellation law:
x⊕κ (−x⊕κ y) = y (37) and the following simple identity stating that orthogonal maps commute with κ-addition
Rx⊕κ Ry = R(x⊕κ y), ∀R ∈ O(d) (38)
Next, we generalize the gyro operator from Möbius gyrovector spaces as defined in Ungar (2008):
gyr[u,v]w := −(u⊕κ v)⊕κ (u⊕κ (v ⊕κ w)) (39)
Note that this definition applies only for u,v,w ∈ stdκ for which the κ-addition is defined (see theorem 1). Following Ungar (2008), we have an alternative formulation (verifiable via computer algebra):
gyr[u,v]w = w + 2 Au +Bv
D . (40)
where the quantities A,B,D have the following closed-form expressions:
A = −κ2〈u,w〉‖v‖2 − κ〈v,w〉+ 2κ2〈u,v〉 · 〈v,w〉, (41)
B = −κ2〈v,w〉‖u‖2 + κ〈u,w〉, (42)
D = 1− 2κ〈u,v〉+ κ2‖u‖2‖v‖2. (43)
We then have the following relations: Lemma 11. For all u,v,w ∈ stdκ for which the κ-addition is defined we have the following relations: i) gyration is a linear map, ii) u ⊕κ v = gyr[u,v](v ⊕κ u), iii) −(z ⊕κ u) ⊕κ (z ⊕κ v) = gyr[z,u](−u⊕κ v), iv) ‖gyr[u,v]w‖ = ‖w‖.
Proof. The proof is similar with the one for negative curvature given in Ungar (2008). The fact that gyration is a linear map can be easily verified from its definition. For the second part, we have
−gyr[u,v](v⊕κu) = gyr[u,v](−(v⊕κu)) = −(u⊕κv)⊕κ(u⊕κ(v⊕κ(−(v⊕κu)))) = −(u⊕κv) (44) where the first equality is a trivial consequence of the fact that gyration is a linear map, while the last equality is the consequence of left-cancellation law.
The third part follows easily from the definition of the gyration and the left-cancellation law. The fourth part can be checked using the alternate form in eq. (40).
We now follow Ungar (2014) and describe all isometries of stdκ spaces:
Theorem 12. Any isometry φ of stdκ can be uniquely written as:
φ(x) = z⊕κ Rx, where z ∈ stdκ,R ∈ O(d) (45)
The proof is exactly the same as in theorems 3.19 and 3.20 of Ungar (2014), so we will skip it.
We can now prove the main theorem. Let φ(x) = z⊕κ Rx be any isometry of stdκ, where R ∈ O(d) is an orthogonal matrix. Let us denote by v := ∑n i=1 αi log κ x(xi). Then, using lemma 11 and the formula of the log map from theorem 2, one obtains the following identity:
n∑ i=1 αi log κ φ(x)(φ(xi)) = λκx λκφ(x) gyr[z,Rx]Rv (46)
and, thus, using the formula of the exp map from theorem 2 we obtain:
tgφ(x)({φ(xi)}; {αi}) = φ(x)⊕κ gyr[z,Rx]R(−x⊕κ expκx(v)) (47) Using eq. (39), we get that
gyr[z,Rx]Rw = −φ(x)⊕κ (z⊕κ (Rx⊕κ Rw)), ∀w ∈ stdκ (48)
giving the desired
tgφ(x)({φ(xi)}; {αi}) = z⊕κ R expκx(v) = φ (tgx({xi}; {αi})) (49)

heading: D LOGITS
text: The final element missing in the κ-GCN is the logit layer, a necessity for any classification task. We here use the formulation of (Ganea et al., 2018a). Denote by {1, . . . ,K} the possible labels and let ak ∈ Rd, bk ∈ R and x ∈ Rd. The output of a feed forward neural network for classification tasks is usually of the form p(y = k|x) = softmax(〈ak,x〉 − bk) In order to generalize this expression to hyperbolic space, the authors of (Ganea et al., 2018a) realized that the term in the softmax can be rewritten as
〈ak,x〉 − bk = sign(〈ak,x〉 − bk)||ak||d(x, Hak,bk)
where Ha,b = {x ∈ Rd : 〈x,a〉 − b = 0} = {x ∈ Rd : 〈−p + x,a〉 = 0} = H̃a,p with p ∈ Rd. As a first step, they define the hyperbolic hyperplane as
H̃κa,p = {x ∈ stdκ〈−p⊕κ x,a〉 = 0}
where now a ∈ Tpstdκ and p ∈ stdκ. They then proceed proving the following formula:
dκ(x, H̃a,p) = 1√ −κ
sinh−1 ( 2 √ −κ|〈−p⊕κ x,a〉|
(1 + κ|| − p⊕κ x||2)||a||
) (50)
Using this equation, they were able to obtain the following expression for the logit layer: p(y = k|x) = softmax ( ||ak||pk√ −κ sinh−1 ( 2 √ −κ〈−pk ⊕κ x,ak〉 (1 + κ|| − pk ⊕κ x||2)||ak|| )) , (51)
where ak ∈ T0stdκ ∼= Rd, x ∈ stdκ and pk ∈ stdκ. Combining all these operations leads to the definition of a hyperbolic feed forward neural network. Notice that the weight matrices W and the normal vectors ak live in Euclidean space and hence can be optimized by standard methods such as ADAM (Kingma & Ba, 2015).
For positive curvature κ > 0 we use in our experiments the following formula for the softmax layer: p(y = k|x) = softmax ( ||ak||pk√
κ sin−1
( 2 √ κ〈−pk ⊕κ x,ak〉
(1 + κ|| − pk ⊕κ x||2)||ak||
)) , (52)
which is inspired from the formula i sin(x) = sinh(ix) where i := √ −1. However, we leave for future work the rigorous proof that the distance to geodesic hyperplanes in the positive curvature setting is given by this formula.

heading: E MORE EXPERIMENTAL DETAILS
text: We here present training details for the node classification experiments.
We closely follow the training and evaluation scheme from previous work, e.g. (Klicpera et al., 2019). We split the data into training, early stopping, validation and test set. Namely we first split the dataset into a known subset of size nknown and an unknown subset consisting of the rest of the nodes. For
all the graphs we use nknown = 1500 except for MS Academics, where we use nknown = 5000. The known subset is further split into a training set consisting of 20 data points per label, an early stopping set of size 500 and a validation set of the remaining nodes. Notice that the whole structure of the graph and all the node features are used in an unsupervised fashion since the embedding of a training node might for instance depend on the embedding of a node from the validation set. But when calculating the loss, we only provide supervision with the training data. The unknown subset serves as the test data and is only used for the final evaluation of the model. Hyperparameter-tuning is performed on the validation set. We further use early stopping in all the experiments. We stop training as soon as the early stopping cross entropy loss has not decreased in the last npatience = 200 epochs or as soon as we have reached nmax = 2000 epochs. The model chosen is the one with the highest accuracy score on the early stopping set. For the final evaluation we test the model on 10 different data splits and report mean accuracy and bootstrapped confidence intervals. We use the described setup for both the Euclidean and non-Euclidean models to ensure a fair comparison.

heading: Learned curvatures .
text: Citeseer: Hyp GCN: Trained curvature, average curvature over all runs: -1.057 +-0.03 Sphr GCN: Trained curvature, average curvature over all runs: 0.951 +-0.019 Prod GCN: Trained curvatures, average curvatures: [1.331, -0.91]
Cora Hyp GCN: Trained curvature, average curvature: -1.127 +-0.011 Sphr GCN: Trained curvature, average curvature: 0.857+-0.013
Prod GCN: Trained curvature, average curvature: [-1.03, -1.01]
Pubmed: Hyp GCN: Trained curvature, average curvature: 1.123 +- 0.01 Sphr GCN: Trained curvature, average curvature: 0.896 +- 0.008 Prod GCN: Not training curvature, fixed to [-1, -1]
Ms-Academic Hyp GCN: Trained curvature, average curvature: 1.26 +- 0.09 Sphr GCN: Trained curvature, average curvature: 0.8 +- 0.07 Prod GCN: Not training curvature, fixed to [-1, -1]

heading: F GRAPH CURVATURE ESTIMATION ALGORITHM
text: We used the following procedure to estimate the curvature of a dataset developed by Gu et al. (2019):
1. Fix a node m ∈ G and sample two neighbouring nodes a, b ∈ G uniformly. Further sample an additional reference node c ∈ G uniformly (again avoiding m = c).
2. Calculate ψ(m; a, b; c) = 12dG(a,b) ( d2G(a,m) + d2G(b,c) 4 − ( d2G(a,b)+d 2 G(a,c) 2 )) 3. Reiterate the above sampling niter times and obtain an average curvature at node m. 4. Do this procedure for every node m ∈ G.

Review summaries:
This paper builds a new graph convolutional network (GCN) based on hyperbolic representations of the graph nodes: all latent representations of the nodes are points on Poincare disks. The authors adapted the Hyperbolic Neural Networks by Ganea et al. (2018) into the Kipf & Welling's (2016) version of GCN. Specifically, the authors variated the right matrix multiplication in GCN with Ganea et al.'s Mobius matrix-vector multiplication (that can be regarded as a transformation between two Poincare disks). Moreover, as a non-trivial adaptation, the author defined the left matrix multiplication in GCN (that can be regarded as a weighted linear combination of several points on the same Poincare disk) with Ungar's (2010) weighted barycenter, which is from computational geometry but not the machine learning community. The resulting method is tested on a toy problem and semi-supervised node classification of commonly used datasets, showing the possibility of improvement.

The authors propose using non-Euclidean spaces for GCNs. This is inspired by the recent work into non-Euclidean, and especially hyperbolic, embeddings. A few papers have recently tried to go past embeddings into building non-Euclidean models, requiring the lifting of standard operations in Euclidean space to non-Euclidean settings. This has been done in particular in hyperbolic space, but some datasets benefit from more complex spaces. The authors combine the mixed-curvature product formalism that uses products of Euclidean, hyperbolic, and spherical spaces for embeddings, but use these for GCN operations. 

In this paper, the authors address representation learning in non-Euclidean spaces.  The authors are motivated by constant curvature geometries, that can  provide a useful trade-off between Euclidean representations and Riemannian manifolds, i.e. arriving at more suitable representations than possible in the Euclidean space, while not sacrifising closed-form formulae for estimating distances, gradients and so on.