Title: LEARNING MIXED-CURVATURE REPRESENTATIONS IN PRODUCTS OF MODEL SPACES

Abstract: The quality of the representations achieved by embeddings is determined by how well the geometry of the embedding space matches the structure of the data. Euclidean space has been the workhorse for embeddings; recently hyperbolic and spherical spaces have gained popularity due to their ability to better embed new types of structured data-such as hierarchical data-but most data is not structured so uniformly. We address this problem by proposing learning embeddings in a product manifold combining multiple copies of these model spaces (spherical, hyperbolic, Euclidean), providing a space of heterogeneous curvature suitable for a wide variety of structures. We introduce a heuristic to estimate the sectional curvature of graph data and directly determine an appropriate signature-the number of component spaces and their dimensions-of the product manifold. Empirically, we jointly learn the curvature and the embedding in the product space via Riemannian optimization. We discuss how to define and compute intrinsic quantities such as means-a challenging notion for product manifolds-and provably learnable optimization functions. On a range of datasets and reconstruction tasks, our product space embeddings outperform single Euclidean or hyperbolic spaces used in previous works, reducing distortion by 32.55% on a Facebook social network dataset. We learn word embeddings and find that a product of hyperbolic spaces in 50 dimensions consistently improves on baseline Euclidean and hyperbolic embeddings, by 2.6 points in Spearman rank correlation on similarity tasks and 3.4 points on analogy accuracy. * . For any K > 0, the hyperboloid H d K is defined on the 1 In this paper, we use the language of graphs; note that any discrete metric space can be identified with a weighted graph, and all of our algorithms operate on weighted graphs.

1. INTRODUCTION:

With four decades of use, Euclidean space is the venerable elder of embedding spaces. Recently, non-Euclidean spaces-hyperbolic (Nickel & Kiela, 2017; Sala et al., 2018) and spherical (Wilson et al., 2014; Liu et al., 2017) -have gained attention by providing better representations for certain types of structured data. The resulting embeddings offer better reconstruction metrics: higher mean average precision (mAP) and lower distortion compared to their Euclidean counterparts. These three spaces are the model spaces of constant curvature (Lee, 1997) , and this improvement in representation fidelity arises from the correspondence between the structure of the data (hierarchical, cyclical) and the geometry of non-Euclidean space (hyperbolic: negatively curved, spherical: positively curved). The notion of curvature plays the key role.

To improve representations for a variety of types of data-beyond hierarchical or cyclical-we seek spaces with heterogeneous curvature. The motivation for such mixed spaces is intuitive: our data may have complicated, varying structure, in some regions tree-like, in others cyclical, and we seek the best of all worlds. We expect mixed spaces to match the geometry of the data and thus provide higher quality representations. However, to employ these spaces, we face several key obstacles. We must perform a challenging manifold optimization to learn both the curvature and the embedding. Afterwards, we also wish to operate on the embedded points. For example, analogy operations for word embeddings in Euclidean vector space (e.g., a -b + c) must be lifted to manifolds. We propose embedding into product spaces in which each component has constant curvature. As we show, this allows us to capture a wider range of curvatures than traditional embeddings, while retaining the ability to globally optimize and operate on the resulting embeddings. Specifically, we form a Riemannian product manifold combining hyperbolic, spherical, and Euclidean components and equip it with a decomposable Riemannian metric. While each component space in the product has constant curvature (positive for spherical, negative for hyperbolic, and zero for Euclidean), the resulting mixed space has non-constant curvature. However, selecting appropriate curvatures for the embedding space is a potential challenge. We directly learn the curvature for each component space along with the embedding (via Riemannian optimization), recovering the correct curvature, and thus the matching geometry, directly from data. We show empirically that we can indeed recover non-uniform curvatures and improve performance on reconstruction metrics.

Another technical challenge is to select the underlying number of components and dimensions of the product space; we call this the signature. This concept is vacuous in Euclidean space: the product of E r1 , . . . , E rn is identical to the single space E r1+...+rn . However, this is not the case with spherical and hyperbolic spaces. For example, the product of the spherical space S 1 (the circle) with itself is the torus S 1 × S 1 , which is topologically distinct from the sphere S 2 . We address this challenge by introducing a theory-guided heuristic estimator for the signature. We do so by matching an empirical notion of discrete curvature in our data with the theoretical distribution of the sectional curvature, a fine-grained measure of curvature on Riemannian manifolds that is amenable to analysis in products. We verify that this approach recovers the correct signature on reconstruction tasks.

Standard techniques such as PCA require centering so that the embedded directions capture components of variation. Centering in turn needs an appropriate generalization of the mean. We develop a formulation of mean for embedded points that exploits the decomposability of the distance and has theoretical guarantees. For T = {p 1 , . . . , p n } in a manifold M with dimension r, the mean is µ(T ) := arg min p i d 2 M (p, p i ). We give a global existence result: under symmetry conditions on the distribution of the points in T on the spherical components, gradient descent recovers µ(T ) with error ε in time O(nr log ε -1 ).

We demonstrate the advantages of product space embeddings through a variety of experiments; products are at least as good as single spaces, but can offer significant improvements when applied to structures not suitable for single spaces. We measure reconstruction quality (via mAP and distortion) for synthetic and real datasets over various allocations of embedding spaces. We observe a 32.55% improvement in distortion versus any single space on a Facebook social network graph. Beyond reconstruction, we apply product spaces to skip-gram word embeddings, a popular technique with numerous downstream applications, which crucially require the use of the manifold structure. We find that products of hyperbolic spaces improve performance on benchmark evaluations-suggesting that words form multiple smaller hierarchies rather than one larger one. We see an improvement of 3.4 points over baseline single spaces on the Google word analogy benchmark and of 2.6 points in Spearman rank correlation on a word similarity task using the WS-353 corpus. Our results and initial exploration suggest that mixed product spaces are a promising area for future study.

2. PRELIMINARIES & BACKGROUND:

Embeddings For metric spaces 1 U, V equipped with distances d U , d V , an embedding is a mapping f : U → V . The quality of an embedding is measured by various fidelity measures. A standard measure is average distortion D avg . The distortion of a pair of points a, b is

(|d V (f (a), f (b)) - d U (a, b)|)/d U (a, b)

, and D avg is the average over all pairs of points.

Distortion is a global metric; it considers the explicit value of all distances. At the other end of the global-local spectrum of fidelity measures is mean average precision (mAP), which applies to unweighted graphs. Let G = (V, E) be a graph and node a ∈ V have neighborhood N a = {b 1 , . . . , b deg(a) }, where deg(a) is the degree of a. In the embedding f , define R a,bi to be the smallest ball around f (a) that contains b i (that is, R a,bi is the smallest set of nearest points required to retrieve the ith neighbor of a in f ). Then, mAP(f

) = 1 |V | a∈V 1 deg(a) |Na| i=1 |N a ∩ R a,bi |/|R a,bi |.

Note that mAP does not track explicit distances; it is a ranking-based measure for local neighborhoods. Observe that mAP(f ) ≤ 1 (higher is better) while d avg ≥ 0 (lower is better).

Riemannian Manifolds:

We briefly review some notions from manifolds and Riemannian geometry. A more in-depth treatment can be found in standard texts (Lee, 2012; do Carmo, 1992) . Let M be a smooth manifold, p ∈ M be a point, and T p M be the tangent space to the point p. If M is equipped with a Riemannian metric g, then the pair (M, g) is called a Riemannian manifold. The shortest-distance paths on manifolds are called geodesics. To compute distance functions on a Riemannian manifold, the metric tensor g is integrated along the geodesic. This is a smoothly varying function (in p) g : T p M × T p M → R that induces geometric notions such as length and angle by defining an inner product on the tangent space. For example, the norm of v ∈ T p M is defined as

v g := g p (v, v) 1 2 .

In Euclidean space R d , each tangent space T p R d is canonically identified with R d , and the metric tensor g E is simply the normal inner product.

Product Manifolds Consider a sequence of smooth manifolds M 1 , M 2 , . . . , M k . The product manifold is defined as the Cartesian product M = M 1 × M 2 × . . . × M k . Notationally, we write points p ∈ M through their coordinates p = (p 1 , . . . , p k ) : p i ∈ M i , and similarly a tangent vector v ∈ T p M can be written (v 1 , . . . , v k ) : v i ∈ T pi M i . If the M i are equipped with metric tensor g i , then the product M is also Riemannian with metric tensor g(u, v) = k i=1 g i (u i , v i ). That is, the product metric decomposes into the sum of the constituent metrics.

