Hydra: A method for strain-minimizing hyperbolic embedding of network- and distance-based data

TL;DR

Hydra is a new hyperbolic embedding method that minimizes hyperbolic strain, recovers points exactly on hyperbolic submanifolds, and offers faster computation with competitive or superior embedding quality.

Contribution

Hydra introduces a mathematically optimal hyperbolic embedding technique with an extended version, hydra+, that improves speed and quality over existing methods.

Findings

Hydra achieves competitive embedding quality with reduced computation time.

Hydra+ outperforms existing methods in both speed and embedding accuracy.

The method guarantees exact recovery for points on hyperbolic submanifolds.

Abstract

We introduce hydra (hyperbolic distance recovery and approximation), a new method for embedding network- or distance-based data into hyperbolic space. We show mathematically that hydra satisfies a certain optimality guarantee: It minimizes the `hyperbolic strain' between original and embedded data points. Moreover, it recovers points exactly, when they are located on a hyperbolic submanifold of the feature space. Testing on real network data we show that the embedding quality of hydra is competitive with existing hyperbolic embedding methods, but achieved at substantially shorter computation time. An extended method, termed hydra+, outperforms existing methods in both computation time and embedding quality.

Figures

Tables

Table 1. Table 1. Networks used for numerical experiments

Network	Description	Source	# Nodes
karate	Social interaction network (‘Zachary’s karate club network’) from [24]	igraphdata [7]	$34$
UKfaculty	Personal friendship network of a UK university faculty from [17]	igraphdata [7]	$81$
opsahl	One-node projection of message Exchange Network from [18]; two isolated nodes have been removed	toreopsahl.com	$897$
facebook	Facebook social circles network from [15]; combined edge sets	snap.stanford.edu	$4039$
collaboration	Co-authorship network from ArXiv submissions to category Hep-Ph (High Energy Physics); largest connected component. From [14]	snap.stanford.edu	$8638$
oregon	Autonomous systems peering information inferred from route-views in Oregon on March 26, 2001. From [13]	snap.stanford.edu	$11174$

Equations

$$\mathbf{x}\circ \mathbf{y}:=x_1{y}_1-\left(x_2{y}_2+\dots +x_{d+1}{y}_{d+1}\right),$$

$${\mathcal{H}}_d=\left\{\mathbf{x}\in {\mathbb{R}}^{1,d}:\mathbf{x}\circ \mathbf{x}=1,x_1>0\right\}.$$

$${\mathrm{d}}_H^{\kappa }(\mathbf{x},\mathbf{y})=\frac{1}{\sqrt{\kappa }}\mathrm{arcosh}\left(\mathbf{x}\circ \mathbf{y}\right),\mathbf{x},\mathbf{y}\in {\mathcal{H}}_d.$$

$$d{s}^2=\frac{1}{\kappa }\left(d\mathbf{x}\circ d\mathbf{x}\right)$$

$${\mathcal{B}}_d:=\left\{\mathbf{z}\in {\mathbb{R}}^d:|\mathbf{z}|<1\right\},$$

$$\xi (\mathbf{x})=\left(\frac{x_2}{1+x_1},\dots ,\frac{x_{d+1}}{1+x_1}\right)$$

$${\mathrm{d}}_B^{\kappa }(\xi ({\mathbf{z}}_1),\xi ({\mathbf{z}}_2))={\mathrm{d}}_H^{\kappa }({\mathbf{z}}_1,{\mathbf{z}}_2),{\mathbf{z}}_1,{\mathbf{z}}_2\in {\mathcal{B}}_d.$$

$$r:=\sqrt{z_1^2+\cdots +z_d^2},\mathbf{u}:=\frac{\mathbf{z}}{r}.$$

$$r$$

$$\mathbf{u}$$

$$z_1=r\cos \theta ,z_2=r\sin \theta .$$

$${\mathrm{d}}_H^{\kappa }({\mathbf{x}}_i,{\mathbf{x}}_j)\approx d_{ij}.$$

$$\mathrm{Stress}({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n{)}^2:=\sum _{i,j=1}^n{\left(d_{ij}-{\mathrm{d}}_H^{\kappa }({\mathbf{x}}_i,{\mathbf{x}}_j)\right)}^2$$

$$\mathrm{Strain}({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n{)}^2:=\sum _{i,j=1}^n{\left(\cosh (\sqrt{\kappa }{d}_{ij})-{\mathbf{x}}_i\circ {\mathbf{x}}_j\right)}^2,$$

$$\underset{{\mathbf{x}}_i\in {\mathbb{R}}^{1,d}}{\min }\sum _{i,j}(\cosh (\sqrt{\kappa }{d}_{ij})-({\mathbf{x}}_i\circ {\mathbf{x}}_j){)}^2.$$

$$r_i={\xi }_r\left(\frac{x_{i1}}{x_{\min }}\right),$$

$$A=[a_{ij}]:=[\cosh (\sqrt{\kappa }{d}_{ij})]$$

$$A=Q\Lambda Q^{\top },$$

$$X:=\left[\sqrt{{\lambda }_1}{\mathbf{q}}_1\sqrt{(-{\lambda }_{n-d+1}{)}^{+}}{\mathbf{q}}_{n-d+1}\cdots \sqrt{(-{\lambda }_n{)}^{+}}{\mathbf{q}}_n\right],$$

$${\mathbf{u}}_i:=\frac{(x_{i2},\dots ,x_{i(d+1)})}{\sqrt{x_{i2}^2+\cdots +x_{i(d+1)}^2}},$$

$$x_{\min }$$

$$r_i$$

$$\mathrm{d}_{B}^{\kappa}\Big{(}(r_{i},\bm{u}_{i}),(r_{j},\bm{u}_{j})\Big{)}=\mathrm{d}_{H}(\bm{x}_{i},\bm{x}_{j})=d_{ij},\qquad i,j=1,\dotsc,n.$$

$${\theta }_i^{\prime }:=\lambda {\theta }_i+(1-\lambda )(\mathrm{ark}({\theta }_i)-1)\frac{2\pi }{n},i=1,\dots ,n.$$

$$X={\left({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n\right)}^{\top }\in {\mathbb{R}}^{n\times (d+1)}$$

$$J=\mathrm{diag}(1,-1,\dots ,-1),$$

$$\underset{X\in {\mathbb{R}}^{n\times (d+1)}}{\min }{\Vert A-X^{\top }JX\Vert }_F^2,$$

$$\mathrm{diag}(X^{\top }JX)=(1,\dots ,1)\text{and}X{\mathbf{e}}_1>0,$$

$$g_{ij}={\mathbf{x}}_i\circ {\mathbf{x}}_j,\forall i,j\in 1,\dots ,n.$$

$$G=XJ{X}^{\top },$$

$$G=\hat{X}\hat{J}\hat{X}^{\top},\quad\text{where}\quad\hat{J}=\operatorname{diag}\Big{(}{+}1,0,\dotsc,0,\underbrace{-1,\dotsc,-1}_{\text{$d$ times}}\Big{)}.$$

Full text

The authors would like to thank Carlo Vittorio Cannistraci and Alessandro Muscoloni for comments and discussions.

1. Introduction

Embeddings of networks and distance-based data into hyperbolic geometry have received substantial interest in recent years. Such embeddings have been used for visualization [23], link prediction [19, 20] and community detection [20, 16]. They offer insight into the tradeoff between popularity and similarity effects in network growth [19] and have interesting implications for routing, network navigability [10, 1] and efficient computation of shortest network paths [25, 5]. Moreover, such embedding methods can be seen as alternatives to classic visualization and dimensionality reduction techniques based on Euclidean geometry, such as principal component analysis or multidimensional scaling. However, the hyperbolic embedding methods as yet proposed in the literature have either been based on specific assumptions about network growth (e.g. [20, 16]), or methods with strong theoretical properties, but requiring costly non-linear numerical optimization procedures (e.g. H-MDS of [23], Rigel of [25] and HyPy of [5]). Here, we introduce hydra (hyperbolic distance recovery and approximation), a novel method for embedding network or distance-based data into hyperbolic space, which has strong mathematical foundations and does not depend on specific assumptions on network growth or structure. At the same time, the method is computationally efficient and based on reduced matrix Eigendecomposition. We show mathematically, that when presented with mutual distances of data points located on a low-dimensional hyperbolic submanifold of the feature space, hydra will recover these points exactly. For general data, the method satisfies a certain optimality property, similar to the strain-minimizing property of multidimensional scaling. Finally, we introduce hydra+, an extension where the result of hydra is used as inital condition for hyperbolic embedding methods based on optimization, such as Rigel/HyPy, substantially improving their efficiency. When tested on real network data, hydra and its variants typically outperform existing hyperbolic embedding methods. All new methods introduced are available in the package hydra [9] for the statistical computing environment R [21].

2. Embeddings into Hyperbolic Space

2.1. Hyperbolic Space

We summarize the key features of the hyperboloid model of hyperbolic geometry (cf. [22, 3]) in dimension $d$. This will provide the mathematical framework in which we formulate our embedding method. To start, we define for $\bm{x},\bm{y}\in\mathbb{R}^{d+1}$ the indefinite inner product

[TABLE]

also called Lorentz product. The real vector space $\mathbb{R}^{d+1}$ equipped with this inner product is called Lorentz space and denoted by $\mathbb{R}^{1,d}$. As nested subsets, it contains the positive Lorentz space $\mathbb{R}^{1,d}_{+}=\left\{\bm{x}\in\mathbb{R}^{1,d}:x_{1}>0\right\}$ and the single-sheet hyperboloid

[TABLE]

The hyperboloid model with curvature $-\kappa$, $(\kappa>0)$, consists of $\mathcal{H}_{d}$ endowed with the hyperbolic distance

[TABLE]

The hyperbolic distance $\mathrm{d}^{\kappa}_{H}$ is a distance on $\mathcal{H}_{d}$ in the usual mathematical sense; in particular it takes only positive values and satisfies the triangle inequality, cf. [22, §3.2]. In fact, it can be shown that $\mathcal{H}_{d}$ endowed with the Riemannian metric tensor

[TABLE]

is a Riemannian manifold and $\mathrm{d}^{\kappa}_{H}(\bm{x},\bm{y})$ is the corresponding Riemannian distance.111That is, $\mathrm{d}^{\kappa}_{H}(\bm{x},\bm{y})$ is the length of the shortest path from $\bm{x}$ to $\bm{y}$, where lengths are measured using the length element $ds=\sqrt{ds^{2}}$. The sectional curvature of this manifold is constant and equal to $-\kappa$, which explains the role of $\kappa$ as curvature parameter. Just as Euclidean space is the canonical model of geometry with zero curvature, hyperbolic space is the canonical model of geometry with negative curvature.

2.2. The Poincaré ball Model

In addition to the hyperboloid model, we introduce the Poincaré ball model of hyperbolic geometry, which is more appealing for visualizations of hyperbolic space and hyperbolic embeddings. In $\mathbb{R}^{d}$, consider the open unit ball

[TABLE]

where $|\bm{z}|=\sqrt{z_{1}^{2}+\dotsm+z_{d}^{2}}$ is the usual Euclidean norm. The stereographic projection $\xi:\mathbb{R}^{1,d}\to\mathbb{R}^{d}$ defined by

[TABLE]

restricts to a bijective mapping from the hyperboloid $\mathcal{H}_{d}$ onto $\mathcal{B}_{d}$ (cf. [22, §4.2]). It transfers the hyperbolic distance from $\mathcal{H}_{d}$ to $\mathcal{B}_{d}$, by setting

[TABLE]

Endowed with this distance, $(\mathcal{B}_{d},\mathrm{d}^{\kappa}_{B})$ is isometric to $(\mathcal{H}_{d},\mathrm{d}^{\kappa}_{H})$ and therefore an equivalent model of hyperbolic geometry.

It will be convenient to parameterize $\mathcal{B}_{d}$ by the radial coordinate $r\in[0,1)$ and the directional coordinate $\bm{u}$ (a unit vector in $\mathbb{R}^{d}$), given by

[TABLE]

An easy calculation shows that the conversion from coordinates in $\mathcal{H}_{d}$ is given by

[TABLE]

In dimension $d=2$, the Poincaré ball becomes the Poincaré disc, and each of its points can be described by the radius $r$ and the unique angle $\theta\in[0,2\pi)$ such that

[TABLE]

2.3. Embedding of Distances and Networks

To formulate the embedding problem, let a symmetric matrix $D=[d_{ij}]\in\mathbb{R}_{\geqslant 0}^{n\times n}$ with zero diagonal be given, which represents the pairwise dissimilarities between some objects $\bm{o}_{1},\dotsc,\bm{o}_{n}$. The basic premise of hyperbolic embedding is that the matrix $D$ can be approximated by a hyperbolic distance matrix $H=[\mathrm{d}^{\kappa}_{H}(\bm{x}_{i},\bm{x}_{j})]$, i.e., that we can find points $\bm{x}_{1},\dotsc,\bm{x}_{n}$ in low-dimensional hyperbolic space $\mathcal{H}_{d}$, such that

[TABLE]

The points $\bm{x}_{1},\dotsc,\bm{x}_{n}$ give a low-dimensional representation in hyperbolic space of the configuration of $\bm{o}_{1},\dotsc,\bm{o}_{n}$ induced by their dissimilarities. In Euclidean space, such approximations are well studied and can be calculated e.g. by multidimensional scaling (MDS), see also Section A.2 and [2].

An important special case is the network embedding problem, where a (unweighted, undirected) graph $G=(V,E)$ is given and $D=[d_{ij}]$ is the graph distance matrix of $G$, i.e., $d_{ij}$ is the length of the shortest path in $G$ from vertex $v_{i}$ to $v_{j}$. In particular for graphs with locally tree-like structure it can be expected that hyperbolic geometry gives a better representation than Euclidean geometry, see e.g [10]. Instead of the shortest-path distance, other dissimilarity measures based on the structure of $G$ can be used, such as the repulsion-attraction (RA) rule or edge-betweenness-centrality (EBC), cf. [16].

2.4. Connection to prior work and innovations

Most existing methods for hyperbolic embedding can be placed into one of two classes: Stress-based methods or network-specific methods.

• •

Stress-based methods aim to solve the embedding problem (7) by minimizing the stress functional

[TABLE]

over all $\bm{x}_{1},\dotsc,\bm{x}_{n}\in\mathcal{H}_{d}$. This minimization problem is a challenging high-dimensional non-convex optimization problem, and methods largely differ in their numerical approach to minimize (8). The H-MDS method proposed in [23] is a gradient descent scheme for minizing (8) based on explicit calculation of the gradient. [4] propose a neural-network-based approach to minimizing (8), while [25] and [5] develop so-called ‘landmark-based’ minimization algorithms (Rigel and HyPy respectively) based on iterative quasi-Newton minimization. Due to the ‘landmark’ heuristic, these methods are able to deal with large-scale instances of (8) and do not require full knowledge of $D$, see [5] for details.

• •

Model-based methods focus on the network embedding problem and rely on underlying assumptions on the generating mechanism of the graph $G$, see e.g. [19] for a model of ‘hyperbolic network growth’. In HyperMap of [20] and the coalescent embedding of [16], the radial coordinate $r_{i}$ of the embedded points in the Poincaré ball model is determined directly from the degree of the vertices $v_{i}$, using the assumption of a power-law relationship. The directional component $\bm{u}_{i}$ of the embedding is then determined by maximizing likelihood in an underlying probabilistic model (cf. [20]) or by applying existing nonlinear dimensionality reduction methods (such as Laplacian Eigenmapping or ISOMAP) to the underlying data (cf. [16]).

Here, our main innovation is to replace the stress functional (8) by the strain functional

[TABLE]

which results from (8) when all distances are transformed by hyperbolic cosine. Furthermore, we introduce a highly efficient method for the minimization of hyperbolic strain, called hydra (hyperbolic distance recovery and approximation). Contrary to stress-minimization, hydra is based on matrix Eigendecomposition, similar to principal component analysis or classic multidimensional scaling.222In fact, the relation between hyperbolic strain- and stress-minimization is similar to the relation between ‘classic’ and ‘metric’ multidimensional scaling in the Euclidean case, cf. [2]. In Theorems 3.1 and 3.2 we show that hydra satisfies important theoretical optimality properties, in particular, it returns a guaranteed global minimum of (9). For instances based on real data, the embedding results of hydra are comparable to those based on pure stress-minimization, even when embedding quality is measured by the stress functional (8); see Section 4 below. This shows, that even when minimization of stress is the final goal, the strain functional (9) is a valuable and useful proxy for stress, as it can be minimized in a highly efficient way. The best results in terms of stress are obtained when strain- and stress-minimization are combined. This is the basis of the hydra+ method, introduced in Section 3.2, where the embedding result of hydra is used as initial condition for a stress-minimization run.

3. A new hyperbolic embedding method

3.1. The hydra algorithm

We introduce the hydra algorithm, displayed as Algorithm 1, which calculates an embedding into the Poincaré ball model of hyperbolic space by efficiently solving the strain-minimization problem

[TABLE]

The algorithm proceeds as follows:

• •

In steps A1 and A2 the strain-minimization problem is solved by means of a matrix Eigendecomposition. These steps return a coordinate matrix $X=[x_{ij}]$, whose rows $\bm{x}_{1},\dotsc,\bm{x}_{n}$ are elements of positive Lorentz space $\mathbb{R}^{1,d}_{+}$ and the optimizers of (10). The optimality of $\bm{x}_{1},\dotsc,\bm{x}_{n}$ is the subject of Theorem 3.2 below.

• •

In steps B1 and B2, the points $\bm{x}_{1},\dotsc,\bm{x}_{n}$ are projected onto the Poincaré ball $\mathcal{B}_{d}$ using the stereographic projection (4) and converted to radial/directional coordinates $(\bm{u}_{i})$ using (5) and (6). No adjustment is necessary for the directional coordinates, which are computed in step B1.

• •

Due to (5), the radial coordinates $(r_{i})$ depend only on the first column $(x_{11},\dotsc,x_{n1})$ of $X$ and can be obtained by applying $\xi_{r}$ elementwise. But $\xi_{r}(x_{1i})$ may be undefined for elements with $x_{1i}\in(0,1)$.333By Theorem 3.2, steps A1 and A2 guarantee that all $x_{1i}$ are positive, but not that they are larger than one. Therefore, $r_{i}$ is calculated in step B2 as

[TABLE]

that is, after rescaling the first column of $X$ by dividing by its smallest element $x_{\min}$.

The key theoretical properties of the hydra algorithm are summarized in the following theorems, whose proofs are given in Appendix A. The first theorem shows that hydra recovers any configuration of points in $d$-dimensional hyperbolic space up to isometry:

Theorem 3.1** (Exact Recovery).**

Let $\bm{a}_{1},\dotsc,\bm{a}_{n}$ be points in hyperbolic $d$-space $\mathcal{H}_{d}$, and let $D=[d_{ij}]=[\mathrm{d}^{\kappa}_{H}(\bm{a}_{i},\bm{a}_{j})]$ be the matrix of their hyperbolic distances with curvature $-\kappa$. Then $\texttt{hydra}(D,d,\kappa)$ recovers the points $\bm{a}_{1},\dotsc,\bm{a}_{n}$ up to isometry. In particular, the rows $\bm{x}_{1},\dotsc,\bm{x}_{n}$ of the matrix $X$ and the points $(r_{1},\bm{u}_{1}),\dotsc,(r_{n},\bm{u}_{n})$ returned by $\texttt{hydra}(D,d,\kappa)$ satisfy

[TABLE]

For applications to real data, exact recovery is an atypical situation. However, hydra enjoys an optimality guarantee for strain minimization, expressed in the following theorem:

Theorem 3.2** (Optimal Approximation).**

The rows $\bm{x}_{1},\dotsc,\bm{x}_{n}$ of the matrix $X$ returned by hydra$(D,d,\kappa)$ are the globally optimal solutions of the strain minimization problem (10). Moreover, the first column of $X$ is strictly positive; equivalently, all $\bm{x}_{i}$ are elements of positive Lorentz space $\mathbb{R}^{1,d}_{+}$.

3.2. Practical guidelines and extensions

While the result of hydra satisfies the theoretical optimality guarantees in Theorem 3.1 and 3.2, it can still be advantageous to adjust the results in order to improve the attractiveness of visualization or the embedding quality in terms of stress (8) (as opposed to strain, which is globally minimal). The so-called equiangular adjustment has been introduced in [16] and can be applied to two-dimensional hyperbolic embeddings. Here we propose a slight modification, which allows to interpolate smoothly between no adjustment and full equiangular adjsutment.

**Equiangular adjustment: **

Let $\lambda\in[0,1]$ be the adjustment parameter and define $\mathrm{ark}(\theta_{i})$ as the angular rank of $\bm{x}_{i}$, i.e. when the embedded points are ordered by increasing angular coordinate $\theta$, then $\mathrm{ark}(\theta_{i})$ is defined as the rank (from $1$ to $n$) of $\bm{x}_{i}$ in this list. The adjusted angular coordinate is then set to

[TABLE]

If $\lambda=0$, no adjustment takes place. If $\lambda=1$ then the angles $\theta^{\prime}_{i}$ are regularly spaced (‘equiangular’) and only the ordering given by $\theta_{i}$ is retained.444This is the equiangular adjustment as proposed in [16]. Values of $\lambda\in(0,1)$ interpolate between these two extremes. We propose a values of $\lambda=1/2$, which typically leads to improvements in both visual appeal and stress value of the embedding; see also method hydra-equi in Figure 1.

**hydra+: **

If minimization of stress is the ultimate objective and strain is used only as a proxy, the result of hydra can be used as an initial condition for a direct minimization of the stress functional (8). This can be seen as a chaining of hydra and HyPy/Rigel [25, 5], where hydra substitutes the random initial condition of HyPy/Rigel. For the minimization of stress, efficient quasi-Newton routines, such as LBFGS [26] can be used and supplied with the explicit gradient of stress, given in [5, Eqs. (3.1),(3.2)].

In terms of efficiency, the following simple improvement can be made to hydra:

**Reduced Eigendecomposition: **

The numerically dominating part of hydra is the Eigendecomposition in (12). Note, however, that in (13) only the single first and the last $d$ Eigenvalues and Eigenvectors of the matrix $A$ are needed. There are efficient numerical routines (see e.g. [12]) to perform such a reduced Eigendecomposition without computing the full Eigendecomposition of $A$. These routines substantially improve efficiency if $n\gg d$ and are used in our implementation of hydra.

Using the reduced Eigendecomposition, we expect the time complexity of hydra to be $\mathcal{O}(n^{\alpha})$ with $\alpha$ slightly above, but close, to $2$, cf. [8, Ch. 55]. For hydra+, the time complexity is harder to estimate, since it is based on iterative minimization of a non-convex objective function. In a single step of LBFGS both the stress functional and its gradient have to be evaluated at a complexity of $\mathcal{O}(n^{2})$. Depending on the number of steps to convergence, we thus also expect a complexity of $\mathcal{O}(n^{\alpha})$, with $\alpha$ strictly larger than $2$. Empirical estimates of $\alpha$ are given in Section 4 below.

3.3. Remarks on strain-minimizing graph embeddings

In the seminal paper [19] it has been argued that the inherent negative curvature in hyperbolic geometry resolves the trade-off between the conflicting attractive forces of popularity and similarity in network growth models. For this reason [19] have proposed to interpret the radial coordinate $r$ in the Poincaré disc as dimension of ‘popularity’ and the angular coordinate $\theta$ as dimension of ‘similarity’. Interestingly, the strain minimization problem (10) and its solution by hydra gives additional mathematical support for this interpretation. More precisely, revisiting Algorithm 1 in the graph embedding context, we observe that:

• •

The radial coordinates $r_{i}$ are determined only from the Perron-Frobenius Eigenvector555The Perron-Frobenius Eigenvector is the Eigenvector associated to the largest Eigenvalue of a positive matrix (i.e. a matrix consisting only of positive entries) and is itself a positive vector, cf. [8, Ch. 10]. $\bm{q}_{1}$ of the matrix $A$. This provides a remarkable connection to the Eigenvector centralities (corresponding to the popularity dimension) of the nodes $v_{i}$, which are determined from the Perron-Frobenius Eigenvector of their adjacency matrix.

• •

The directional coordinates $\bm{u}_{i}$ are determined only through the Eigenvectors $\bm{q}_{n-d+1},\dotsc,\bm{q}_{n}$ (and corresponding Eigenvalues) at the low end of the spectrum of $A$. This provides a remarkable connection to Cheeger’s inequality (cf. [6, Ch. 9]), which shows that the low end of the spectrum of the graph Laplacian matrix encodes the separability of the graph into sparsely connected ‘communities’ (corresponding to the similarity dimension).

We remark that while the connections described above are a first step towards a mathematization of the popularity-similarity paradigm in hyperbolic network geometry, the matrix $A=[\cosh\left(\sqrt{\kappa}d_{ij}\right)]$ is in general neither identical to the adjacency nor to the Laplacian matrix of a given graph, and thus further research into the rigorous mathematical underpinning of these connections is warranted.

4. Numerical Results

4.1. Methods and Data

In our numerical experiments, we evaluate different variants of hydra and compare them to existing hyperbolic embedding methods, using stress as performance criterion. We focus on small to medium sized networks (see Table 1), for which it is still feasible to compute the full distance matrix as input to our methods. Edge weights (when available) were discarded, i.e., all networks were treated as unweighted undirected graphs. This network data was used as input for the following methods:

**hydra: **

The hydra method (without equi-angular adjustment) as described in Algorithm 1

**hydra-equi: **

The hydra method with equiangular adjustment $\lambda=0.5$, as described in Section 3.2

**hydra+: **

The hydra+ method as described in Section 3.2 and using the result of hydra-equi as initial condition.

**HyPy/Rigel: **

The HyPy algorithm from [5], which is based on Rigel from [25]. Both methods are based on direct minimization of the stress functional (8). Landmark selection, as proposed in [5] was not implemented, since it serves to reduce runtime and memory use for large networks, but is not expected to improve embedding results. As in [5], the initial condition for minimization was chosen at random and we repeated the embedding $100$ times.

**CE-LE: **

The coalescent embedding (CE) using Laplacian Eigenmapping (LE) as dimension-reduction method, full equiangular adjustment and repulsion-attraction (RA) pre-weighting; see [16] for details. Among the methods developed in [16], this was the best performing method to invert the PSO generating model of [20] for hyperbolic networks.

For the methods hydra, hydra-equi, hydra+ and HyPy/Rigel we used our own implementations in R, which are available in the R-package hydra, [9]. For CE-LE we used the MATLAB implementation of the methods of [16] available from github.666https://github.com/biomedical-cybernetics/coalescent_embedding The stress-optimization in hydra+ and HyPy/Rigel was performed using the LBFGS method (see [26]) as implemented in the R-function optim and using the analytic form of the gradient of the stress functional (8) from [5]. Note that all methods except CE-LE use the shortest-path matrix as input dissimilarities; CE-LE uses repulsion-attraction (RA) weights as input dissimilarities, see [16]. For all methods hyperbolic curvature was fixed to $-\kappa=-1$ and we embed into dimension $d=2$.

4.2. Results and Discussion

Results on embedding quality (measured by stress) for all networks and methods (except CE-LE) are shown in Figure 1. Note that stress values are normalized by using the average result of HyPy/Rigel as a reference. This facilitates the comparison of results between different networks. As HyPy/Rigel depends on randomized initial conditions we indicate the $5\%$- and $95\%$-quantiles (over 100 runs) in addition to its average result. The stress-values of the embeddings produced by CE-LE were substantially larger (by a factor from $13$ to $26$) than the reference method and we have therefore excluded this method from the plot and from the further analysis of computation times.

Computation times of the different methods is shown in Figure 2 in doubly logarithmic coordinates. As should be expected, the methods split into two groups with computation times for hydra and hydra-equi being shorter than for hydra+ and HyPy/Rigel by two orders of magnitude. The seemingly small gap between hydra+ and HyPy/Rigel still corresponds to a difference of about 50% in runtime. Based on the discussion in Section 3.2 we have added regression lines to estimate the exponent $\alpha$ in the conjectured complexity $\mathcal{O}(n^{\alpha})$. To avoid clutter, regression lines are only shown for hydra-equi and hydra+; the estimates for all methods are $\alpha\approx 2.0$ for hydra, $\alpha\approx 2.1$ for hydra-equi and $\alpha\approx 2.3$ for both hydra+ and HyPy/Rigel. As setup costs seem to dominate the computation times for the smallest network, we have excluded it from the regression analysis.

Finally, an exemplary embedding result produced by the hydra-equi method for the facebook network is shown in Figure 3. Nodes are placed into the Poincaré disc model of hyperbolic geometry according to their embedding coordinates (the full disc is indicated in grey) and a random subsample of links is drawn as hyperbolic geodesics. Visually, the embedding conforms well with the popularity-vs-similarity paradigm of [19] for hyperbolic networks: Nodes with a function as hubs between communities or individuals (popularity dimension) are placed closer to the center of the hyperbolic disc. Communities are identified along the angular coordinate (similarity dimension) with the effective distance between communities indicated by angular separation.

Summarizing our numerical experiments, we conclude the following:

• •

In general, the strain-minimization performed by hydra seems to be a good proxy for stress-minimization, but is faster by a factor of 100 or more in comparison to stress-minimization from a random initial condition (HyPy/Rigel). Note that hydra also eliminated the uncertainty associated with the randomized nature of HyPy/Rigel, which can lead to large variations in embedding quality in some instances (e.g., the facebook network).

• •

The simple equi-angular adjustment performed in hydra-equi consistently improves embedding quality in terms of stress at negligible numerical costs. The returned embeddings outperform HyPy/Rigel for two networks (karate, facebook) and are competitive for all others, with a largest observed difference of 23% in terms of stress.

• •

Using the results of hydra as a starting value for stress-minimization, instead of a random initial condition, i.e., replacing HyPy/Rigel by hydra+ reduces computation time by approx. 30% - 50% and leads to better (average) embedding quality in all cases. The reduction in stress is considerable in the facebook network, where stress is reduced by approx. 40%.

• •

The CE-LE method, based on the PSO network growth model of [20], is not competitive with the other methods in terms of embedding quality. This suggests that the structure of the real networks that we have considered deviates from the theoretical growth model (PSO-model) of [20] upon which CE-LE is build.

As a next step, we plan to make strain-minimizing hyperbolic embedding methods feasible for large and very large networks. For such networks the computational complexity of $\mathcal{O}(n^{\alpha})$ (with $\alpha>2$) of the proposed methods, but also of the graph distance calculation itself, are prohibitive. For this reason, heuristics such as the landmark heuristic of [5] will have to be adapted to strain-minimizing embedding methods.

Appendix A Theoretical Results

To prove the theoretical properties of the hydra method, it is convenient to reformulate the strain minimization problem (10) in matrix form. To this end, let $D=[d_{ij}]$ be the given dissimilarity matrix, set $A=[\cosh(\sqrt{\kappa}\,d_{ij})]$ and write

[TABLE]

for the coordinate matrix of some points $\bm{x}_{1},\dotsc,\bm{x}_{n}$ in $\mathbb{R}^{d+1}$. Finally, let $J$ be the $(d+1)\times(d+1)$ diagonal matrix

[TABLE]

cf. [22, §3.1]. The strain minimization problem (10) can now be written in compact form as

[TABLE]

where $\left\|.\right\|_{F}$ denotes the Frobenius norm. Imposing the constraint that all $\bm{x}_{i}$ are elements of the hyperboloid $\mathcal{H}_{d}$ is equivalent to requiring that

[TABLE]

where $\bm{e}_{1}$ is the first standard unit vector. In particular, the first condition guarantees $\bm{x}_{i}\circ\bm{x}_{i}=1$, and the second one selects the upper sheet of the two-sheet hyperboloid thus described.

A.1. Hyperbolic strain minimization and exact recovery

For a real symmetric matrix $A$, denote by $n_{+}(A)$ and $n_{-}(A)$ the number of positive and negative Eigenvalues of $A$. The following Lemma characterizes matrices that can be written as inner product matrices (‘Gram matrices’) with respect to the Lorentz product (1):

Lemma A.1**.**

Let $G=[g_{ij}]\in\mathbb{R}_{\geqslant 0}^{n\times n}$ be positive and symmetric, and let $d\leq n-1$. The following are equivalent

1. a)

$G$* satisfies $n_{+}(G)=1$ and $n_{-}(G)\leq d$.* 2. b)

$G$* is a ‘Lorentzian Gram matrix’, i.e., there exist $\bm{x}_{1},\dotsc,\bm{x}_{n}$ in $\mathbb{R}^{1,d}$, such that*

[TABLE] 3. c)

There exists $X\in\mathbb{R}^{n\times(d+1)}$, such that

[TABLE]

where $J$ is given by (15).

In addition,

• •

The first column of $X$ is positive if and only if $\bm{x}_{1},\dotsc,\bm{x}_{n}$ are in the positive Lorentz space $\mathbb{R}^{1,d}_{+}$;

• •

The points $\bm{x}_{1},\dotsc,\bm{x}_{n}$ are in $\mathcal{H}_{d}$ if and only if $\operatorname{diag}(G)=(1,\dotsc,1)$ and the first column of $X$ is positive.

Proof.

The equivalence of (b) and (c) follows directly from the definition of the Lorentz product in (1). Next, we show that (c) implies (a): From [11, Ch. 10.3] it follows from (17) that $n_{+}(G)\leq n_{+}(J)=1$ and $n_{-}(G)\leq n_{-}(J)=d$. But $G$ is a positive matrix and Perron’s theorem (cf. [11, Ch. 16]) guarantees that its leading Eigenvalue is positive, i.e., $n_{+}(G)\geq 1$, and we conclude (a). To show that (a) implies (c), assume first that $n_{-}(G)=d$. By Sylvester’s law of inertia, there exists a decomposition

[TABLE]

This decomposition can be reduced to (17), by simply dropping all rows and columns containing only zeroes from $\hat{J}$ and by also dropping the corresponding columns from $\hat{X}$. If $n_{-}(G)=d^{\prime}<d$, the same procedure yields a decomposition with $X$ of dimension $n\times(d^{\prime}+1)$ and $J$ of dimension $(d^{\prime}+1)\times(d^{\prime}+1)$. Padding $X$ with zero columns and $J$’s diagonal with $-1$s, (17) also follows in this case.

The additional statements follow directly from the following observations: The first column of $X$ contains exactly the first coordinate of all points $\bm{x}_{1},\dotsc,\bm{x}_{n}$. If the first coordinate of a point $\bm{x}$ is positive, it is an element of positive Lorentz space and vice versa. The diagonal of $G$ contains the values $\bm{x}_{i}\circ\bm{x}_{i},i=1,\dotsc,n$. If $\bm{x}_{i}\circ\bm{x}_{i}=1$ and $\bm{x}_{i}\in\mathbb{R}^{1,d}_{+}$ then $\bm{x}_{i}$ is an element of the hyperboloid $\mathcal{H}_{d}$ and vice versa. ∎

Proof of Theorem 3.2.

Let $A=[a_{ij}]=[\cosh(\sqrt{\kappa}\,d_{ij})]$ and let $B=[b_{ij}]$ be another symmetric matrix in $\mathbb{R}^{n\times n}$. Let $(\lambda_{i}(A))_{i=1,\dotsc,n}$ and $(\lambda_{i}(B))_{i=1,\dotsc,n}$ be their Eigenvalues in descending order, and denote by $\left\|.\right\|_{F}$ the Frobenius norm. By a result of Wielandt-Hoffmann, cf. [11, Ch. 10, Thm. 18],

[TABLE]

Assume now that $B$ is a ‘Lorentzian Gram matrix’ with elements given by

[TABLE]

for some $\bm{b}_{1},\dotsc,\bm{b}_{n}\in\mathbb{R}^{1,d}$. By Lemma A.1 this implies that $n_{+}(B)=1$ and $n_{-}(B)\leq d$. Hence all Eigenvalues of $B$ with index $2,\dotsc,n-d$ are zero, and we obtain

[TABLE]

For the first summand on the right hand side we have the trivial lower bound [math]. In the last sum, all $\lambda_{i}(B)$ are negative or zero, and hence, for any $i=(n-d+1),\dotsc,n$, we can estimate

[TABLE]

which is the same as $(\lambda_{i}(A)^{+})^{2}$. Together, we obtain that

[TABLE]

Denote by $A=Q\Lambda_{A}Q^{\top}$ the Eigendecomposition of $A$ with $\Lambda_{A}=\operatorname{diag}(\lambda_{1}(A),\dotsc,\lambda_{n}(A))$. Let $X$ be the matrix returned by hydra($D$, $d$, $\kappa$) and $\bm{x}_{1},\dotsc,\bm{x}_{n}$ the rows of $X$. By (13) the associated Lorentzian Gram matrix $G=XJX^{\top}$ has the Eigendecomposition $G=Q\Lambda_{G}Q^{\top}$ with

[TABLE]

Using the unitary invariance of the Frobenius norm and the trivial identity $x-(-x)^{+}=x^{+}$, we obtain

[TABLE]

This shows that setting $\bm{b}_{i}:=\bm{x}_{i}$ for all $i\in 1,\dotsc,n$ achieves equality in (19) and hence that the points $\bm{x}_{i}$ minimize (16). ∎

Proof of Theorem 3.1.

Let $D=[d_{ij}]$ be the hyperbolic distance matrix of $\bm{a}_{1},\dotsc,\bm{a}_{n}$ in $\mathcal{H}^{d}$. Then $A=[a_{ij}]=[\cosh(\sqrt{\kappa}d_{ij})]$ is the associated Lorentzian Gram matrix with elements

[TABLE]

By Lemma A.1 $A$ satisfies $n_{+}(A)=1$ and $n_{-}(A)\leq d$, i.e. the Eigenvalues of $A$ satsify $\lambda_{i}(A)=0$ for $i=2,\dotsc,n-d$ and $\lambda_{i}(A)\leq 0$ for $i=n-d+1,\dotsc,n$. Hence, it follows from (20) that $\sum_{i,j}(a_{ij}-\bm{x}_{i}\circ\bm{x}_{j})^{2}=0$ or, equivalently, that

[TABLE]

for all $i,j\in 1,\dotsc,n$. Applying $\cosh(\sqrt{\kappa}\,\cdot)$ to both sides, we see that

[TABLE]

and hence that $(\bm{x}_{i})$ and $(\bm{a}_{i})$ are isometric. ∎

A.2. Comparison to classic multidimensional scaling

In several aspects, the hydra method can be seen as the ‘hyperbolic analogue’ of classic multidimensional scaling (MDS), cf. [2], which is based on Euclidean geometry. Below, we summarize the classical MDS method and point out parallels to (and differences from) hydra. Classical MDS also takes a matrix $D=[\mathrm{d}_{ij}]\in\mathbb{R}_{\geqslant 0}^{n\times n}$ with zero diagonal as input. Using the centering matrix $C=I-\frac{1}{n}\bm{1}\in\mathbb{R}^{n\times n}$, where $\bm{1}$ denotes a matrix of ones of matching dimension, the ‘doubly centered’ matrix

[TABLE]

is derived from $D$, and its Eigendecomposition

[TABLE]

computed. Again, $\Lambda$ is the diagonal matrix of the Eigenvalues $\lambda_{1}\geq\dotsm\geq\lambda_{n}$ and the columns of $Q$ are the Eigenvectors $\bm{q}_{1},\dotsc,\bm{q}_{n}$. MDS then returns the (Euclidean) coordinate matrix

[TABLE]

whose rows $\bm{x}_{i}$ are interpreted as points in Euclidean space $\mathbb{R}^{d}$. This coordinate matrix $X$ solves the strain minimization problem

[TABLE]

cf. [2, Ch. 12]. Moreover, if the input matrix $D$ is a matrix of squared Euclidean distances, i.e., $\mathrm{d}_{ij}=|\bm{x}_{i}-\bm{x}_{j}|^{2}$ then MDS recovers the points $\bm{x}_{i}$ exactly (up to Euclidean isometry). Note that $X^{\top}X$ appearing above is the Gram matrix of the points $\bm{x}_{1},\dotsc,\bm{x}_{n}$, i.e. the matrix of their scalar products $\bm{x}_{i}^{\top}\bm{x}_{j}$, whereas the matrix $X^{\top}JX$ in (16) is the ‘Lorentzian Gram matrix’ of the Lorentz products $\bm{x}_{i}\circ\bm{x}_{j}$.