Geodesics and Distances:

Optimization on manifolds requires a notion of taking a step. This step can be performed in the tangent space and transferred to the manifold via the exponential map Exp p : T p M → M . In a product manifold P, for tangent vectors v = (v 1 , . . . , v k ) at p = (p 1 , . . . , p k ) ∈ M , the exponential map simply decomposes, as do squared distances (Ficken, 1939; Turaga & Srivastava, 2016) :

EQUATION

In other words, the shortest path between points in the product travels along the shortest paths in each component simultaneously. Note the analogy to Euclidean products R d ≡ (R 1 ) d .

Hyperbolic and Spherical Model:

We use the hyperboloid model of hyperbolic space, with points in R d+1 . Let J ∈ R (d+1)×(d+1) be the diagonal matrix with J 00 = -1 and J ii = 1 : i > 0. For p, q ∈ R d+1 , the Minkowski inner product is p, q * := p T Jq = -p 0 q 0 + p 1 q 1 + . . . + p d q d , and the corresponding norm is p * = p, p

1 2 subset {p ∈ R d+1 : p * = -K 1/2 , p 0 > 0}.

When the subscript K is omitted, it is taken to be 1.

The hyperbolic distance on H d is d H (p, q) = acosh(-p, q * ).

Similarly, spherical space S d K is most easily defined when embedded in R d+1 . The manifold is defined on the subset {p ∈ R d+1 : p 2 = K 1/2 }, with metric g S induced by the Euclidean metric on R d+1 . The spherical distance on S d is d S (p, q) = arccos( p, q ) .

3. PRODUCT SPACES AND CONSTRUCTIONS:

We now tackle the challenges of mixed spaces. First, we introduce a product manifold embedding space P composed of multiple copies of simple model spaces, providing heterogeneous curvature. Next, in Section 3.1, given the signature of P (the number of components of each type and their dimensions), we describe how to simultaneously learn an embedding and the curvature for each component through optimization. In Section 3.2, we provide a heuristic to choose the signature by estimating a discrete notion of curvature for given data. Finally, in Section 3.3, given an embedding in P, we introduce a Karcher-style mean which can be recovered efficiently.

Let S d:

K and H d K be the spherical and hyperbolic spaces of dimension d and curvature K, -K, respectively, and E d the Euclidean space of dimension d. 2 We describe our main embedding space: for sequences of dimensions s 1 , s 2 , . . . , s m , h 1 , . . . , h n , and e, we write

P = S s1 × S s2 × • • • × S sm × H h1 × H h2 × • • • × H hn × E e

, a product manifold with m + n + 1 component spaces and total dimension i s i + j h j + e. We refer to each S si , H hi , E e as components or factors. We refer to the decomposition, e.g., (H 2 ) 2 = H 2 ×H 2 , as the signature. For convenience, let M 1 , . . . , M m+n+1 refer to the factors in the product.

Distances on P As discussed in Section 2, the product P is a Riemannian manifold defined by the structure of its components. For p, q ∈ P, we write d Mi (p, q) for the distance d Mi restricted to the appropriate components of p and q in the product. In particular, the squared distance in the product decomposes via (1). In other words, d P is simply the 2 norm of the component distances d Mi .

We note that P can also be equipped with different distances (ignoring the Riemannian structure), leading to a different embedding space. Without the underlying manifold structure, we cannot freely operate on the embedded points such as taking geodesics and means, but some simple applications only interact through distances. For such settings, we consider the 1 distance d P, 1 (p, q) = sm i=1 d Si (p, q) + hn i=1 d Hi (p, q) + d E (p, q) and the min distance d P,min (p, q) = min {d S1 (p, q), . . . , d H1 (p, q), . . . , d E (p, q)}. These distances provide simple and interpretable embedding spaces using P, enabling us to introduce combinatorial constructions that allow for embeddings without the need for optimization. We give an example below and discuss further in the Appendix. We then focus on the Riemannian distance, which allows Riemannian optimization directly on the manifold, and enables full use of the manifold structure in generic downstream applications.

Example Consider the graph G shown on the right of Figure 2 . This graph has a backbone cycle with 9 nodes, each attached to a tree; such topologies are common in networking. If a single edge (a, b) is removed from the cycle, the result is a tree embeddable arbitrarily well into hyperbolic space (Sala et al., 2018) . However, a, b (and their subtrees) would then incur an additional distance of 8 -1 = 7, being forced to go the other way around the cycle. But using the 1 distance, we can embed G tree into H 2 and G cycle into S 1 , yielding arbitrarily low distortion for G. We give the full details and another combinatorial construction for the min-distance in the Appendix.

3.1. OPTIMIZATION & COMPONENT CURVATURES:

To compute embeddings, we optimize the placement of points through an auxiliary loss function. Given graph distances {d G (X i , X j )} ij , our loss function of choice is

EQUATION

Algorithm 1 R-SGD in products 1: Input: Loss function L : P → R 2: Initialize x (0) ∈ P randomly 3: for t = 0, . . . , T -1 do 4:

h ← ∇L(x (t) )

5:

for i = 1, . . . , m do 6:

v i ← proj S x (t) i (h i ) 7:

for i = m + 1, . . . , m + n do 8:

v i ← proj H x (t) i (h i ) 9: v i ← Jv i 10: v m+n+1 ← h m+n+1

11: (2) depends on hyperbolic distance d H (for which the gradient is unstable) only through the square d 2 H , which is continuously differentiable (Sala et al., 2018) . In any Riemannian manifold, a loss function can be optimized through standard Riemannian optimization methods such as RSGD (Bonnabel, 2013) and RSVRG (Zhang et al., 2016) . We write down the full RSGD specialized to our product spaces in Algorithm 1. This proceeds by first computing the Euclidean gradient ∇L(x) with respect to the ambient space of the embedding (Step 4), and then converting it to the Riemannian gradient by applying the Riemannian correction (multiply by the inverse of the metric tensor g -1 P ). This overall strategy has been detailed in previous work in the hyperboloid model (Nickel & Kiela, 2018; Wilson & Leimeister, 2018) , and the same calculations apply to our hyperbolic components.

for i = 1, . . . , m + n + 1 do 12: x (t+1) i ← Exp x (t) i (v i ) 13: return x (T ) G G tree G cycle

Since g P is block diagonal on a product manifold, it suffices to apply the correction and perform the gradient step in each component M i independently. In the spherical and hyperboloid models, which have smaller dimension than the ambient space, this is performed by first projecting the gradient vector h onto the tangent space T x M via proj S x (h) = h -h, x x (Step 6) and proj H x (h) = h + h, x * x (Step 8). In the hyperboloid model, a final rescaling by the inverse of the metric J is needed (Step 9). This is not required in the spherical model since it inherits the same metric from the ambient Euclidean space.

Learning the Curvature There exists a spherical model for every curvature K > 0 (for example, the sphere S d K of radius K -1/2 ) and a hyperbolic model for every K < 0 (the hyperboloid H d -K ). We jointly optimize the curvature K i of every non-Euclidean factor M i along with the embeddings.

The idea is that distances on the spherical and hyperboloid models of arbitrary curvature can be emulated through distances on the standard models S, H of curvature 1. For example, given p, q on the sphere S 1/R 2 of radius R, then d(p, q) = R • d S1 (p/R, q/R) where p/R, q/R lie on the unit sphere. Therefore the radius R, which is monotone in the curvature K, can be treated as a parameter as well, so that we can optimize K and implicitly represent points lying on the manifold of curvature K, while explicitly only needing to store and optimize points in the standard model of curvature 1 via Algorithm 1. The hyperboloid model is analogous. Moreover, the loss (2) depends only on squared distances on the product manifold, which are simple functions of distances in the components through (1), so we can optimize the curvature of each factor in P. , the median am is longer in cycle-like positively curved space (Left), and shorter in tree-like negatively curved space (Right). The relative length of am can be used as a heuristic to estimate discrete curvature.

3.2. ESTIMATING THE SIGNATURE:

To choose the signature of an appropriate space P corresponding to given data, we again turn to curvature. We use the sectional curvature, a finer-grained notion defined over all two-dimensional subspaces passing through a point. Unlike coarser notions like scalar curvature, this is not constant in a product of basic spaces. Given linearly independent u, v ∈ T p M spanning a two-dimensional subspace U , the sectional curvature K p (u, v) or K p (U ) is defined as the Gaussian curvature of the surface Exp(U ) ⊆ M . Intuitively, this captures the rate that geodesics on the surface emanating from p spread apart, which relates to volume growth. In Appendix C.2, we show that the sectional curvature of P interpolates between the sectional curvatures of the factors, enabling us to better capture a wider range of structures in our embeddings:

Lemma 1. Let M = M 1 × M 2 where M i has constant curvature K i . For any u, v ∈ T p M , if K 1 , K 2 are both non-negative, the sectional curvature satisfies K(u, v) ∈ [0, max{K 1 , K 2 }]. If K 1 , K 2 are both non-positive, the sectional curvature satisfies K(u, v) ∈ [min{K 1 , K 2 }, 0]. If K i < 0 and K j > 0 for i = j, then K(u, v) ∈ [K i , K j ].

Our estimation technique employs a triangle comparison theorem following from Toponogov's theorem and the law of cosines, which characterizes sectional curvature through the behavior of small triangles (note that a triangle determines a 2-dimensional submanifold). Let abc be a geodesic triangle in manifold (or metric space) M and m be the (geodesic) midpoint of bc, and consider the quantity

EQUATION

)

This is non-negative (resp. non-positive) when the curvature is non-negative (resp. non-positive).

Note that consequently the equality case occurs exactly when the curvature is 0, and equation 3 becomes the parallelogram law of Euclidean geometry (Figure 3 ).

Analogous to sectional curvature, which is a function of a point p and two directions x, y from p, in an undirected graph G we define an analog for every node m and two neighbors b, c. Given a reference node a we set:

ξ G (m; b, c; a) = 1 2d G (a,m) ξ G (a, b, c

). This is exactly the expression from equation 3, normalized suitably so as to yield the correct scaling for trees and cycles. Our curvature estimation is then a simple average ξ

G (m; b, c) = 1 |V |-1 a =m ξ G (m; b, c; a).

Importantly, ξ G recovers the right curvature for graph atoms such as lines, cycles, and trees (Appendix C.2, Lemma 4,5), and the correct sign for other special discrete objects like polyhedra (Thurston, 1998) . The curvature is zero for lines, positive for cycles, and negative for trees.

For a generic graph G, we use this to generate a potential product manifold to embed in. An empirical sectional curvature of G is estimated via Algorithm 3, which is based off the homogeneity of product manifolds (i.e. isometries act transitively), implying that it suffices to analyze the curvature at a random point. In particular, we moment-match the distributions of sectional curvature through uniformly random 2-planes in the graph and in the manifold through Algorithms 3,2 (Appendix C.2).

3.3. MEANS IN THE PRODUCT MANIFOLD:

A critical operation on manifolds is that of taking the mean; it is necessary for many downstream applications, including, for example, analogy tasks with word embeddings, for clustering, and for centering before applying PCA. Even in simple settings like the circle S 1 , defining a mean is nontrivial. A classic approach is to take the Euclidean mean (in E 2 ) of the points and to project back onto S 1 -but this operation fails in the case where the points are uniformly spaced on S 1 . A further roadblock is the varying curvature of P. Fortunately, we can exploit the decomposability of the distance on P, reducing the challenge to breaking symmetries in the component spaces. To do so, we introduce the following Karcher-style weighted mean. Let T = {p 1 , p 2 , . . . , p n } be a set of points in P and w 1 , . . . , w n be positive weights satisfying

n i=1 w i = 1. Then the mean µ(T ) is arg min p∈P n i=1 w i d 2 P (p, p i ).

In special cases, this matches commonly used means (the centroid in the Euclidean case E d , the spherical average for S 2 in Buss & Fillmore (2001) ). We further note that when w i ≥ 0, the squared-distance components above are individually convex: this is trivial in the Euclidean term, holds in the hyperbolic case (cf. Theorem 4.1 (Bishop & O'Neill, 1969) ), and holds in the spherical case under certain restrictions, e.g., when the points in T lie entirely in one hemisphere of S r (Buss & Fillmore, 2001) . Moreover, in this case, peforming the optimization on the mean with gradient descent via the exponential map offers linear rate convergence:

Lemma 2. Let P be a product of model spaces of total dimension r, T = {p 1 , . . . , p n } points in P and w 1 , . . . , w n weights satisfying w i ≥ 0 and n i=1 w i = 1. Moreover, let the components of the points in P, p i|S j restricted to each spherical component space S j fall in one hemisphere of S j . Then, Riemannian gradient descent recovers the mean µ(T ) within distance in time O(nr log -1 ). This is a global result; with weaker assumptions, we can derive local results; for example, in the case where some of the w i are negative, which is useful for analogy operations.

In summary, we offer the following key takeaways of our development:

• Product manifolds of model spaces capture heterogeneous curvature while providing tractable optimization,

• Each component's curvature can be learned empirically through a reparametrization,

• A signature for the product can be found by matching discrete notions of curvature on graphs with sectional curvature on manifolds,

• There exists an easily-computed formulation of mean with theoretical guarantees.

4. EXPERIMENTS:

We evaluate the proposed approach, comparing the representation quality of synthetic graphs and real datasets among different embedding spaces by measuring the reconstruction fidelity (through average distortion and mAP). We expect that mixed product spaces perform better for nonhomogeneous data. We consider the curvature of graphs, reporting the curvatures learned through optimization as well as the theoretical allocation from Section 3.2. Beyond reconstruction, we evaluate the intrinsic performance of product space embeddings in a skip-gram word embedding model, by defining tasks with generic manifold operations such as means.

4.1. GRAPH RECONSTRUCTION:

Datasets We examine synthetic datasets-trees, cycles, the ring of trees shown in Figure 1 , confirming that each matches its theoretically optimal embedding space. We then compare on several real-world datasets with describable structure, including the USCA312 dataset of distances between North American cities (Burkardt); a tree-like graph of computer science Ph.D. advisor-advisee relationships (De Nooy et al., 2011) reported in previous hyperbolics work (Sala et al., 2018) ; a powergrid distribution network with backbone structure (Watts & Strogatz, 1998) ; and a dense social network from Facebook (McAuley & Leskovec, 2012) . For the former two graphs with well-defined structure, we expect optimal embeddings in spaces of positive and negative curvature, respectively. We hypothesize that the backbone network embeds well into simple products of hyperbolic and spherical spaces as in Figure 2 , and the dense graph also benefits from a mixture of spaces.

Approaches We minimize the loss (2) using Algorithm 1. We fix a total dimension d and consider the most natural ways to construct product manifolds of the given dimension, through iteratively comprising factors of dimension 2.3 For a given signature, the curvatures are initialized to the appropriate value in {-1, 0, 1} and then learned using the technique in Section 3.1. We additionally compare to the outputs of Algorithms 2,3 for heuristically selecting a combination of spaces in which to embed these datasets.

Quality We focus on the average distortion-which our loss function (2) optimizes-as our main metric for reconstruction, and additionally report the mAP metric for the unweighted graphs. As expected, for the synthetic graphs (tree, cycle, ring of trees), the matching geometries (hyperbolic, spherical, product of hyperbolic and spherical) yield the best distortion (Table 1 ). Next, we report in Table 2 the quality of embedding different graphs across a variety of allocations of spaces, fixing total dimension d = 10 following previous work (Nickel & Kiela, 2018) . We confirm that the structure of each graph informs the best allocation of spaces. In particular, the cities graph-which has intrinsic structure close to S 2 -embeds well into any space with a spherical component, and the treelike Ph.D.s graph embeds well into hyperbolic products. We emphasize that even for such datasets that theoretically match a single constant-curvature space, the products thereof perform no worse. In general, the product construction achieves high quality reconstruction: the traditional Euclidean approach is often well below several other signatures. We additionally report the learned curvatures associated with the optimal signature, finding that the resulting curvatures are non-uniform even for products of identical spaces (cf. Ph.D.s). Finally, Table 3 reports the signature estimations of Algorithms 2, 3 for the unweighted graphs. Among the signatures over two components, the estimated curvature signs agree with best distortion results from Table 2 .

4.2. WORD EMBEDDINGS:

To investigate the performance of product space embeddings in applications requiring the underlying manifold structure, we learned word embeddings and evaluated them on benchmark datasets for word similarity and analogy. In particular, we extend results on hyperbolic skip-gram embeddings from Leimeister & Wilson (2018) (LW), who found that hyperbolic embeddings perform favorably against Euclidean word vectors in low dimensions (d = 5, 20), but less so in higher dimensions (d = 50, 100). Building on these results, we hypothesize that in high dimensions, a product of multiple smaller-dimension hyperbolic spaces will substantially improve performance.

Setup We use the standard skip-gram model (Mikolov et al., 2013) and extend the loss function to a generic objective suitable for arbitrary manifolds, which is a variant of the objective proposed by LW. Concretely, given a word u and target w, with label y = 1 if w is a context word for u and y = 0 if it is a negative sample, the model is P (y|w, u) = σ (-1) 1-y (-cosh(d(α u , γ w )) + θ) .

Training followed the setup of LW, building on the fastText skip-gram implementation. Euclidean results are reported directly from fastText. Aside from choice of model, the training setup including hyperparameters (window size, negative samples, etc.) is identical to LW for all models.

Word Similarity We measure the Spearman rank correlation ρ between our scores and annotated ratings on the word similarity datasets WS-353 (Finkelstein et al., 2001) , Simlex-999 (Hill et al., Table 2 : Graph reconstruction: fidelity measures for graph embeddings using d = 10 total dimensions, with varying allocations of spaces and dimensions. Our loss function (2) targets distortion, and for each dataset the best model reflects the structure of the data. Even on near-perfectly spherical or hierarchical data, products of S (resp. H) perform no worse than the single copy. H 5 0.9 × S 5 1.6

2015) and MEN (Bruni et al., 2014) . The results are in Table 4 . Notably, we find that hyperbolic word embeddings are consistently competitive with or better than the Euclidean embeddings, and the improvement increases with more factors in the product. This suggests that word embeddings implicitly contain multiple distinct but smaller hierarchies rather than forming a single larger one.

Analogies In manifolds, there is no exact analog of the "word arithmetic" of conventional word embeddings arising from vector space structure. However, analogies can still be defined via intrinsic product manifold operations. In particular, note that the loss function depends on the embeddings solely through their pairwise distances. We thus define analogies a : b :: c : d by matching the distances

d 2 (a, b) = d 2 (c, d) and d 2 (a, c) = d 2 (b, d)

through constructing an analog of the parallelogram, by geodesically reflecting a through the geodesic midpoint (i.e. mean) m of b, c. Note that this defines both the loss function and the intrinsic tasks purely in terms of distances and manifold operations. Hence, unlike traditional word embeddings, this formulation is generic to any space.

Our evaluation, shown in Table 5 , uses the standard Google word analogy benchmark (Mikolov et al., 2013) . We observe a 22% accuracy improvement over single-space hyperbolic embeddings in 50 dimensions and similar improvements over a single hyperbolic space in 100 dimensions. As with similarity, accuracy on the analogy task consistently improves as the number of factors increases.

5. CONCLUSION:

Product spaces enable improved representations by better matching the geometry of the embedding space to the structure of the data. We introduced a tractable Riemannian product manifold class that combines Euclidean, spherical, and hyperbolic spaces. We showed how to learn embeddings and curvatures, estimate the product signature, and defined a tractable formulation of mean. We hope that our techniques encourage further research on non-Euclidean embedding spaces. The Appendix starts with a glossary of symbols and a discussion of related work. Afterwards, we provide the proof of Lemma 2. We continue with a more in-depth treatment of the curvature estimation algorithm. We then introduce two combinatorial constructions-embedding techniques that do not require optimization-that rely on the alternative product distances. We give additional details on our experimental setup. Finally, we additionally evaluate the interpretability of these embeddings (i.e., do the separate components in the embedding manifold capture intrinsic qualities of the data?) through visualizations of the synthetic example from Figure 1 .

A GLOSSARY OF SYMBOLS:

We provide a glossary of commonly-used terms in our paper. 

Symbol:



B RELATED WORK:

Hyperbolic space has recently been proposed as an alternative to Euclidean space to learn embeddings in cases where there is a (possibly latent) hierarchical structure. In fact, many types of data (from various domains) such as social networks, word frequencies, metabolic-mass relationships, and phylogenetic trees of DNA sequences exhibit a non-Euclidean latent structure, as shown in Bronstein et al. (2017) .

Initial works on hyperbolic embeddings include Nickel & Kiela (2017) and Chamberlain et al. (2017) . In Chamberlain et al. (2017) , neural graph embeddings are performed in hyperbolic space and used to classify the vertices of complex networks. A similar application is link prediction in Nickel & Kiela (2017) for the lexical database WordNet; this work also measured predicted lexical entailment on the HyperLex benchmark dataset. The follow-up work Nickel & Kiela (2018) performs optimizations in the hyperboloid (i.e. Lorentz) model instead of the Poincaré model. Tay et al. (2018) proposed a neural ranking based question answering (Q/A) system in hyperbolic space that outperformed many state-of-the-art models using fewer parameters compared to competitor learning models. Ganea et al. (2018a) proposed hyperbolic embeddings of entailment relations, described by directed acyclic graphs by applying hyperbolic cones as a heuristic and showed improvements over baselines in terms of representational capacity and generalization. Sala et al. (2018) developed a combinatorial construction for efficiently embedding trees and tree-like graphs without optimization, studied the fundamental tradeoffs of hyperbolic embeddings, and explored PCA-like algorithms in hyperbolic space.

Unlike Euclidean space, most Riemannian manifolds are not vector spaces, and thus even basic operations such as vector addition, vector translation and matrix multiplication do not have universal interpretations. In more complex geometries, closed form expressions for basic objects like distances, geodesics, and parallel transport do not exist. As a result, standard machine learning or deep learning tools, such as convolutional neural networks, long short term memory networks (LSTMs), logistic regression, support vector machines, and attention mechanisms, do not have exact correspondences in these complex geometries.

A pair of recent approaches seek to formulate standard machine learning methods in hyperbolic space. Gulcehre et al. (2018) introduces a hyperbolic version of the attention mechanism using the hyperboloid model. This work shows improvements in terms of generalization on several downstream applications including neural machine translation, learning on graphs and visual question answering tasks, while having compact neural representations. Ganea et al. (2018b) formulates basic machine learning tools in hyperbolic space including multinomial logistic regression, feed-forward and recurrent neural networks like gated recurrent units and LSTMs in order to embed sequential data and perform classification in hyperbolic space. They demonstrate empirical improvements on textual entailment and noisy-prefix recognition tasks using hyperbolic sentence embeddings. Cho et al. (2018) introduced a hyperbolic formulation for support vector machine classifiers and demonstrated performance improvements for multi-class prediction tasks on real-world complex networks as well as simulated datasets.

Zipf's law states that word-frequency distributions obey a power law, which defines a hierarchy based on semantic specificities. Concretely, semantically general words that occur in a wider range of contexts are closer to the root of the hierarchy while rarer words are further down in the hierarchy. In order to capture the latent hierarchy in the natural language, there has been several proposals for training word embeddings in hyperbolic space. Dhingra et al. (2018) trains word embeddings using the algorithm from Nickel & Kiela (2017) . They show that resulting hyperbolic word embeddings perform better on inferring lexical entailment relation than Euclidean embeddings trained with skip-gram model which is a standard method for training word embeddings, initially proposed by Mikolov et al. (2013) . Leimeister & Wilson (2018) formulated the skip-gram loss function in hyperboloid model of hyperbolic space and evaluated on the standard the intrinsic evaluation tasks for word embeddings such as similarity and analogy in hyperbolic space.

Finally, the popularity of hyperbolic embeddings has stimulated interest in descent methods suitable for hyperbolic space optimization. In addition to tools like Bonnabel (2013) and Zhang et al. (2016) , Zhang & Sra (2016) offers convergence rate analysis for a variety of algorithms and settings for Hadamard manifolds. Enokida et al. (2018) proposes an explicit update rule along geodesics in a hyperbolic space with a theoretical guarantee on convergence, and Zhang & Sra (2018) introduces an accelerated Riemannian gradient methods.

Our work also touches on previous work on maximum distance scaling (MDS) and PCA-like algorithms in hyperbolic, spherical, and more general manifolds. MDS-like algorithms in hyperbolic space are developed for visualization in Walter (2004) and Lamping & Rao (1994) . Embeddings into spherical or into hyperbolic space with a PCA-like loss function were developed in Wilson et al. (2014) . General forms of PCA include Geodesic PCA (Huckemann et al., 2010) and principal geodesics analysis (PGA) (Fletcher et al., 2004) . A very general study of PCA-like algorithms is found in Pennec (2018) .

C MANIFOLD CONCEPTS AND PROOFS:

Below, we include proofs of our results and further discuss manifold notions such as curvature.

C.1 MEANS IN PRODUCT SPACES:

We begin with Lemma 2, restated below for convenience. Lemma 2. Let P be a product of model spaces of total dimension r, T = {p 1 , . . . , p n } points in P and w 1 , . . . , w n weights satisfying w i ≥ 0 and n i=1 w i = 1. Moreover, let the components of the points in P, p i|S j restricted to each spherical component space S j fall in one hemisphere of S j . Then, Riemannian gradient descent recovers the mean µ(T ) within distance in time O(nr log -1 ).

Proof. Consider the squared distance d 2 (p, q) for p, q ∈ M for a manifold M . Fix q. We denote the Hessian in p by H p,M (q). Then, we have the following expressions for the Hessian of the squared distance of a sphere, derived in Pennec

H p,S r (q) = 2uu T + 2 θ cos θ sin θ (I r -pp T -uu t ),

where I r is the identity matrix, θ = acos(p • q) is the distance d S r (p, q) and u = (I r -pp T )q/ sin θ.

In Pennec, it is shown that the eigenvectors of H p,S r (q) are 0, 1, and θ cot(θ); thus the Hessian is bounded and if θ ∈ [0, π/2], it is also positive definite (PD).

For hyperbolic space (under the hyperboloid model), the Hessian is

H p,H r (q) = 2uu T J + 2θ coth θ(J + pp T -uu T )J.

Here θ = acosh(-p, q * ), J is the matrix associated with the Minkowski inner product, i.e., p, q * = p T Jq, and u = log p (q)/θ. The log here refers to the logarithmic map. That is, if q = exp p (v), then v = log p (q). Moreover, exact expressions for the eigenvalues of H p,H r (q) in terms of θ imply that it is always bounded and PD.

The Hessian for Euclidean space is H p,E r (q) = 2I r , which is also PD. Now we can express the Hessian of the weighted mean. We write H p,P for the Hessian of the weighted variance n i=1 w i d 2 P (p, p i ) (recall that µ(p 1 , . . . , p n ) = arg min p n i=1 w i d 2 P (p, p i )). We have, by the decomposability of the distance, that n i=1

w i d 2 P (p, p i ) = s j=1 n i=1 w i d 2 sj (p, p i ) + h j=1 n i=1 w i d 2 hj (p, p i ) + n i=1 w i d 2 E (p, p i ).

Taking the Hessian,

H p,P = s j=1 n i=1 w i H p,S s j (p i ) + h j=1 n i=1 w i H p,H h j (p i ) + 2nI e .

Now, by assumption, the spherical components for our points in each of the spheres, p i|S j , fall within one hemisphere, and we may initialize our gradient descent (that is, our p 0 ) within this hemisphere.

Then, the angle θ in each of the spherical distances is in [0, π/2], so that the corresponding Hessians are PD.

Since each term in the sum is PD and the weights satisfy w i ≥ 0, with at least one positive weight, H p,P is also PD. Moreover, these Hessians are bounded. Then we apply Theorem 4.2 (Chap. 7.4) in Udriste (1994) ), which shows linear rate convergence, as desired.

C.2 CURVATURE ESTIMATION:

We discuss the notions of curvature relevant to our product manifold in more depth. We start with a high-level overview of various definitions of curvature. Afterwards, we introduce the formal definitions for curvature and apply them to the product construction.

Definitions of Curvature There are multiple notions of curvature, with varying granularity. Some of these notions are suitable for working with manifolds abstractly (without reference to an ambient space, that is, intrinsic). Others, in particular older definitions pre-dating the development of the formal mechanisms underpinning differential geometry, require the use of the ambient space. Gauss defined the first intrinsic notion of curvature, Gaussian curvature. It is the product of the principal curvatures, which can be thought of as the smallest and largest curvature in different directions. 4 Below we consider several such notions.

Scalar curvature is a single value associated with a point p ∈ M and intuitively relates to the area of geodesic balls. Negative curvature means volumes grow faster than in Euclidean space, positive means volumes grow slower.

A more fine-grained notion of curvature is that of sectional curvature: it varies over all "sheets" passing through p. Note that curvature is inherently a notion of two-dimensional surfaces, and the sectional curvature fully captures the most general notion of curvature (the Riemannian curvature tensor). More formally, for every two dimensional subspace U of the tangent space T p M , the sectional curvature K(U ) is equal to the Gaussian curvature of the sheet Exp p (U ). Intuitively, it measures how far apart two geodesics emanating from p diverge. In positively curved spaces like the sphere, they diverge more slowly than in flat Euclidean space.

The Ricci curvature of a tangent vector v at p is the average of the sectional curvature K(U ) over all planes U containing v. Geometrically the Ricci curvature measures how much the volume of a small cone around direction v compares to the corresponding Euclidean cone. Positive curvature implies smaller volumes, and negative implies larger. Note that this is natural from the way geodesics bend in various curvatures. The scalar curvature is in fact defined as an average over the Ricci curvature, giving the intuitive relation between scalar curvature and volume. It is thus also an average over the sectional curvature.

Discrete Analogs of Curvature Discrete data such as graphs do not have manifold structure. The goal of curvature analogs such as ξ is to provide a discrete analog of curvature which satisfies similar properties to curvature; we use this to facilitate choosing an appropriate Riemannian manifold to embed discrete data into. In this work, we focus on the sectional curvature, but discrete versions of other curvatures have been proposed such as the the Forman-Ricci (Weber et al., 2017) and Ollivier-Ricci (Ollivier, 2009) curvatures.

The input to the discrete curvature estimation from Section 3.2 is analogous to other discrete curvature analogs. For example, the Ricci curvature is defined for a point p and a tangent vector u, and the coarse Ricci curvature is defined for a node p and neighbor x (Ollivier, 2011) . Similarly, the sectional curvature is defined for a point and two tangent vectors, and ξ is defined for a a node and two neighbors.

Sectional Curvature in Product Spaces Now we are ready to tackle the question of curvature in our proposed product space. Let M be our Riemannian manifold and X (M ) be the set of vector fields on M . The curvature R of M assigns a function R(X, Y ) : X (M ) → X (M ) to each pair of vector fields (X, Y ) ∈ X × X . For a vector field Z in X (M ), the function R(X, Y ) can be written

R(X, Y )Z = ∇ Y ∇ X Z -∇ X ∇ Y Z + ∇ [X,Y ] Z.

Here ∇ is the Riemannian connection for the manifold M , and [X, Y ] is the Lie bracket of the vector fields X, Y .

For convenience, we shall write the inner product R(X, Y )Z, T as (X, Y, Z, T ); this is the Riemannian curvature tensor. Then, the sectional curvature is defined as follows. Let us take V to be a two-dimensional subspace of T p M and x, y ∈ V be linearly independent (so that they span V).

Then, the sectional curvature at p for subspace V is K(x, y) := (x, y, x, y) x 2 y 2 -x, y 2 .

(4)

The model spaces S, H, E are the spaces of constant curvature, where K is constant for all points p and 2-subspaces V.

For simplicity, suppose we are working with M = M 1 × M 2 ; the approach extends easily for larger products. We write x = (x 1 , x 2 ) for x ∈ T p M . Similarly, let R 1 , R 2 be the curvatures and K 1 , K 2 be the sectional curvatures of M 1 , M 2 at p, respectively. Then the curvature tensor decomposes as

EQUATION

)

Our goal is to evaluate the sectional curvature K((x 1 , x 2 ), (y 1 , y 2 )) for the product manifold M . We show the following, re-stated for convenience:

Lemma 1. Let M = M 1 × M 2 where M i has constant curvature K i . For any u, v ∈ T p M , if K 1 , K 2 are both non-negative, the sectional curvature satisfies K(u, v) ∈ [0, max{K 1 , K 2 }]. If K 1 , K 2 are both non-positive, the sectional curvature satisfies K(u, v) ∈ [min{K 1 , K 2 }, 0]. If K i < 0 and K j > 0 for i = j, then K(u, v) ∈ [K i , K j ].

Proof. Let us start with the numerator of equation ( 4):

(x, y, x, y) = ((x 1 , x 2 ), (y 1 , y 2 ), (x 1 , x 2 ), (y 1 , y 2 )) = R((x 1 , x 2 ), (y 1 , y 2 ))((x 1 , x 2 )), (y 1 , y 2 ) = (R 1 (x 1 , y 1 )x 1 ), (R 2 (x 2 , y 2 )x 2 ) (y 1 , y 2 ) = R 1 (x 1 , y 1 )x 1 , y 1 + R 2 (x 2 , y 2 )x 2 , y 2

Here, we used equation 5 in the third line.

Note that when x 1 , y 1 are linearly independent, then R

EQUATION

2 ) by ( 4). Otherwise, this still holds since it is 0. So we can relate the above to K 1 , K 2 :

(x, y, x, y) = K 1 ( x 1 2 y 1 2 -x 1 , y 1 2 ) + K 2 ( x 2 2 y 2 2 -x 2 , y 2 2 ).

For convenience, we write α i = x i 2 y i 2 -x i , y i 2 for i = 1, 2. Then the numerator is simply K 1 α 1 + K 2 α 2 . Next, we consider the denominator of (equation 4):

x 2 y 2 -x, y 2 = (x 1 , x 2 ) 2 (y 1 , y 2 ) 2 -(x 1 , x 2 ), (y 1 , y 2 ) 2 = ( x 1 2 + x 2 2 )( y 1 2 + y 2 2 ) -( x 1 , y 1 + x 2 , y 2 ) = α 1 + α 2 + x 1 2 y 2 2 + x 2 2 y 1 2 = α 1 + α 2 + β,

where we set β = x 1 2 y 2 2 + x 2 2 y 1 2 . Thus, we have that

EQUATION

Now, note that β > 0, since we assumed that x 1 , y 1 and x 2 , y 2 are linearly independent. By Cauchy-Schwarz, α i ≥ 0. Then, if K i ≥ 0, we have that

EQUATION

)

Thus, we relate the product sectional curvature to a convex combination of the factor sectional curvatures K 1 , K 2 . We have for non-negative K 1 , K 2 (e.g., Euclidean and spherical spaces) that K((x 1 , x 2 ), (y 1 , y 2 )) ∈ [0, max{K 1 , K 2 }]. A similar result holds for the non-positive (Euclidean and hyperbolic) case. The last case (one negative, one positive space) follows along the same lines.

Distribution of K The range of curvatures from Lemma 1 can be easily extended to a more refined distributional analysis. In particular, consider sampling any point p and a random plane V ⊆ T p M . By homogeneity, we can equivalently fix p. The 2-subspaces of T p M R d forms the Grassmannian manifold Gr(2, T p M ). The uniform measure on this (i.e. invariant to multiplication by an orthogonal matrix) can be recovered from the Haar measure on the orthogonal group O(d), which Algorithm 2 Sectional curvature distribution

1: Input: Dimensions d 1 , d 2 2: a 1 ← χ 2 (d 1 -1) 3: b 1 ← χ 2 (d 1 -1) 4: t 1 ← Beta((d 1 -1)/2, (d 1 -1)/2) 5: c 1 ← a 1/2 1 b 1/2 1 (2t 1 -1) 6: a 2 ← χ 2 (d 2 -1) 7: b 2 ← χ 2 (d 2 -1) 8: t 2 ← Beta((d 2 -1)/2, (d 2 -1)/2) 9: c 2 ← a 1/2 2 b 1/2 2 (2t 2 -1) 10: α 1 ← a 1 b 1 -c 2 1 11: α 2 ← a 2 b 2 -c 2 2 12: β ← a 1 b 2 + a 2 b 1 13: return α1 α1+α2+β K 1 + α2 α1+α2+β K 2

itself can be constructed by orthonormalizing independent random normal vectors. In particular, it suffices to consider V spanned by independent Gaussians x, y ∼ N (0, I). Furthermore, we do not actually need to sample d-dimensional vectors x, y to compute the relevant curvature in equation 6. It suffices to sample the quantitities x 1 , y 1 , x 1 , x 1 , y 1 , y 1 and x 2 , y 2 , x 2 , x 2 , y 2 , y 2 directly. Note that α := x 1 , x 1 and β := y 1 , y 1 are χ 2 -distributed, while x 1 , y 1 ∼ √ αβγ, where γ is distributed as the dot product of two uniformly random unit vectors. By rotational invariance, this is the same as the first coordinate of a random unit vector, which in turn is distributed as X

2 1 /(X 2 1 + • • • + X 2 d )

for independent normal X i , and therefore (γ + 1)/2 ∼ Beta((d -1)/2, (d -1)/2).

Thus a random K(V) can be computed by sampling from well known distributions in constant time, via Algorithm 2. Furthermore, without knowing K 1 , K 2 a priori, an estimate for these curvatures can be found by matching the distribution of sectional curvature from Algorithm 2 to the empirical curvature computed from Algorithm 3. In particular, Algorithm 2 can be used to generate samples for α1 α1+α2+β and α2 α1+α2+β . The overall moments are then simple functions of K 1 , K 2 , and the sample moments of the above quantities, so that K 1 , K 2 can then be found by matching moments.

Curvature Estimation:

We prove the facts we mentioned in the main body of the paper relating to the evaluation of ξ over fundamental pieces of graphs: lines, cycles, and trees. 

; a) = 1 2 • n/2 n 2 2 + 1 - 1 2 (n/2 -1) 2 + (n/2 -1) 2 = 1.

The case when n is odd is similar, where we find that two nodes a satisfy ξ(m; b, c; a) = n/(n -1) and the rest are 0. Given a graph, the distribution of ξ(m; b, c) over random "planes" (i.e. pairs of neighbors b, c) is easily calculable. This yields a distribution that can then be averaged over m to obtain an average sectional curvature distribution. To simplify this, we find the distribution via sampling (Algorithm 3) in the calculations for Table 3 , before being fed into Algorithm 2 to estimate K i .

As a corollary of Lemma 5, note that the ξ becomes more negative for trees of higher degree, matching the intuition that higher degree trees are more appropriate for hyperbolic space embeddings (Sala et al., 2018) .

We additionally note that a line of work has studied discrete analogs of curvature for regular objects such as triangular planar tessellations (including polyhedra) (Thurston, 1998) . For example, these notions assign positive curvature to regular polyhedra (the tetrahedron, octahedron, icosahedron), zero curvature to the flat planar tessellation, and positive curvature to the order-7 triangular tiling of the hyperbolic plane. It is easily checked that Algorithm 3 assigns the right curvature sign to each of these objects.

C.3 COMBINATORIAL CONSTRUCTIONS:

The 1 and min-based distances are suitable for combinatorial constructions where we do not learn embeddings by optimizing a surrogate loss function, but rather by directly placing points in the product space, often via recursive procedures. Such constructions offer superior speed and other benefits. On the other hand,they are only available for certain classes of graphs. Additionally, since the constructions rely on the 1 and min distances, they do not take advantage of the Riemannian structure, and thus do not have the same applicability downstream.

We often exploit the combinatorial construction for trees from Sala et al. (2018) as a building block; it offers worst-case distortion5 1 + ε when embedding a tree into H r for all r ≥ 2, where we can control ε.

Hanging Tree Construction Consider the class of graphs G = (V, E) where V = B ∪ T 1 ∪ T 2 . . . ∪ T |B| , so that B is a base set of nodes and for each node a ∈ B, there is a tree T a connected to a (the hanging trees). We show how to use P with the 1 distance to reduce the cost of embedding G to that of embedding the subgraph induced by B.

We embed G into the product space P = P × H r equipped with the 1 distance. Here, P is some product manifold. We do the embedding in two steps:

1. Embed the subgraph induced by B into P by any method; let the resulting worst-case distortion be 1 + δ. Embed every node in T i into the embedded image of node i, 2. Form the tree T by connecting each of the T 1 , . . . , T |B| to a single node (equivalent to crushing all the nodes in B into a single node), and embed T into H r by using the combinatorial construction. Additionally, all of the nodes in B are embedded into the image of the single central node in H r .

We can check the distortion. For nodes x a , y b in subtrees hanging off a, b ∈ B, the distance is d G (x a , y b ) = d T (x a , y b ) + d B (x, y). Since the distortion for the two embeddings are given by 1 + δ and 1 + ε, it is easy to check that the overall distortion is at most max{1 + δ, 1 + ε}.

As a concrete example, consider the ring of trees in Figure 1 . Then, B = C r , the cycle on r nodes.

In this case, we can embed B into P = S . Let the nodes of B be indexed a 1 , . . . , a r . We embed a i into A i = (cos( 2πi d ), sin( 2πi d )). Then, for i < j, Thus indeed, the embedding has worst-case distortion 1. Thus, the overall distortion for the ring of trees is 1 + ε. Since we control ε, we can achieve arbitrarily good distortion for the ring of trees. The complexity of this algorithm is linear in the number of nodes, since embedding the trees and ring is linear time.

General Graph Construction Now we use the min distance space to construct an embedding of any graph G on r nodes with arbitrarily low distortion via the space P = H 2 × H 2 × . . . × H 2 with r -1 copies. As we shall see, this construction is ideal (arbitrarily low distortion, any graph) other than requiring O(r) spaces.

Let the nodes of G be V = {a 1 , . . . , a r }. Now, for each a i , 1 ≤ i ≤ r -1, form the minimum distance tree T i rooted at a i . Then, embed T i into the ith copy of H 2 via the combinatorial construction. Then, for any nodes a i , a j ∈ V , the distance d G (a i , a j ) is attained by d Ti (a i , a j ), or d Tj (a j , a i ) in the case i = r. Since at least one of T i or T j , say T i , is embedded in H 2 with distortion 1 + ε, if we make ε small enough, the smallest distance among the embedded copies is indeed that for T i , so our overall distortion is still 1 + ε.

D VISUALIZATIONS AND INTERPRETABILITY:

The combinatorial construction using the 1 distance (Section C.3) can embed the hanging tree graph arbitrarily well, unlike any single type of space. Unlike a single space, this also lends more interpretability to the embedding, as each component displays different qualitative aspects of the underlying graph structure. Figure 4 shows that this phenomenon does in fact happen empirically, even using the optimization approach over the 2 (Riemannian) instead of 1 distance.

E EXPERIMENTAL DETAILS:

We provide some additional details for our experimental setups.

Graph Reconstruction:

The optimization framework was implemented in PyTorch. The loss function (2) was optimized with SGD using minibatches of 65536 edges for the real-world datasets, and ran for 2000 epochs. For the Cities graph, the learning rate was chosen among {0.001, 0.003, 0.01}. For the rest of the datasets, the learning rate was chosen from a grid search among {10, 30, 100, 300, 1000} for each method.6 

Each point in the embedding is initialized randomly according to a uniform or Normal distribution in each coordinate with standard deviation 10 -3 (In the hyperboloid and spherical models, all but the first coordinate is chosen randomly, and the first coordinate is a function of the rest.)

Table 2 uses only Algorithm 1, and initializes the curvatures to -1 for hyperbolic components and 1 for spherical components. These curvatures are learned using the method described in Section 3.1, and the "Best model" row reports the final curvatures of the best signature.

Word Embeddings Following LW, the input corpus is a 2013 dump of Wikipedia that has been preprocessed by lower casing and removing punctuation, and filtered to remove articles few page views. All other hyperparameters are chosen exactly as in as LW, including their numbers for Euclidean embeddings from fastText. The datasets used for similarity (WS-353, Simlex-999, MEN) and analogy (Google) are also identical to the previous setup.

:

We write E for our Euclidean embedding space component to distinguish it from R, since our models of hyperbolic and spherical geometry also use R as an ambient space.

Note that S 1 and H 1 are metrically equivalent to R, so these are not considered.

The curvature of a curve can be found by considering the osculating circles which match it to second order.

The worst-case distortion Dwc is a commonly-considered variant of distortionDwc(f ) = max u,v∈U :u =v dV (f (u), f (v))/dU (u, v) min u,v∈U :u =v dV (f (u), f (v))/dU (u, v) .The worst-case distortion is the ratio of the worst expansion and the worst contraction of distances; note that it is scale-invariant. Here, the best worst-case distortion is Dwc(f ) = 1.

Note that the high LR stems from the particular choice of normalization for (2) in our implementation.

Figures and Tables:

Figure 1: Three component spaces: sphere S 2 , Euclidean plane E 2 , and hyperboloid H 2 . Thick lines are geodesics; these get closer in positively curved (K = +1) space S 2 , remain equidistant in flat (K = 0) space E 2 , and get farther apart in negatively curved (K = -1) space H 2 .

Figure 2: Left: Riemannian SGD decomposes per component. Subscripts i index components in the product. Right: Ring of trees graph G. Neither hyperbolic nor spherical space is suitable for G, but the product H × S captures it with low distortion. Note the decomposition into tree and cycle.

Figure3: Geodesic triangles in differently curved spaces: compared to Euclidean geometry in which it satisfies the parallelogram law (Center), the median am is longer in cycle-like positively curved space (Left), and shorter in tree-like negatively curved space (Right). The relative length of am can be used as a heuristic to estimate discrete curvature.

|V | = 1025, |E| = 1043 |V | = 4941, |E| = 6594 |V | = 4039, |E| = 88234

Suppose a lies on the same geodesic line as b, m, c; in other words, WLOG d G (a, b) ≤ d G (a, c) and suppose d G (a, c) = d G (a, b) + d G (b, m) + d G (m, c). Then ξ(m; b, c; a) = 0. Lemma 4. Consider a cycle graph C with nodes b, m, c such that (m, b) and (m, c) are neighbors. Then for all a ∈ C, ξ(m; b, c; a) is either 0 or positive. Proof. Without loss of generality, let the cycle have an even number of vertices n. Let k be the node diametrically opposite from m. Note that for any vertex a = n, a, b, m, c lie on a geodesic line, and therefore ξ(m; b, c; a) = 0. On the other hand, ξ(m; b, c

Consider a tree graph T with nodes b, m, c such that (m, b) and (m, c) are neighbors. Then for all a ∈ T , ξ(m; b, c; a) is either 0 or negative. Algorithm 3 Empirical estimation of sectional curvature distribution 1:Input: Graph G = (V, E) 2: m ← Uniform(V ) 3: b ← Uniform(N (m)) {N (v)is the neighbor set of v} 4: c ← Uniform(N (m)) 5: a ← Uniform(V ) 6: K ← ξ(m; b, c; a) 7: return K Proof. Due to the tree structure, a is either geodesic with b, m, c, or a is connected to m with a path that does not pass through b or c. In the former case, ξ(m; b, c; a) = 0. In the latter case, ξ(m; b, c; a) = 1 2d d 2 + 1 -(d + 1) 2 = -1, where d = d G (a, m).

S (A i , A j ) = acos(A i • A j )

Figure 4: Ring of trees graph embedding into (H 2 ) 1 × (S 1 ) 1 ; left: early epoch, right: completion. Only accessing graph distances, the optimization separates the different intrinsic structures of the underlying graph-the cycle and the trees-into interpretable components. Compare to Figure 2.

Matching geometries: Average distortion on canonical graphs (tree, cycle, ring of trees) with 40 nodes, comparing four spaces with total dimension 3. The best distortion is achieved by the space with matching geometry.

Table content: <table><tr><td/><td>Cycle</td><td>Tree</td><td>Ring of Trees</td></tr><tr><td/><td colspan="3">|V | = 40, |E| = 40 |V | = 40, |E| = 39 |V | = 40, |E| = 40</td></tr><tr><td>(E 3 ) 1</td><td>0.1064</td><td>0.1483</td><td>0.0997</td></tr><tr><td>(H 3 ) 1</td><td>0.1638</td><td>0.0321</td><td>0.0774</td></tr><tr><td>(S 3 ) 1</td><td>0.0007</td><td>0.1605</td><td>0.1106</td></tr><tr><td>(H 2 ) 1 × (S 1 ) 1</td><td>0.1108</td><td>0.0538</td><td>0.0616</td></tr><tr><td colspan="4">doubling the number of factors. These models include the products consisting of only a constant-</td></tr><tr><td colspan="3">curvature base space, ranging to various combinations of S d/2 2 , H</td><td>d/2 2</td></tr></table>

Heuristic allocation: estimated signatures for embedding unweighted graphs from Table2into two factors, using Algorithms 2,3 to match the empirical distribution of graph curvature. The resulting curvature signs agree with results from Table2for choosing among two-component spaces.

Table content: <table><tr><td>CS PhDs</td><td>Power</td><td>Facebook</td></tr><tr><td>Estimated Signature H 5 1.3 × H 5 0.2</td><td>H 5 1.8 × S 5 1.7</td><td/></tr></table>

Spearman rank correlation on similarity datasets. Top: Previous results from embeddings into spaces of fixed curvature. Bottom: Embeddings into products of H with fixed total dimension.

Table content: <table><tr><td/><td/><td>Dim 50</td><td/><td>Dim 100</td></tr><tr><td/><td colspan="2">WS-353 Simlex MEN</td><td colspan="2">WS-353 Simlex</td><td>MEN</td></tr><tr><td>Euclidean</td><td>0.6628</td><td colspan="2">0.2738 0.7217 0.6986</td><td>0.2923</td><td>0.7473</td></tr><tr><td>Hyperbolic</td><td>0.6787</td><td colspan="2">0.2784 0.7117 0.6846</td><td>0.2832</td><td>0.7217</td></tr><tr><td colspan="2">2 Hyperbolics 0.6955</td><td colspan="2">0.2870 0.7246 0.7297</td><td>0.3168</td><td>0.7450</td></tr><tr><td colspan="2">5 Hyperbolics 0.7048</td><td colspan="2">0.2837 0.7270 0.7379</td><td>0.3212</td><td>0.7530</td></tr></table>

Accuracy on the Google word analogy dataset. Taking products of smaller hyperbolic spaces significantly improves performance. Unlike conventional embeddings, the operations in hyperbolic and product spaces are defined solely through distances and manifold operations.

Table content: <table><tr><td colspan="2">Total Dim d / Model R d</td><td>(H d ) 1</td><td colspan="2">(H d/2 ) 2 (H d/5 ) 5 (H 2 ) d/2</td></tr><tr><td>50</td><td colspan="3">0.3866 0.3424 0.3928</td><td>0.4181</td><td>0.4209</td></tr><tr><td>100</td><td colspan="3">0.5513 0.3738 0.4310</td><td>0.4731</td><td>0.5216</td></tr></table>

Glossary of variables and symbols used in this paper.

Table content: <table><tr><td/><td>Used for</td></tr><tr><td colspan="2">mAP(f ) the mean average precision fidelity measure of the embedding f</td></tr><tr><td>D(f )</td><td>the distortion fidelity measure of the embedding f</td></tr><tr><td>D wc (f )</td><td>the worst-case distortion fidelity measure of the embedding f</td></tr><tr><td>G</td><td>a graph, typically with node set V and edge set E</td></tr><tr><td>T</td><td>a tree</td></tr><tr><td>a, b, c</td><td>nodes in a graph or tree</td></tr><tr><td>f</td><td>an embedding</td></tr><tr><td>N a</td><td>neighborhood around node a in a graph</td></tr><tr><td>R a,b</td><td/></tr><tr><td>E d</td><td>d-dimensional Euclidean space</td></tr><tr><td>S d</td><td>d-dimensional spherical space</td></tr><tr><td>H d</td><td>d-dimensional hyperbolic space</td></tr><tr><td>P</td><td>product manifold consisting of spherical, Euclidean, hyperbolic factors</td></tr><tr><td colspan="2">Exp x (v) the exponential map for tangent vector v at point x</td></tr><tr><td>R</td><td>the Riemannian curvature tensor</td></tr><tr><td>K(x, y)</td><td>the sectional curvature for a subspace spanned by linearly independent x, y ∈ T p M</td></tr><tr><td>d E</td><td>metric distance between two points in Euclidean space</td></tr><tr><td>d S</td><td>metric distance between two points in spherical space</td></tr><tr><td>d H</td><td>metric distance between two points in hyperbolic space</td></tr><tr><td>d U</td><td>metric distance between two points in metric space U</td></tr><tr><td>d</td><td/></tr></table>

Peer review (iclr2019):
 Page 2: What are p_i, i=1,2,...,n, their set T and .mathcal{P}? What is | | used to compute distortion between a and b? Please fix the definition of the Riemannian manifold, such that M is not just any manifold, but should be a smooth manifold or a particular differentiable manifold. Please update your definition more precisely, by checking page 328 in J.M. Lee, Introduction to Smooth Manifolds, 2012, or Page 38 in do Cormo, Riemannian Geometry, 1992. Please define .mathcal{P} in equation (1). Define K used in the definition of the hyperboloid more precisely. Please provide proofs of these statements for product of manifolds with nonnegative and nonpositive curvatures: “In particular, the squared distance in the product decomposes via (1). In other words, dP is simply the l2 norm of the component distances dMi.” Please explain what you mean by “without the need for optimization” in “These distances provide simple and interpretable embedding spaces using P, enabling us to introduce combinatorial constructions that allow for embeddings without the need for optimization.” In addition, how can you compute geodesic etc. if you use l1 distance for the embedded space? By equation (2), the paper focuses on embedding graphs, which is indeed the main goal of the paper. Therefore, first, the novelty and claims of the paper should be revised for graph embedding. Second, three particular spaces are considered in this work, which are the sphere, hyperbolic manifold, and Euclidean space. Therefore, you cannot simply state your novelty for a general class of product spaces. Thus, the title, novelty, claims and other parts of the paper should be revised and updated according to the particular input and output spaces of embeddings considered in the paper. Please explain how you compute the metric tensor g_P and apply the Riemannian correction (multiply by the inverse of the metric tensor g_P) to determine the Riemannian gradient in the Algorithm 1, more precisely. Step (9) of the Algorithm 1 is either wrong, or you compute v_i without projecting the Riemannian gradient. Please check your theoretical/experimental results and code according to this step. What is h_i used in the Algorithm 1? Can we suppose that it is the ith component of h? In step (6) and step (8), do you project individual components of the Riemannian gradient to the product manifold? Since their dimensions are different, how do you perform these projections, since definitions of the projections given on Page 5 cannot be applied? Please check your theoretical/experimental results and code accordingly. Please define exp_{x^(t)_i}(vi) and Exp(U) more precisely. I suppose that they denote exponential maps. How do you initialize x^(0) randomly? The notation is pretty confusing and ambiguous. First, does x belong to an embedded Riemannian manifold P or a point on the graph, which will be embedded? According to equation (2), they are on the graph and they will be embedded. According to Algorithm 1, x^0 belongs to P, which is a Riemannian manifold as defined before. So, if x^(0) belongs to P, then L is already defined from P to R (in input of the Algorithm 1). Thereby, gradient .nabla L(x) is already a Riemannian gradient, not the Euclidean gradient, while you claim that .nabla L(x) is the Euclidean gradient in the text. Overall, Algorithm 1 just performs a projection of Riemannian or Euclidean gradient .nabla L(x) onto a point v_i for each ith individual manifold. Then, each v_i is projected back to a point on an individual component of the product manifold by an exponential map. What do you mean by “sectional curvature, which is a function of a point p and two directions x; y from p”? Are x and y not points on a manifold? You define .xi_G(m;b,c) for curvature estimation for a graph G. However, the goal was to map G to a Riemannian manifold. Then, do you also consider that G is itself a Riemannian manifold, or a submanifold? What is P in the statement “the components of the points in P” in Lemma 2? What is .epsilon in Lemma 2? How do you optimize positive w_i, i=1,2,...,n? What is the “gradient descent” refered to in Lemma 2? Please provide computational complexity and running time of the methods. Please define .mathbb{I}_r. At the third line of the first equation of the proof of Lemma 1, there is no x_2. Is this equation correct? If at least of two of x1, y1, x2 and y2 are linearly dependents, then how does the result of Lemma 1 change? Statements and results given in Lemma 1 are confusing. According to the result, e.g. for K=1, curvature of product manifold of sphere S and Euclidean space E is 1, and that of E and hyperbolic H is 0. Then, could you please explain this result for the product of S, E and H, that is, explain the statement “The last case (one negative, one positive space) follows along the same lines.”? If the curvature of the product manifold is non-negative, then does it mean that the curvature of H is ignored in the computations? What is .gamma more precisely? Is it a distribution or density function? If it is, then what does (.gamma+1)/2 denote? The statements related to use of Algorithm 1 and SGD to optimize equation (2) are confusing. Please explain how you employed them together in detail. Could you please clarify estimation of K_1 and K_2, if they are unknown. More precisely, the following statements are not clear; - “Furthermore, without knowing K1, K2 a priori, an estimate for these curvatures can be found by matching the distribution of sectional curvature from Algorithm 2 to the empirical curvature computed from Algorithm 3. In particular, Algorithm 2 can be used to generate distributions, and K1, K2 can then be found by matching moments.” Please explain how in more detail? What is matching moments? - “we find the distribution via sampling (Algorithm 3) in the calculations for Table 3, before being fed into Algorithm 2 to estimate Ki” How do you estimation K_1 and K_2 using Algorithm 3? - Please define, “random (V)”, “random neighbor m” and “.delta_K/s” used in Algorithm 3 more precisely.