Link prediction with hyperbolic geometry

TL;DR

This paper systematically investigates hyperbolic geometry's effectiveness in link prediction, showing it excels particularly in challenging scenarios with difficult-to-predict links, and introduces a new coordinate inference method for better accuracy.

Contribution

The paper develops a new hyperbolic coordinate inference method and demonstrates its effectiveness in improving link prediction for complex, hard-to-predict links in networks.

Findings

Hyperbolic link prediction is competitive for easy links but excels for difficult links.

A new coordinate inference method improves hyperbolic model accuracy.

Hyperbolic geometry is especially useful for predicting hard-to-predict links.

Abstract

Link prediction is a paradigmatic problem in network science with a variety of applications. In latent space network models this problem boils down to ranking pairs of nodes in the order of increasing latent distances between them. The network model with hyperbolic latent spaces has a number of attractive properties suggesting it must be a powerful tool to predict links, but the past work in this direction reported mixed results. Here we perform systematic investigation of the utility of latent hyperbolic geometry for link prediction in networks. We first show that some measures of link prediction accuracy are extremely sensitive with respect to inaccuracies in the inference of latent hyperbolic coordinates of nodes, so that we develop a new coordinate inference method that maximizes the accuracy of such inference. Applying this method to synthetic and real networks, we then find that while there exists a multitude of competitive methods to predict obvious easy-to-predict links, among which hyperbolic link prediction is rarely the best but often competitive, it is the best, often by far, when the task is to predict less obvious missing links that are really hard to predict. These links include missing links in incomplete networks with large fractions of missing links, missing links between nodes that do not have any common neighbors, and missing links between dissimilar nodes at large latent distances. Overall these results suggest that the harder a specific link prediction task is, the more seriously one should consider using hyperbolic geometry.

Tables

Table 1. Table 1: The summary of the results in Sec.  III : hyperlink ’s measures of accuracy of link prediction in RHGs with known node coordinates as functions of the parameters in Sec.  III .

Parameter	AUC	AUPR, Precision
Exponent $\gamma \in (2,3)$	$\approx \mathrm{const}$	$\approx \mathrm{const}$
Temperature $T\in (0,1)$	$\approx \mathrm{const}$	decreasing
Fraction of missing links $1-q$	increasing	increasing
Noise amplitude $a$	$\approx \mathrm{const}$	decreasing

Table 2. Table 2: The summary of the results in Sec.  IV : hyperlink ’s measures of accuracy of link prediction in RHGs with inferred coordinates and in real networks, as well as those for nonlocal links, compared to other methods. The parameters are the same as in Table  1 .

Scenario	AUC	AUPR, Precision
RHGs with inferred coordinates	Winner	Winner if $T$, $\gamma$, or $1-q$ is large
Real networks	Winner/ Runner-up	The more competitive, the larger the $1-q$
Nonlocal links in RHGs and real networks	Winner/ Runner-up	Winner/ Runner-up

Table 3. Table 3: Basic properties of the considered real networks. N 𝑁 N is the number of top nodes; E 𝐸 E is the number of edges; k ¯ ¯ 𝑘 \overline{k} is the average degree; γ 𝛾 \gamma is the degree distribution exponent, which we estimated using methods from Ref.  Voitalov  et al. ( 2019 ) ; c ¯ ¯ 𝑐 \overline{c} is the average degree-dependent clustering coefficient; and T 𝑇 T is the corresponding RHG temperature.

Network name	$N$	$E$	$\overline{k}$	$\gamma$	$\overline{c}$	$T$
Internet	$6,474$	$13,234$	$4.09$	$2.1$	$0.51$	$0.7$
Metabolic network	$2,732$	$4,040$	$2.96$	$2.9$	$0.29$	$0.6$
PGP web of trust	$14,138$	$160,080$	$22.65$	$2.1$	$0.66$	$0.8$

Equations

$$\tilde{p}\left(x_{ij}\right)=\frac{p\left(x_{ij}\right)\left(1-q_{ij}\right)}{1-p\left(x_{ij}\right)+p\left(x_{ij}\right)\left(1-q_{ij}\right)}.$$

$$\tilde{p}\left(x_{ij}\right)=\frac{\left(1-q\right)p\left(x_{ij}\right)}{1-qp\left(x_{ij}\right)}$$

$$\cosh x=\cosh r\cosh r^{\prime }-\sinh r\sinh r^{\prime }\cos \Delta \theta ,$$

$${\theta }_i\leftarrow \rho (\theta )$$

$$r_i\leftarrow \rho (r)$$

$$p\left(x_{ij}\right)=\frac{1}{1+e^{\frac{x_{ij}-R}{2T}}}.$$

$$\mathcal{L}\left(\{{\mathbf{x}}_i\}|a_{ij},\mathcal{P},q\right)=\frac{\mathcal{L}\left(a_{ij}|\{{\mathbf{x}}_i\},\mathcal{P},q\right)\mathrm{Prob}({\mathbf{x}}_i)}{\mathcal{L}\left(a_{ij}|\mathcal{P},q\right)},$$

$$\mathcal{L}\left(a_{ij}|\{{\mathbf{x}}_i\},\mathcal{P},q\right)=\prod _{i<j}{\left[\tilde{p}\left(x_{ij}\right)\right]}^{a_{ij}}{\left[1-\tilde{p}\left(x_{ij}\right)\right]}^{1-a_{ij}},$$

$$\tilde{p}\left(x\right)\equiv qp(x),$$

$$\mathrm{Prob}({\mathbf{x}}_i)=\frac{1}{{\left(2\pi \right)}^N}\prod _{i=1}^N\rho (r_i),$$

$$\ln \mathcal{L}\left(\{{\mathbf{x}}_i\}|a_{ij},\mathcal{P},q\right)=K+\sum _{i=1}^N\ln \rho (r_i)+$$

$$+$$

$${\hat{\theta }}_i$$

$$\mathrm{Precision}=\frac{|{\Omega }_R\cap {\Omega }_M({\lambda }^{\ast })|}{|{\Omega }_R|},$$

$$\mathrm{TP}\left(\lambda \right)$$

$$\mathrm{FN}\left(\lambda \right)$$

$$\mathrm{FP}\left(\lambda \right)$$

$$\mathrm{TN}\left(\lambda \right)$$

$$\mathrm{tpr}(\lambda )$$

$$\mathrm{fnr}(\lambda )$$

$$\mathrm{fpr}(\lambda )$$

$$\mathrm{tnr}(\lambda )$$

$$\mathrm{AUC}={\int }_0^1\mathrm{tpr}(\lambda ){\mathrm{fpr}}^{\prime }(\lambda )\mathrm{d}\lambda .$$

$$\mathrm{pr}(\lambda )$$

$$\mathrm{rc}(\lambda )$$

$$\mathrm{AUPR}={\int }_0^1\mathrm{pr}(\lambda ){\mathrm{rc}}^{\prime }(\lambda )\mathrm{d}\lambda .$$

$$\mathrm{AUC}={\int }_0^{2R}\mathrm{tpr}(x){\mathrm{fpr}}^{\prime }(x)\mathrm{d}x,$$

$$n(x)=\frac{4{\alpha }^2}{\pi {\left(2\alpha -1\right)}^2}{e}^{x/2-R}$$

$$1-AUC\{\begin{array}{ll}\sim N^{-1}& \text{if T∈[0,21) },\\ =\mathcal{O}\left(\frac{\ln N}{N}\right)& \text{if T=21 },\\ =\mathcal{O}\left(N^{1-\frac{1}{T}}\right)& \text{if T∈(21,1). }\end{array}$$

$$\mathrm{AUPR}={\int }_0^{2R}\mathrm{pr}(x){\mathrm{rc}}^{\prime }(x)\mathrm{d}x,$$

$${\hat{\theta }}_i$$

$$X_i$$

$$\mathrm{tp}(x|a)\sim \{\begin{array}{ll}{a}^{1-2\gamma }& \text{if x≤R},\\ a^{1-2\gamma }\left(R+2\ln \frac{a}{2}\right)& \text{if x>R},\end{array}$$

$$\mathrm{AUPR}(a)\sim a^{2-4\gamma }{\left(R+2\ln \frac{a}{2}\right)}^2.$$

$$s_{ij}^{CN}=|\Gamma (i)\cap \Gamma (j)|.$$

Full text

I Introduction

Link prediction is a paradigmatic example of forecasting network dynamics Peng et al. (2015); Lü and Zhou (2011); Menon and Elkan (2011); Peixoto (2018), with diverse applications including the reconstruction of networks based on partial data Marchette and Priebe (2008); Guimerà and Sales-Pardo (2009); Kim and Leskovec (2011) and prediction of future social ties Peng et al. (2015); Adamic and Adar (2003); Newman and Clauset (2016), protein interactions von Mering et al. (2002); Yu et al. (2008); Kovács et al. (2019), and user ratings in recommender systems Zhou et al. (2007); Lü et al. (2012); Bobadilla et al. (2013); Schafer et al. (1999).

Latent space network models Gilbert (1961); McFarland and Brown (1973); McPherson et al. (2001); Krioukov et al. (2010); Newman and Peixoto (2015) offer an intuitive and simple approach to link prediction. In these models, network nodes are points in a latent space, while connections are established with probabilities that decrease with latent distances between nodes. Latent distances model similarity between nodes, and the main idea behind these models is to model homophily: more similar nodes are more likely to be linked. Link prediction then reduces to ranking unconnected node pairs in the order of increasing latent distances between them: the closer the two unlinked nodes in the latent space, the higher the probability of a missing link Brew and Salter-Townshend (2010); Zhu et al. (2016); Peixoto (2018); García-Pérez et al. (2020).

Among many latent space models considered in literature, only the one that assumes that the latent space is hyperbolic reproduces sparsity, self-similarity, scale-free degree distribution, strong clustering, the small-world property, and community structure Serrano et al. (2008); Krioukov et al. (2010); Papadopoulos et al. (2012); Zuev et al. (2015). All these properties are often observed in many real networks Lazega et al. (2006); Newman (2010); Barabási and Pósfai (2016), and hyperbolic geometry captures them all. In addition, the hyperbolic network model is likely to be the simplest or parsimonious with respect to these properties, as in some of its limiting regimes it has been proven to be statistically unbiased, satisfying the maximum entropy principle van der Hoorn et al. (2018); Krioukov (2016).

Given the combination of these attractive properties, one could naturally expect that the hyperbolic latent space model must be a powerful tool in link prediction. Yet the previous studies on this subject reported mixed results Serrano et al. (2012); Papadopoulos et al. (2015a, b); Muscoloni and Cannistraci (2018a, b); García-Pérez et al. (2020).

Here we perform systematic investigation of the efficiency of link prediction using latent hyperbolic geometry. We organize the presentation of the results as follows.

In Sec. II we recall the definitions of the hyperbolic latent space network model, which for short we call random hyperbolic graphs (RHGs), and outline the basic idea behind link prediction based on this model. We also recall the definitions of the main measures of link prediction accuracy—AUC (area under receiver-operating characteristic), AUPR (area under precision-recall curve), and Precision—and discuss what these measures actually measure: while AUPR cares mostly about most obvious easy-to-predict missing links, AUC puts more weight on less obvious and harder-to-predict missing links between more dissimilar nodes, albeit with the cost of not caring that much about false positives.

Our main results are then given in Secs. III and IV. In Sec. III we calculate analytically the AUC and AUPR on RHGs with known hyperbolic coordinates of all nodes. That is, the same coordinates are used both to generate RHGs and to predict missing links in them, an ideal situation yielding the upper bound for the link prediction accuracy using hyperbolic geometry. To understand the robustness of link prediction in the case where coordinates are inferred (Sec. IV), so that they are not equal exactly to the true coordinates, we add uniform noise to the true coordinates, and analyze the AUC, AUPR, and Precision as functions of the noise amplitude to find that: (1) AUC is not that sensitive to noise, but (2) AUPR and Precision decrease quickly as noise grows. The latter result implies that the AUPR and Precision scores of link prediction using hyperbolic geometry in real networks can be high only if node coordinates are inferred with sufficiently high accuracy. This is because the most likely missing links candidates are those between similar nodes at small hyperbolic distances, which are most sensitive to coordinate inaccuracies.

To predict missing links in networks with unknown coordinates one first needs to infer these coordinates. Motivated by the results in Sec. III calling for high-accuracy coordinate inference, and given that no existing hyperbolic coordinate inference algorithm is sufficiently accurate, in Sec. IV we develop an alternative one, which we call the hyperlink embedder, the focus of which is on high precision in coordinate inference. We present its overview in Sec. IV, while all the details are delegated to Appendix F, where we also compare it to some existing inference algorithms to show that its accuracy is indeed higher. A software package implementing the hyperlink embedder is hosted by the Bitbucket repository cod .

We then apply the hyperlink embedder to a collection of RHGs with “forgotten” coordinates, and to real networks, calling the overall link prediction procedure the hyperlink method, and comparing it to a representative collection of other link prediction methods.

Section V contains both high-level (Tables 1 and 2) and more detailed summaries of all the results. The results are definitely not that the hyperlink or any other method is a clear winner in all the considered scenarios according to all the considered link prediction accuracy measures. We discuss what methods are strong in what scenarios. The hyperlink appears to be the strongest in the most difficult link prediction tasks. That is, the more challenging a particular link prediction task/scenario, the better off is the hyperlink compared to other methods. We conclude the paper with an outline of open problems at the end of Sec. V.

These results emphasize that the hyperlink is definitely not the universally best link prediction method, which simply cannot exist as was recently shown in Vallès-Català et al. (2018); Ghasemian et al. (2019, 2020); García-Pérez et al. (2020). That is, there can exist no one size fits all solution for the link prediction problem. Different methods are good at predicting different types of links. Therefore, as far as a particular link prediction method is concerned, the best one can do is to document what particular link prediction scenarios the method is good at; that is, what types of links the method is good at predicting, which is exactly the subject of this paper.

II Methods

We begin the exposition by discussing the latent geometric link prediction framework and the null model that we utilize to predict missing links.

II.1 Link prediction with latent geometry

Link prediction with hyperbolic geometry is a two-step procedure. First, one needs to infer node coordinates in the hyperbolic space and calculate hyperbolic distances between node pairs. This coordinate inference procedure is often referred to as network mapping or embedding. The second step of the procedure is to identify most likely missing link candidates. This subsection focuses on the second step of this procedure. The technical details of the null geometric model and the network mapping algorithm constituting the first step, are provided in Secs. II.2 and II.3 and Appendix (F). We refer to the network mapping algorithm and the entire hyperbolic link prediction framework as the hyperlink embedder and the hyperlink, respectively.

The latent geometric link prediction framework is applicable to all latent geometric models, where connections are established independently with decreasing connection probability function $p(x)$. Intuitively, the smaller the latent distance between two nodes, the higher the probability of a link between them. Then, if two nodes located close to each other in the latent space are not connected, it is likely that there is a missing link between them.

Specifically, consider a latent geometric model where nodes are assigned positions $\{\mathbf{x}_{i}\}$ in a certain latent space $\mathcal{M}$, and every node pair $\{ij\}$ is connected with probability $p_{ij}=p\left(x_{ij}\right)$, where $x_{ij}=d\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)$ is the latent distance between the nodes, and $p:\mathbb{R}_{+}\to[0,1]$ is the decreasing connection probability function specified by the model. After all connections are established, some links are removed with probabilities $1-q_{ij}$. These pairs of nodes are referred to as missing links.

Any unconnected node pair $\{ij\}$ in the resulting network is either not connected in the network formation process, or connected in the network formation and later removed with probability $1-q_{ij}$. Therefore, the probability for an unconnected pair of nodes $\{ij\}$ separated by $x_{ij}$ to be a missing link, is

[TABLE]

In the particular case of a decreasing connection probability function $p(x)$ and the random link removal process, $q_{ij}=q$

[TABLE]

is the decreasing function of $x_{ij}$ for any $q>0$. Thus, the most probable candidates for missing links are indeed unconnected node pairs located at small latent distances, as stated, and the latent geometric link prediction algorithm only needs to rank unconnected node pairs in the increasing order of latent distance between them.

It is important to note, however, that this approach is only guaranteed to work in the case the links are removed uniformly at random. In the general case, missing link probabilities in Eq. (1) depend both on latent distances $\{x_{ij}\}$ and missing link rates $\{1-q_{ij}\}$ and further information on the nature of $\{q_{ij}\}$ is needed to rank missing link candidates properly.

II.2 Random hyperbolic graphs

While the latent geometric framework described above is applicable to all latent space models, in our paper we use the RHG as a null model for link prediction.

RHGs have been extensively studied in the literature Krioukov et al. (2009, 2010); Boguñá et al. (2010); Papadopoulos et al. (2015a, b); Kitsak et al. (2017); Aldecoa et al. (2015); García-Pérez et al. (2018) and have been shown to reproduce common properties of many real networks including heterogeneous distributions of node degrees, strong clustering, as well as community structure Krioukov et al. (2010); Zuev et al. (2015); Muscoloni and Cannistraci (2018c).

The latent space of the RHG model is the two-dimensional hyperbolic disk of constant negative curvature $K=-1$ and radius $R$. The hyperbolic distance $x$ between any two points in the hyperbolic disk is given by the hyperbolic law of cosines:

[TABLE]

where $(r,\theta)$ and $(r^{\prime},\theta^{\prime})$ are the hyperbolic coordinates of the two points within the disk and $\Delta\theta=\pi-|\pi-|\theta-\theta^{\prime}||$ is the angle between them.

The RHG has three parameters — hyperbolic disk radius $R>0$, temperature $T\in[0,1)$ and node density parameter $\alpha>1/2$ — and is defined as follows:

Draw node coordinates $\{r_{i},\theta_{i}\}$, $i=1,2,\ldots,N,$ from probability density functions:

[TABLE] 2. 2.

Compute distances $\{x_{ij}\}$ between all node pairs using Eq. (3). 3. 3.

Connect node pairs with probability

[TABLE]

We summarize basic RHG properties in Appendices C: parameter $\alpha$ controls the exponent $\gamma=2\alpha+1$ of the power-law degree distribution, while clustering is a decreasing function of temperature $T$ approaching zero in the $N\to\infty$ limit as $T\to 1$. In this limit, clustering is zero for any $T\geq 1$.

II.3 HYPERLINK embedder in a nutshell

To infer hyperbolic coordinates of nodes in a given network with random links removed, we aim to find the set of node coordinates $\{\mathbf{x}_{i}\}\equiv\{(r_{i},\theta_{i})\}$, $i=1,2,\ldots,N$, that maximize the posterior probability $\mathcal{L}\left(\{\mathbf{x}_{i}\}|a_{ij},\mathcal{P},q\right)$, which is the probability density function of coordinates $\{\mathbf{x}_{i}\}$ in an RHG with adjacency matrix $a_{ij}$, parameters $\mathcal{P}=\{\alpha,T,R\}$, and link removal probability $1-q$. By the Bayes’ rule this probability is

[TABLE]

where $\mathcal{L}\left(a_{ij}|\{\mathbf{x}_{i}\},\mathcal{P},q\right)$ is the likelihood that network $a_{ij}$ is generated as an RHG with subsequent random link removal with probability $1-q$, ${\rm Prob}(\mathbf{x}_{i})$ is the prior probability of node coordinates generated by the RHG, and $\mathcal{L}\left(a_{ij}|\mathcal{P},q\right)$ is the probability that the network has been generated as the RHG with random link removal. Since node pairs are connected independently, this likelihood is

[TABLE]

where $\tilde{p}\left(x_{ij}\right)$ is the effective connection probability in the RHG generation process with subsequent random link removal:

[TABLE]

where $p(x)$ is the RHG connection probability function in Eq. (6). Finally, in RHGs node coordinates $\{\mathbf{x}_{i}\}\equiv\{r_{i},\theta_{i}\}$, and the prior probability is given by

[TABLE]

where $\rho(r_{i})$ is as in Eq. (5).

The hyperlink embedder aims to find node coordinates $\hat{\mathbf{x}}_{i}$ that maximize the likelihood $\mathcal{L}\left(\{\mathbf{x}_{i}\}|a_{ij},\mathcal{P},q\right)$, or equivalently its logarithm

[TABLE]

where constant $K$ absorbs all terms independent of $\{\mathbf{x}_{i}\}$.

Similar to other maximum-likelihood estimation (MLE) based embedders Boguñá et al. (2010); Papadopoulos et al. (2015a, b); García-Pérez et al. (2019), node coordinates $\hat{\mathbf{x}}_{i}$ are computed iteratively: starting with initial random coordinate configuration, the hyperlink embedder updates node coordinates at each iteration step to increase $\ln\mathcal{L}\left(\{\mathbf{x}_{i}\}|a_{ij}\right)$ and stops when we arrive to a stable configuration. One feature of the hyperlink embedder which is different from other MLE-based embedders is that at each iteration step $\ell$ the embedder adds synthetic noise of variable magnitude $a(\ell)$ to angular node coordinates:

[TABLE]

where $X_{i}$ is a random number drawn from the uniform distribution on the circle $[0,2\pi]$. These coordinate perturbations allow the hyperlink embedder to avoid getting trapped for long time in local maxima of the log-likelihood function and to find (nearly) optimal solutions much faster, thus increasing the coordinate inference accuracy given the same amount of computational resources (see Appendix F for details).

II.4 Link prediction accuracy

We evaluate the accuracy of the hyperlink as well as other link prediction methods through random link removal experiments. To this end, we first remove existing links uniformly at random with probability $1-q$ from the network of interest $G$. We refer to the remaining network as the pruned network and denote it by $\tilde{G}$. We refer to removed links as missing links and denote them by $\Omega_{R}$. The set of remaining links in $\tilde{G}$ is referred to as $\Omega_{E}$.

To test the link prediction method of interest we compute likelihood scores for all unconnected node pairs in $\tilde{G}$, $\overline{\Omega}_{E}$, which include both missing links $\Omega_{R}$ and true nonlinks $\Omega_{N}$, so that $\overline{\Omega}_{E}=\Omega_{R}\cup\Omega_{N}$. We then rely on these scores to rank unconnected node pairs in the decreasing order of missing link likelihood and refer to them as missing link candidates. We denote the fraction $\lambda\in[0,1]$ of most likely missing link candidates as set $\Omega_{M}(\lambda)$. In the case $\lambda=0$, $\Omega_{M}(\lambda)$ is the empty set, while in the $\lambda=1$ case $\Omega_{M}(\lambda)=\Omega_{R}\cup\Omega_{N}=\overline{\Omega}_{E}$.

In the case the exact number of missing links is known, the most direct way to assess link prediction accuracy is to consider the same number of the most likely missing link candidates and evaluate its intersection with the set of missing links. This metric is known as Precision and is formally defined as

[TABLE]

where fraction $\lambda^{*}=1-q$ is chosen such that $|\Omega_{M}(\lambda^{*})|=|\Omega_{R}|$. The Precision score is bounded by [math] and $1$ with the upper bound corresponding to the ideal link predictor ranking all missing links in $\Omega_{R}$ higher than nonlinks in $\Omega_{N}$.

In practical circumstances, however, the exact number of missing links is often unknown. Further, depending on the application, one might be interested to minimize the number of false positives in the prediction set, possibly by the expense of false negatives, or vice versa, minimize the number of false negatives by the expense of false positives. One example of the former case where one is interested to minimize the number of false positives, i.e., good citizens misclassified as criminals, is the criminal justice system. This example is in contrast to cancer screening, where the number of false negatives, or not-identified cancer cases, should be minimized. In both cases one is interested to explore the performance of the link predictor for a range of $\Omega_{M}(\lambda)$ sizes.

A number of link prediction metrics have been developed to this end with the receiver operating characteristic (ROC) and the precision-recall (PR) being the most popular.

To formally introduce ROC and PR curves we first define the confusion matrix. The latter consists of four values — the numbers of true positives (TP), false positives (FP), false negatives (FN), and true negatives (TN), Fig. 1 — and is extensively used in statistical classification problems. Link prediction is not a genuine classification problem since one is only interested to predict links and not their absence. Nonlink node pairs are predicted implicitly as unconnected node pairs that are not part of $\Omega_{R}$.

In the context of link prediction, the number of true positives is the number of correctly identified missing links from $\Omega_{M}(\lambda)$, Eq. (14). The number of false negatives is the remaining number of missing links that are not part of the $\Omega_{M}(\lambda)$, Eq. (15). The number of false positives is the number of missing link candidates in $\Omega_{M}(\lambda)$ that are not correctly identified, Eq. (16). Finally, the number of true negatives is the number of unconnected node pairs that are neither true positives nor false positives nor false negatives [see Eq. (17) and Fig. 1]:

[TABLE]

Since network sizes vary, it is common to normalize confusion matrix elements, obtaining true positive, false positive, false negative, and true negative rates, formally defined as

[TABLE]

An ROC statistics or curve is defined as the parametric plot of the true positive rate ${\rm tpr}(\lambda)$ as a function of the false positive rate ${\rm fpr}(\lambda)$ obtained by varying the fraction of considered link candidates $\lambda\in[0,1]$. The ideal predictor is expected to rank all node pairs corresponding to missing links, $\Omega_{R}$, higher than nonlinks, $\Omega_{N}$, resulting in unit true positive rate and zero false positive rate for $\lambda=1-q$, ${\rm tpr}(1-q)=1$, ${\rm fpr}(1-q)=0$. The corresponding ROC curve of the ideal predictor is thus a rectangle going through the upper left corner $(0,1)$ of the ROC space. A fully random link predictor, on the other hand, will guess missing links at random from $\overline{\Omega}_{E}$ and is expected to yield equal true positive and false positive rates, ${\rm tpr}(\lambda)={\rm fpr}(\lambda)$ for all $\lambda$ values, resulting in the diagonal ROC curve, Fig. 2a.

The standard way to quantify ROC-based prediction accuracy is through the AUC:

[TABLE]

AUC values vary in between [math] and $1$ with ${\rm AUC}=0.5$ corresponding to a fully random predictor and ${\rm AUC}=1.0$ corresponding to the perfect predictor.

The AUC score can be interpreted as the probability that a randomly chosen missing link is assigned a higher link prediction score than a randomly chosen unconnected node pair. ROC curves are easy to read and interpret, which is arguably the basic reason behind their popularity.

At the same time, there is a growing consensus that ROC curves and corresponding AUC scores are insensitive in class imbalance problems, where the size of the positives is disproportional to that of the negatives Davis and Goadrich (2006). Link prediction in sparse networks is one example of class imbalance. Here the number of missing links is of the order of $N$ and is significantly smaller than the number of nonlinks, which is of the order of $N^{2}$. Intuitively, in this situation the ${\rm tpr}(\lambda)$ rate grows much faster than the false positive rate since the latter is normalized by $|\Omega_{N}|$ and, as a result, most ROC curves tend to be substantially above the random baseline, yielding AUC scores close to $1.0$, regardless of the link prediction method.

An alternative to the ROC curve is the PR characteristic, defined as the parametric plot of the precision rate ${\rm pr}(\lambda)$ as a function of the recall rate ${\rm rc}(\lambda)$ obtained by varying $\lambda\in[0,1]$, where the two rates are defined by

[TABLE]

That is, the recall rate is identical to the true positive rate, while the precision rate differs from the latter by a different normalization — to the number of predicted links versus the number of removed links.

In the case of an ideal predictor, the precision rate is maximized, ${\rm pr}(\lambda)=1.0$ for $\lambda\leq 1-q$, while the recall is growing from ${\rm rc}(0)=0$ to ${\rm rc}(1-q)=1$, resulting in the rectangular PR curve going through the upper right corner $(1,1)$ of the PR space. A fully random predictor, on the other hand, maintains constant precision rate equal to the ratio of the number of true missing links to the total number of unconnected node pairs, ${\rm pr}^{rand}(\lambda)=\frac{|\Omega_{R}|}{|\Omega_{R}|+|\Omega_{N}|}$ for all $\lambda$ values, Fig. 2b. The standard metric quantifying PR-based prediction accuracy is the AUPR:

[TABLE]

AUPR values vary between $\frac{|\Omega_{R}|}{|\Omega_{R}|+|\Omega_{N}|}$ and $1$ with the unit score corresponding to the ideal predictor. In the case of sparse networks $\Omega_{R}\ll\Omega_{N}$, leading to ${\rm AUPR}\ll 1$ in the case of a random predictor. Unlike ROC curves, PR characteristics do not directly depend on the number of true negatives and, as a result, do not suffer from the class imbalance problem in case of sparse networks.

II.5 AUC versus AUPR

While both AUC and AUPR quantify link prediction accuracy, they tend to weigh missing link candidates differently. AUPR scores tends to emphasize highly ranked missing links candidates, i.e., those corresponding to small $\lambda$ values. AUC scores, on the other hand, put more weight on missing links candidates corresponding to larger $\lambda$ values.

Indeed, AUPR averages precision rate ${\rm pr(\lambda)}$ over the recall rate ${\rm rc(\lambda)}$. Since the recall rate is given by ${\rm rc}(\lambda)=\frac{|\Omega_{R}\cap\Omega_{M}(\lambda)|}{|\Omega_{R}|}$, Eq. (24), good link predictors tend to reach ${\rm rc(\lambda)}=1$ values when the size of missing link candidates set $\Omega_{M}(\lambda)$ becomes comparable to that of $\Omega_{R}$: $|\Omega_{M}(\lambda)|\approx|\Omega_{R}|\ll|\Omega_{N}|$. The latter inequality holds in the case of sparse networks, where the number of links is much smaller than the number of nonlinks. Thus, $|\Omega_{M}(\lambda)|\ll|\Omega_{N}|$, which corresponds to $\lambda\ll 1$ values. Thus, AUPR link prediction scores are dominated by small $\lambda$ fractions, i.e., by the most likely and, typically, most obvious missing link candidates in $\Omega_{M}$.

AUC scores, on the other hand, average true positive rate ${\rm tpr}(\lambda)$ over false positive rate ${\rm fpr}(\lambda)$. The latter takes large values when $|\Omega_{M}(\lambda)|$ becomes comparable to $|\Omega_{N}|$, i.e. for $\lambda$ values close to $1$. AUC scores, thus, are emphasizing not only easy-to-predict links at small $\lambda$ values but also harder to predict links in $\Omega_{M}$ at intermediate and large $\lambda$ values.

In summary, AUC and AUPR scores complement each other by weighing missing link candidates in $\Omega_{M}$ differently. Thus, in our paper we compute both metrics to obtain a comprehensive view on the utility of hyperbolic geometry in link prediction. In addition to AUPR and AUC scores, we also compute Precision scores, which are the scores to use if the number of missing links is known exactly, although such knowledge is rarely the case in practice.

III Link Prediction with Known Coordinates

Before investigating link prediction accuracy in real networks, we conduct link prediction experiments with RHGs with known coordinates. In doing so we pursue several goals. The RHGs provide the upper bound for link prediction accuracy of the hyperlink if the same node coordinates are used both for the graph construction and for link prediction García-Pérez et al. (2020), so that we want to quantify this upper bound. Second, we want to measure link prediction accuracies of other methods, listed in Appendix B, and compare them to that of the hyperlink. Establishing these results provides a baseline for interpreting link prediction results on real networks. To achieve these goals, we first calculate analytically the AUC and AUPR in RHGs with known coordinates and with coordinates disturbed by noise of varying magnitude. The latter result allows us to quantify in a controlled environment the level of coordinate inaccuracy beyond which the hyperlink becomes essentially impuissant.

We start with the analysis of hyperlink accuracy in the case of randomly missing links in RHGs. After the generation of an RHG we visit each of its links and remove it with probability $1-q$, arriving at a pruned network. We then rank missing link candidates using distances between all unconnected node pairs calculated with coordinates from which the network was originally generated.

As seen in Fig. 3, the predictive power of the hyperlink is maximized as $T\to 0$ and decreases as $T$ increases. This result is expected. In the $T\to 0$ limit the RHG is deterministic since the connection probability in Eq. (6) becomes the Heaviside step function, $p(x)\to\Theta(R-x)$. As a result, all node pairs with $x<R$ are connected and other node pairs are not. Then, an unconnected pair of nodes at distance $x<R$ is guaranteed to be a true positive and all unconnected pairs at $x\geq R$ are true negatives. As $T$ increases, connections are allowed at distances $x>R$ with increasing probability and, as a result, underlying geometry plays a smaller role in the formation of links, explaining the decreasing link prediction accuracy as a function of $T$, as quantified by all scores in Fig. 3.

Even though all scores, AUC, AUPR and Precision, are decreasing functions of $T$, they behave differently. AUC scores remain constant in the $T\in\left(0,\frac{1}{2}\right)$ interval and then exhibit a slow decay to ${\rm AUC=0.95}$ at $T=0.9$. AUPR and Precision scores, on the other hand, decrease rapidly in the entire testing interval of $T\in[0,0.9]$ from ${\rm AUPR}=1$ (${\rm Precision}=1$) at $T=0$ to ${\rm AUPR}=0.34$ (${\rm AUPR}=0.29$) at $T=0.9$.

We can predict these results analytically as we explain next.

III.1 AUC

The AUC score in RHGs is

[TABLE]

where ${\rm tpr}(x)$ and ${\rm fpr}(x)$ are, respectively, distance-dependent true positive and false-positive rates among node pairs separated by distances up to $x$. As seen from Fig. 4a, the true positive rate grows exponentially for $x<R$ and saturates to ${\rm tpr}(x)=1$ as $x$ approaches $2R$. This observation is easy to predict analytically. Let $n(x)$ be the distribution of hyperbolic distances $x$ between node pairs in the RHG. It follows from the results in Alanis-Lobato and Andrade-Navarro (2016) that $n(x)$ can be approximated as

[TABLE]

for $\alpha>\frac{1}{2}$. To be more specific, we note that $R$ in the RHG is a function of network size $N$, given by Eq. (47) and the approximation in Eq. (27) holds in the large $N$ limit for any $x=cR$, where constant $c\in(0,2)$, $\lim_{N\to\infty}\frac{n(x)}{n^{\rm true}(x)}=1$. Henceforth, we say $f(x)\approx g(x)$ if $\lim_{N\to\infty}\frac{f(x)}{g(x)}=1$, and, more generally, $f(x)\sim g(x)$ if $\lim_{N\to\infty}\frac{f(x)}{g(x)}=K\neq 1$.

The connection probability $p(x)$ is close to unity for $x<R$, so that, the number of true positives for $x<R$ grows proportional to the number of node pairs $N(x)\equiv\int_{0}^{x}n(y){\rm d}y$ in the hyperbolic disk, ${\rm tpr}(x)\sim N(x)\sim e^{\frac{x}{2}}$ for $x<R$. In the $x>R$ regime connection probability $p(x)$ decays exponentially as $p(x)\sim e^{-\frac{x}{2T}}$ faster than the exponential growth of $n(x)$, leading to the saturation of ${\rm tpr}(x)=1$, Sec. D.

The false positive rate remains negligible for $x<R$ and grows exponentially for $x\in(R,2R)$, Fig. 4b. We explain this observation using similar arguments. Since $p(x)$ is close to unity for $x<R$, and all unconnected node pairs with $x<R$ are almost guaranteed to be true positives, the false positive rate is negligible for $x<R$. In the $x>R$ regime, $p(x)\sim e^{-\frac{x}{2T}}$, and the number of unconnected node pairs is proportional to $N(x)$, resulting in ${\rm fpr}(x)\sim e^{\frac{x}{2}}$ for $x>R$, Sec. D.

Taken together, ${\rm tpr}(x)$ and ${\rm fpr}(x)$ rates provide a qualitative explanation for nearly perfect ${\rm AUC}$ scores observed in Fig. 3a. The false positive rate ${\rm fpr}(x)$ takes large values only when $x$ approaches $2R$. At the same time, as $x$ approaches $2R$ the ${\rm tpr}(x)$ approaches 1.

Supporting this rough estimation, our more detailed analytical calculations in Sec. D show that the AUC scores for RHGs with known coordinates converge to $1$ in the large $N$ limit as

[TABLE]

III.2 AUPR

To calculate the AUPR score we need to calculate the distance-dependent precision and recall rates ${\rm pr}(x)$ and ${\rm rc}(x)$ because

[TABLE]

Since $p(x)$ is close to $1$ for $x<R$, all unconnected node pairs at $x<R$ are true positives, resulting in ${\rm pr}(x)=1$ [see Fig. 4c and Sec. D]. The precision rate decays exponentially for $x>R$ since the true positive rate ${\rm tpr}(x)$ approaches $1$ for $x>R$, while the number of unconnected node pairs $N_{d}(x)$ grows exponentially, $\binom{N}{2}\int_{0}^{x}n(y)\left[1-qp(y)\right]{\rm d}y\sim e^{\frac{x}{2}}$ [see Fig. 4c and Sec. D].

The dependence of ${\rm AUPR}$ on $T$ arises from the recall function or its derivative, $rc^{\prime}(x)$, quantifying the expected distance-dependent link density and, consequently, the density of missing links. As $T$ increases, the missing links are more likely to be located at larger distances, Fig. 4d, where precision ${\rm pr}(x)$ is smaller, resulting in lower AUPR scores, consistent with our observations in Fig. 3.

We also note that the AUPR score depends weakly on the node density parameter $\alpha$ and consequently on the degree distribution exponent $\gamma=2\alpha+1$. Indeed, the precision and recall rates depend on $\alpha$ only via the node pair distribution $n(x)$, Sec. D, which depends on $\alpha$ only in subleading terms, as shown in Ref. Alanis-Lobato and Andrade-Navarro (2016).

III.3 Coordinate uncertainty and link prediction accuracy

While the hyperlink provides the upper bound for link prediction on RHGs, it is important to note that its accuracy is comparable to that of other link prediction methods, in particular, resource allocation (RA), Adamic-Adar (AA), and stochastic block models SBM(d,n), Fig. 3. This observation motivates the question: How accurately does one need to infer node coordinates to ensure the superior performance of the hyperlink?

To answer this question we analyze the impact of node coordinate uncertainty on the hyperlink accuracy. To this end, we add synthetic noise to original angular node coordinates, while keeping radial node coordinates unchanged:

[TABLE]

where $a>0$ is the noise amplitude. The effects of synthetic noise on the hyperlink accuracy are depicted in Fig. 5. Our results indicate that AUPR and Precision scores, Fig. 5b,c, decrease rapidly as a function of noise amplitude, while AUC scores remain largely unchanged even at $a>1~{}{\rm radians}$ values.

To better understand the effects of noise on link prediction accuracy we juxtapose hyperlink prediction results to those of the RA method, which is one of its leading competitors according to Fig. 3. We show RA accuracy with dashed lines of matching color in Fig. 5. Consistent with our earlier observations we find that hyperlink AUC scores are robust to noise, preserving its leading ranking among other link prediction methods, Fig. 5a.

In contrast, as quantified by AUPR and Precision scores, the hyperlink is superior to the RA method only if coordinate uncertainty is sufficiently small. The maximum tolerable noise amplitude value $a_{c}$ increases as $T$ increases [see the inset of Fig. 5b,c]. While noise amplitude $a$ does not exceed $10^{-2}$ radians in the case of $T=0.1$, the noise tolerance in the case of $T=0.9$ is significantly higher, $a_{c}\approx 0.5$ radians, suggesting, somewhat surprisingly, that the hyperlink is better off on networks characterized by larger $T$ values or, equivalently, smaller clustering coefficient.

Qualitatively, the observed fast degradation of the AUPR and Precision scores is due to the sensitivity of the hyperbolic distance to the angular distance between the nodes $\Delta\theta$. It follows from Eq. (3) that even a small change in $\Delta\theta$ may significantly change the corresponding hyperbolic distance, adversely affecting the ranking of missing link candidates at small distances $x$, Appendix E. Since AUPR and Precision emphasize link prediction accuracy of most likely candidates, proper ranking of unconnected node pairs at small $x$ values is crucial. AUC scores, on the other hand, place more emphasis on less obvious link candidates and are less affected by coordinate uncertainty. We find that the uniform synthetic noise adversely affects distance dependent true positive rate ${\rm tp}(x|a)$, which scales as

[TABLE]

see Appendix E, leading to

[TABLE]

The robustness of the AUC scores to synthetic noise in RHGs can be qualitatively explained by the fact that AUC scores emphasize the prediction of missing links at large $x$ distances. Large hyperbolic distances are affected by synthetic noise to a lesser extent than small hyperbolic distances. This effect follows directly from Eq. (3) and can be observed in Fig. 14a, displaying the saturation of ${\rm tp}(x|a)\to 1$ as $x$ approaches $2R$, regardless of noise amplitude $a$.

Our conclusions in this section are different for AUC and AUPR/Precision metrics.

The AUPR and Precision metrics emphasize prediction of the most likely missing link candidates and are highly sensitive to the accuracy of node coordinate inference. Synthetic noise added to original node coordinates smears hyperbolic distances among missing link candidates, adversely affecting the hyperlink accuracy. Our results suggest that one needs to maximize the accuracy of the network mapping in order to efficiently predict missing links. We also find that as temperature $T$ increases, the performance of other link prediction methods, as measured by AUPR and Precision, decreases faster than that of the hyperlink, suggesting that the latter has a competitive advantage on networks characterized by large $T$ values.

AUC scores, on the other hand, emphasize less obvious link candidates that correspond to node pairs at larger hyperbolic distances. Since larger hyperbolic distances are affected by coordinate uncertainty to a lesser extent, the AUC scores of the hyperlink are robust to synthetic noise, suggesting that hyperlink is capable of predicting less obvious missing links even under less accurate mapping conditions.

IV Link Prediction with Inferred Coordinates

In this section we build upon our results obtained in the previous section to analyze the hyperlink accuracy on networks with unknown node coordinates. We first conduct systematic analysis of hyperlink accuracy on RHGs with unknown node coordinates and then apply hyperlink to several real networks. In both cases network coordinates are unknown and in order to predict missing links we first infer node coordinates by mapping networks of interest to the two-dimensional hyperbolic disk. To this end, we developed a mapping algorithm, which is tailored to the link prediction problem. This algorithm is referred to as the hyperlink* embedder* and is fully described in Appendix F.

IV.1 Tests on RHGs with inferred coordinates

To evaluate the hyperlink accuracy on RHGs with unknown node coordinates we perform the following experiments. After generating an RHG we remove a fraction of existing missing links. As before, each existing link is removed with probability $1-q$. Occasionally, after links are removed, the remaining network splits into several components. If this is the case, we limit our consideration to the largest connected component of the pruned network. We refer to the resulting connected component of the pruned network as the training network. To predict missing links we erase our knowledge of the true node coordinates and then infer node coordinates by mapping the training network to the hyperbolic disk using the hyperlink embedder (see Appendix F for details on the mapping procedure). After the mapping is complete, we use the inferred node coordinates to calculate distances between all unconnected node pairs in the training network and rank these pairs in the increasing order of distance.

Figures 6, 7, and 8 show the results for the AUPR, Precision, and AUC scores, respectively. Each panel in these figures is a heatmap, aggregating the link prediction accuracy scores for RHGs with different $\gamma\in[2.1,2.9]$ and $T\in[0.1,0.9]$ values, which we change with an increment of $0.1$ each. We compare the hyperlink to the RA method, which is its leading competitor in these experiments [cf. Figs. 3(a) and  3(b)].

The results for the AUPR and Precision scores are similar. Quantified by these scores, the hyperlink accuracy is nearly independent of the degree distribution exponent $\gamma$, and at the same time decreases rapidly as temperature $T$ increases, see Figs. (6,7)a,d,g. This observation is consistent with our theoretical analysis in Sec. III, where we establish that AUPR scores decrease as $T$ increases and do not strongly depend on $\gamma$.

Even though RA performs similar to hyperlink, Figs. (6,7)b,e,h, we note that RA is more accurate at lower $T$ values and less accurate than hyperlink for higher $T$ values. To obtain the direct comparison of the two methods we plot the difference between their AUPR (Precision) scores in Figs. 6c,f,i(7c,f,i). In agreement with our theoretical considerations in Fig. 5, we find that the hyperlink is superior to RA in the region of $\gamma$-$T$ phase space corresponding to higher $T$ values; these regions are denoted with dashed lines in Figs. (6,7)c,f,i.

Compared to RA, the hyperlink yields better link prediction accuracy for larger fractions of missing links. In the case $1-q=0.1$, for instance, hyperlink is better than RA in a small upper right corner region of the $\gamma$-$T$ phase space, Fig. (6,7)c. On the other hand, in the case $50\%$ of links are missing, $1-q=0.5$, the hyperlink outperforms RA for the majority of $\gamma$-$T$ values with the exception of smallest, $T=0.1$, and largest, $T=0.9$, temperature values, Fig. (6,7)i.

The better, compared to RA, performance of the hyperlink in Fig. (6,7)i is the result of two effects. On one hand, the hyperlink accuracy appears to increase as $1-q$ increases. This effect is consistent with a recent observation in Ref. García-Pérez et al. (2020) that the upper bound of link predictability in edge-independent graphs increases with $1-q$. On the other hand, as $1-q$ increases, the accuracy of RA decreases. RA, as well as other similarity-based methods, e.g., RA, Cannistraci resource allocation (CRA), AA, common neighbors (CN), and Jaccard’s index (JC), predict missing links based on the similarity of node neighborhoods, e.g, the number of common neighbors; the higher the similarity the higher the probability of a missing link, Appendix B. Neighborhood similarities are local measures, reflecting network structure in the network-based vicinity of the node pair of interest, and ignoring the structure of the remaining network. The larger the fraction of missing links, the smaller the fraction of links in the training network and, as a result, the poorer the link prediction results. While this is true for all link prediction methods, the similarity-based methods are the ones that suffer most. Since links are established independently in RHGs, and each link is removed with probability $p=1-q$, the number of common neighbors between any node pairs on average decreases proportionally to $p^{2}$. All extensive RHG properties, on the other hand, depend on $p$ linearly. hyperlink as a global method uses the structure of the entire network to map it, so that it is less sensitive to network incompleteness.

An attractive feature of a global method is that it is capable of predicting nonlocal missing links, i.e., links between node pairs with no common neighbors. To quantify hyperlink accuracy for nonlocal links we consider the subset of nonlocal links within the set of links removed with probability $1-q=0.5$, Fig. (6,7)j, which comprise from $20\%$ (for $\gamma=2.1$, $T=0.9$) to $86\%$ (for $\gamma=2.9$, $T=0.1$) of all removed links.

Similarity-based methods, RA, AA, CN, and JC, cannot predict nonlocal missing links since corresponding node pairs have no common neighbors at all, and, consequently, have zero similarity. Therefore, in nonlocal link prediction experiments we compare hyperlink to the structural perturbation method (SPM) index, which is a global method and the leading competitor to hyperlink for nonlocal links. As seen in Figs. (6,7)k,l, the SPM index yields substantially lower link prediction accuracy than hyperlink for all the considered values of $\gamma$ and $T$.

Overall, we observe that according to the AUPR and Precision scores hyperlink’s competitive advantage is higher the more incomplete the network is, and the hyperlink is particularly strong in prediction of nonlocal links.

According to AUC scores, the hyperlink offers superior link prediction accuracy across the entire $\gamma$-$T$ parameter space, surpassing its leading competitors—RA for all links, and SPM for nonlocal links, Fig. 8. This result is again consistent with our calculations in Sec. III showing that RHG-based AUC scores are robust with respect to coordinate uncertainty.

IV.2 Tests on real networks

Finally, we apply the hyperlink to real networks: the network of human metabolism Ma and Zeng (2003), the Internet at the autonomous system level rou , and the Pretty-Good-Privacy (PGP) web of trust ope . Basic properties of these networks as well as the data curation steps are documented in Appendix A.

Our link prediction experiments on real networks are performed identically to those on RHGs with inferred coordinates, and the results are shown in Figs. 9-11.

According to AUPR and Precision metrics, the hyperlink offers competitive performance in random link removal experiments, Figs. 9-11(a-f), but, at the same time, is not the most accurate. We do note that the relative performance of the hyperlink is better in cases of higher missing link rate, $1-q=0.5$, which is consistent with our results in Sec. III.

We also note that the hyperlink offers superior performance in prediction of nonlocal links where it is either the winner or runner-up, with the SBM methods being its leading competitors, Figs. 9-11(g-i). This observation comes in sharp contrast with nearly random performance of similarity based methods, RA, AA, CN, JC, and CRA, in nonlocal link prediction.

In contrast to AUPR-based rankings where the hyperlink is rarely the most accurate method, it is either the winner or runner-up in all the experiments according to the AUC metric, in agreement with all the AUC-related results above. In particular, it is the winner in predicting nonlocal links in the most challenging human metabolic network. This network is the most challenging because it is the sparsest and has the lowest clustering, Appendix A, thus providing the least amount of local information for link prediction.

V Summary, Discussion, and Conclusion

Tables 1 and 2 summarize the results in Secs. III and IV, respectively. We see that when it comes to predicting obvious missing links that are easy to predict employing hyperbolic geometry may be an overkill. In fact, one should consider using much simpler local methods instead of any global ones, according to the AUPR or Precision results presented here. This is because according to these results the local methods appear to be nearly as good as the global ones at predicting easy links. In particular, the hyperlink method cannot be the best at predicting the most obvious missing links because such links are the links between closest nodes in the latent hyperbolic space, and to rank them exactly at the top of the disconnected node pair list one has to infer the coordinates nearly exactly, Sec. III.

However, if the task is to identify missing links that are really hard to predict, then this is the situation where one should consider using global methods in general and hyperbolic geometry in particular. The most striking example is the prediction of missing links between the nodes that do not share any common neighbors. Here the hyperlink is either the winner or runner-up to the SBM methods, according to all the AUC, AUPR, and Precision measures, in all the considered real and synthetic networks. It is not surprising that local methods do a poor job in predicting such links—they are simply not designed to do so. In contrast, the hyperlink, SBM, and SPM are global methods that base their decisions on the global structure of the whole network, which helps enormously to predict nonlocal and other hard-to-predict links. The SBM and SPM methods were reported to outperform a vast collection of other methods Guimerà and Sales-Pardo (2009); Peixoto (2018); Lü et al. (2015). Here we see that the hyperlink outperforms even these powerful methods in many cases. In particular, the hyperlink is the winner according to all the scores in the most challenging considered case, which is nonlocal links in the sparsest lowest-clustering network of human metabolic reactions.

We also see that according to the AUC measure, the hyperlink is either the winner or runner-up in all the considered situations. This is because the AUC does not care that much about false positives, and hyperlink achieves (nearly) the best balance between the true and false positive rates by finding missing links between highly dissimilar nodes located at large distances in the latent hyperbolic space.

We have also shown that the hyperlink is better off the weaker the clustering (the higher the $T$) is, and the larger the fraction of missing links $1-q$ in RHGs with inferred coordinates. This does not mean that hyperlink’s link prediction accuracy scores are getting better in these more difficult conditions; its scores do degrade. But the speeds of the degradation of these that the other methods experience are higher than hyperlink’s.

Our results also resolve the controversy among earlier reports on link prediction using hyperbolic geometry Serrano et al. (2012); Papadopoulos et al. (2015a, b); Muscoloni and Cannistraci (2018a, b); García-Pérez et al. (2020). These reports approached link prediction using different measures of link prediction accuracy. To reiterate, if applied to sparse networks, the AUPR emphasizes the prediction of a small fraction of the most likely missing links and, as a result, is extremely sensitive to inaccuracies in the node coordinate inference. On the other hand, the AUC is more robust to coordinate uncertainties as it emphasizes the prediction of less likely missing links between dissimilar nodes at large latent distances.

To maximize hyperlink’s link prediction accuracy, we have developed a hyperbolic network mapping method, the hyperlink embedder, that maximizes the accuracy of coordinate inference. Its accuracy comes at the computational complexity cost of $O\left(n^{2}\right)$. While faster methods for hyperbolic mapping have been developed recently Bläsius et al. (2016); Wang et al. (2016); Alanis-Lobato et al. (2016a, b); Muscoloni et al. (2017), an optimal balance between the accuracy and speed of hyperbolic mapping is still to be found. Ideally, it would be highly desirable to have a method that would be as accurate as at least the hyperlink embedder, and that would run in $O\left(n\right)$ time.

We emphasize that link prediction using latent hyperbolic geometry is expected to yield good results only if this geometry is there in a given network. That is, the network structure must be consistent with the existence of this geometry. It is well known that RHGs are characterized by sparsity, self-similarity, scale-free degree distributions, and strong clustering, meaning that these properties are necessary conditions for hyperbolic geometry presence. It is also well known that many real networks do possess these properties as well. The results in Krioukov (2016) suggest that clustering is also a sufficient condition for network geometricity, but these results apply only to homogeneous large-world networks, and ignore coordinate entropy. That is, in theory, the detailed sufficient conditions for the presence of latent hyperbolic geometry are currently unknown, remaining a subject of ongoing research. Experimentally it is known however that random hyperbolic graphs are good descriptors of the structure of many real networks. In particular, we are not aware of any other model capable of reproducing self-similarity of real networks, a highly nontrivial property Serrano et al. (2008). As far as the more standard structural properties of real networks are concerned, the adequacy of hyperbolic geometry to model them has been documented many times, as early as in Papadopoulos et al. (2012). Here we report similar results in Fig. 12.

Overall, it appears that the harder a specific link prediction task, the better the hyperlink is at this task. Yet the hyperlink is not always the winner even at such hard tasks. In particular, in application to real networks it is often a close runner-up to the stochastic block model methods. These results are consistent with the findings in Ref. Faqeeh et al. (2018), where the RHG and SBM were compared across a variety of properties. In the SBM the connection probability has a block structure, while in RHGs it is a function of the latent distance, Fig. 13a,b. Clearly neither model can pretend to describe the connection probability in real networks exactly—at least because the RHGs have no communities, while the SBM has no clustering in the large-network limit. In view of the results in Refs. Colomer-de Simón et al. (2013); Faqeeh et al. (2018), the connection probability in real networks is likely to be some nontrivial mixture of the two pictures, Fig. 13c, with geometry appearing as a mesoscopic structure gluing community blocks together. In short, the RHGs and SBM are complementary models capturing different aspects of the structure of real networks, and the link prediction accuracy of a model-based method depends on how prominent and prevalent the model’s features are in a given real network.

VI Acknowledgements

We thank R. Aldecoa, F. Papadopoulos, C. V. Cannistraci, M. Á. Serrano, T. Peixoto, and A. Clauset for useful discussions and suggestions. This work was supported by Army Research Office (ARO) Grant No. W911NF-17-1-0491 and National Science Foundation (NSF) Grant No. IIS-1741355. M. Kitsak was additionally supported by the NExTWORKx project.

Appendix A Real Networks

A.1 Metabolic network

The metabolic network is based on the dataset of metabolic interactions of $107$ organisms constructed by Ma and Zeng Ma and Zeng (2003). The original network is bipartite and consists of metabolites (top domain) connected to chemical reactions (bottom domain). We consider the unipartite projection of the network on the top domain. Basic properties of the metabolic network are summarized in Table 3.

A.2 Internet

The Internet network is a snapshot of the autonomous system level Internet taken from the University of Oregon Route Views Project rou . The full dataset contains $733$ daily instances which span an interval of $785$ days from November 8 1997 to January 2 2000. Here we use a network instance as of January 2, 2000 met .

A.3 PGP web of trust

PGP is a data encryption and decryption computer program that provides cryptographic privacy and authentication for data communication ope . The data is collected and maintained by Cederlöf ced . In the paper we use the PGP snapshot taken in April of 2003. The PGP web of trust is a directed network where nodes are certificates consisting of public PGP keys and owner information. A directed link in the web of trust pointing from certificate A to certificate B represents a digital signature by the owner of A endorsing the owner/public key association of B. We construct the undirected PGP graph by taking into account only bi-directional trust links between the certificates. Further, we only consider the giant connected component of the resulting undirected PGP web of trust network. Basic properties of the PGP network are summarized in Table 3.

Appendix B Link Prediction: Alternative Methods and Scoring Techniques

We compare the accuracy of the hyperlink link prediction method against the following set of link prediction methods: CN Liben-Nowell and Kleinberg (2003), AA Adamic and Adar (2003), RA Zhou et al. (2009)), CRA Cannistraci et al. (2013), JC Jaccard (1901), SPM Lü et al. (2015), and SBM Guimerà and Sales-Pardo (2009); Peixoto (2018) methods.

All these methods, as well as the hyperlink, assign scores to (a subset of) all not directly connected pairs of nodes (nonlinks), and all such pairs are then ranked according to these scores from the most to least likely interaction prediction. To briefly describe these methods, it is thus sufficient to tell how these scores are calculated, for which we use the following notations: $k_{i}$ is the degree of node $i$; $\Gamma(i)$ is the set of $i$’s neighbors (directly connected nodes); $\gamma_{ij}(s)$ is the subset of all $\Gamma(s)$ that are neighbors of both $i$ and $j$; $e_{i}^{j}$ is $i$’s $j$-external degree, the number of $i$’s neighbors that are not $j$’s neighbors; $\boldsymbol{A}$ is the network adjacency matrix.

B.1 Common Neighbors (CN)

The score for a pair of nodes $i$ and $j$ is defined as the cardinality of the intersection of their sets of neighbors,

[TABLE]

B.2 Jaccard’s index (JC)

The score is a normalized measure of the overlap of $i$’s and $j$’s sets of neighbors,

[TABLE]

B.3 Adamic-Adar index (AA)

The score assigns more weight to the less-connected neighbors,

[TABLE]

B.4 Resource Allocation index (RA)

The score is similar to the AA score, but punishes high-degree nodes more strongly,

[TABLE]

B.5 Cannistraci Resource Allocation index (CRA)

The score is similar to the RA score, but takes into account the subset of nodes shared between nodes $i$, $j$ and their common neighbors $s$:

[TABLE]

B.6 Structural Perturbation Method (SPM)

This method is based on repetitive perturbations of the adjacency matrix $\boldsymbol{A}$ by removals of small fractions of links that we denote by $\Delta E$. The original adjacency matrix can then be written as $\boldsymbol{A}=\boldsymbol{A^{\prime}}+\boldsymbol{\Delta A}$, where $\boldsymbol{A^{\prime}}$ is the adjacency matrix of the network after removal of links $\Delta E$, and $\boldsymbol{\Delta A}$ is the adjacency matrix constructed on the set of removed links $\Delta E$. Denoting eigenvectors and eigenvalues of $\boldsymbol{A^{\prime}}$ by $x_{k}$ and $\lambda_{k}$, the perturbations of the original eigenvalues $\lambda_{k}$ using the perturbation matrix $\boldsymbol{\Delta A}$, are

[TABLE]

so that the perturbed adjacency matrix is

[TABLE]

All nonlinks $i,j$ are then ranked by $\widetilde{A}_{ij}$. In our experiments we repeat this perturbation procedure ten times, and then average perturbed matrices over these trials, thus obtaining an averaged perturbed matrix $\langle\boldsymbol{\widetilde{A}}\rangle$, so that the SPM score is

[TABLE]

B.7 Stochastic Block Model (SBM)

The stochastic block model is a generative network model designed to model community structure. Nodes are partitioned into groups (blocks) forming a node partition $\boldsymbol{b}$. The number of links between blocks is given by a matrix $\boldsymbol{e}$ the elements $e_{rs}$ of which are the numbers of links between blocks $r$ and $s$. If the observed degree sequence of a network $\boldsymbol{k}$ is used as an additional model parameter, the model is called degree-corrected SBM Peixoto (2017). Moreover, if node blocks are themselves clustered into groups, and these groups are organized into higher-level groups, and so on recursively up to some nestedness level $l$, the model is called nested SBM. It can capture hierarchical and fine-grained structural properties of a given network Peixoto (2014a). In our experiments, we use both degree-corrected and nested SBMs denoted as SBM(d) and SBM(n) in the figures. We rely on the graph-tool library Peixoto (2014b) in the procedures below. Given the observed data (network adjacency matrix) $\boldsymbol{\mathcal{D}}$, and the prior probability density $P(\boldsymbol{A},\boldsymbol{b})$ given by the network block structure $\boldsymbol{b}$ that produces a network with adjacency matrix $\boldsymbol{A}$ in the model, we reconstruct the full network using the posterior distribution:

[TABLE]

where $P(\boldsymbol{\mathcal{D}}|\boldsymbol{A})$ describes the measurement process of a network. We avoid computing the normalization factor $P(\boldsymbol{\mathcal{D}})$ by Markov Chain Monte Carlo (MCMC) sampling from the joint posterior distribution $P(\boldsymbol{A},\boldsymbol{b}|\boldsymbol{\mathcal{D}})$ as described in Peixoto (2018). In our experiments, we assume that each network link is observed and measured once. A possible link $(i,j)$ then has marginal probability

[TABLE]

To sample over $(\boldsymbol{A},\boldsymbol{b})$ configurations, the MCMC algorithm is initialized with the block structure obtained by the procedure from Peixoto (2014c). The MCMC is equilibrated using $10|E|$ equilibration steps, where $E$ is the set of links in a given network. For the PGP network, due to its large size, we use only $2|E|$ MCMC equilibration steps. Then $T=10$ epochs of MCMC iterations, $1,000$ swaps each, are performed to sample different block-network configurations. After each epoch, marginal link probabilities from Eq. (43) are collected. These probabilities are then averaged over the epochs to obtain a single score used for link prediction:

[TABLE]

Appendix C Basic properties of the RHG

The hyperbolic geometry inference algorithm relies on several properties of the RHG, which we review in this section.

Degree distribution. RHGs are characterized by scale-free degree distributions, $P(k)\sim k^{-\gamma}$, where $\gamma=2\alpha+1$. Indeed, the expected degree of a node located at $(r,\theta)$ is independent of its angular coordinate $\theta$, $\overline{k}(r,\theta)=\overline{k}(r,0)=\overline{k}(r)$, and is given by

[TABLE]

see Krioukov et al. (2010). The average degree of the model is given by

[TABLE]

As seen from Eq. (46), $\overline{k}$ in the most general case depends on the network size $N$.

To achieve sparse models with $\overline{k}$ independent of $N$ one sets the radius of the hyperbolic disk to

[TABLE]

where $\nu>0$ is the tuning parameter, directly related to $\overline{k}$. Indeed, with $R(N)$ given by (47)

[TABLE]

prescribing the value of $\nu$ for the target values of $\overline{k}$, $\alpha$, and $T$.

It has been shown in Boguñá and Pastor-Satorras (2003) that in the sparse limit the probability of a node located at $(r,\theta)$ to have $k$ connections can be approximated with the Poisson distribution with the mean of $\overline{k}(r)$:

[TABLE]

Then the degree distribution of the RHG is

[TABLE]

It follows from Eqs. (45) and (48) that model parameters $\alpha$ and $R$ can be used to control degree distribution exponent $\gamma$ and the average degree of the model, respectively.

Clustering coefficient. As seen from Eq. (6), connection probability $p(x)$ decreases exponentially for distances $x>R$ with the rate of $\frac{1}{2}T$. Thus, the temperature parameter $T$ tunes the role of large distances in the formation of links: the higher the $T$ the more likely are long-distance connections. As a result, $T$ controls the clustering coefficient of the RHG. In the $T\to 0$ limit connections are only possible at hyperbolic distances $x<R$ and the clustering coefficient is maximized. Conversely, the clustering coefficient decreases as $T$ increases and vanishes asymptotically in the $T\geq 1$ case Krioukov et al. (2010).

Appendix D HYPERLINK Accuracy

In this section we calculate analytically the hyperlink accuracy, in terms of AUC and AUPR, on RHGs with known coordinates. Our results in this section are confirmed by the numerical experiments in Sec. III and build our intuition for Sec. IV, where we analyze hyperlink on RHGs and real networks with inferred coordinates.

D.1 AUC

To understand the behavior of AUC scores as a function of RHG parameters we define distance-dependent true positive ${\rm tpr}(x)$ and false positive ${\rm fpr}(x)$ rates as the fractions of true and false positives, respectively, contained among unconnected node pairs separated by distances up to $x$:

[TABLE]

where $E$ is the true number of links in the network, $E=|\Omega_{E}\cup\Omega_{R}|$, $p(y)$ is the connection probability in the RHG given by Eq. (6) and $n(y)$ is the distance distribution for node pairs in the RHG, given by Eq. (27).

It is seen from Eqs. (52) and (53) that in the $T\to 0$ limit $p(y)=\Theta(R-y)$, resulting in ${\rm fpr}(x)=0$ for $x\leq R$ and ${\rm tpr}(x)=1$ for $x\geq R$, resulting in the ideal ROC curve, Fig. 2a, and ${\rm AUC}=1$.

Using the expression for $n(y)$ from Eq. (27), we can evaluate true and false positive rates, up to the proportionality coefficient, as:

[TABLE]

where

[TABLE]

and$~{}_{2}F_{1}$ is the Gaussian hypergeometric function. In the $z\ll 1$ regime $I\left(z;T\right)\approx z$ and, thus, ${\rm tpr}(x)\sim e^{\frac{x-R}{2}}$ and ${\rm fpr}(x)\approx 0$ for $x<R$, Fig. 4a,b.

In the $z\gg 1$ regime $I\left(z;T\right)\approx I(T)$, where $I(T)=\frac{\pi}{T\sin\left(\pi/T\right)}$, explaining the saturation of the true positive rate, ${\rm tpr}(x)\to 1$ as $x$ approaches $2R$, and the exponential growth of the false positive rate, ${\rm fpr}(x)\sim e^{\frac{x}{2}}$ for $x>R$, Fig. 4a,b.

To obtain the analytical estimate of the ${\rm AUC}$ as a function of RHG parameters we represent it as

[TABLE]

By making use of Eqs. (52) and (53) we arrive at

[TABLE]

where $n^{c}(x)\equiv\int_{0}^{x}n(y){\rm d}y$.

In the case of sparse networks the first correction term $\Delta_{1}\sim N^{-1}$ and can be ignored in the large $N$ limit. The second correction term requires further analysis. It is straightforward to verify that in the $T\to 0$ limit $\Delta_{2}\sim N^{-1}$ and can also be ignored. Indeed, in this case $p^{\prime}(x)=-\delta(x-R)$, and

[TABLE]

Since $\binom{N}{2}\left[n^{c}(R)\right]$ equals the number of node pairs in the hyperbolic disk with distances up to $R$ and all these node pairs are connected in the $T\to 0$ case, $\binom{N}{2}\left[n^{c}(R)\right]=E$, resulting in $\Delta_{2}(T=0)=\frac{E}{\binom{N}{2}}\sim N^{-1}$.

To estimate the behavior of $\Delta_{2}$ in the case of $T>0$ we need to understand the behavior of its integrand in Eq. (60). Since $n^{c}(x)\sim e^{\frac{x}{2}}$ and $-p^{\prime}(x)=\frac{1}{2T}\exp\left(\frac{x-R}{2T}\right)\left[p(x)\right]^{2}$, the integrand is sharply peaked at $x=R+2T\ln\left(\frac{1+T}{1-T}\right)$ in the case of $T\in\left(0,\frac{1}{2}\right)$, resulting in $\Delta_{2}\sim N^{-1}$, similar to the $T\to 0$ case.

Conversely, the integrand in Eq. (60) grows monotonously as a function of $x$ in the case of $T\in\left(\frac{1}{2},1\right)$. The evaluation of $\Delta_{2}$ in this regime is quite involved and is not informative. Instead, we elect to compute the upper bound for $\Delta_{2}$, which also provides the lower bound for AUC scores. In doing so we note that the leading term behavior of $n(x)$ given by Eq. (27) is also its upper bound, see Ref. Alanis-Lobato and Andrade-Navarro (2016). Then

[TABLE]

since $e^{\frac{R}{2}}\sim N$ in the case of sparse RHGs, see Eq. (47). In the case of $T=\frac{1}{2}$ Eq. (62) simplifies to

[TABLE]

Taken together, the results above show that the AUC scores for RHGs with known coordinates converge to $1$ in the large $N$ limit as

[TABLE]

D.2 AUPR

AUPR scores can be evaluated in a similar fashion:

[TABLE]

where ${\rm pr}(x)$ and ${\rm rc}(x)$ are, respectively, distance-dependent precision and recall functions for hyperbolic distances up to $x$:

[TABLE]

where $N_{d}(x)$ is the number of disconnected node pairs with distances up to $x$:

[TABLE]

Using Eqs. (52) and (68) we obtain

[TABLE]

In the $T\to 0$ limit ${\rm pr}(x)=1$ for all $x<R$, while ${\rm rc}(x)=E(x)/E$, resulting, as expected, in ${\rm AUPR}=1$. Here $E(x)$ is the cumulative number of links between the node pairs with distances up to $x$.

In the $T>0$ case we rely on Eqs.  (52), (27), and (68) to obtain

[TABLE]

where $I\left(z;T\right)$ is given by Eq. (56).

In the $x\ll R$ regime $I\left(e^{\frac{x-R}{2}};T\right)\sim e^{\frac{x-R}{2}}$ and ${\rm pr}(x)\to 1$. In the $x\gg R$ case $I\left(e^{\frac{x-R}{2}};T\right)\sim\frac{\pi}{T\sin\left(\pi/T\right)}$, and as a result, precision decays exponentially, ${\rm pr}(x)\sim e^{-x/2}$, independent of $T$, Fig. 4c.

The dependence of ${\rm AUPR}$ on $T$ arises from the recall function or its derivative, $rc^{\prime}(x)$, quantifying the expected distance-dependent link density and, consequently, the density of missing links.

[TABLE]

$E_{c}(x)$ grows exponentially as $e^{x/2}$ for $x\ll R$ values and decays as $e^{x\left(1-\frac{1}{T}\right)}$ for $x\gg R$, reaching the maximum at $x^{*}=R-2T\ln\left(\frac{1}{T}-1\right)$, Fig. 4d. Thus, as $T$ increases, the missing links are more likely to be located at larger distances where precision ${\rm pr}(x)$ is smaller, resulting in lower AUPR scores, consistent with the observations in Fig. 3.

Appendix E Effects of Coordinate Uncertainty on HYPERLINK Accuracy

To understand the effects of coordinate uncertainties on hyperlink accuracy we model coordinate inference uncertainty as synthetic noise that we add to true angular coordinates of the RHG. In the following we first generate RHG as described in Sec. II.2 and then simulate uncertainties of angular coordinates by adding synthetic noise to original angular coordinates:

[TABLE]

where $a>0$ is the noise amplitude. Further, we conduct link prediction experiments by calculating latent distances with uncertain coordinates:

[TABLE]

where $x$ is calculated according to the hyperbolic law of cosines, Eq. (3).

E.1 Link prediction with noise

In the case of synthetic noise, the AUPR scores are still given by Eq. (29) with effective precision and recall rates ${\rm pr}(y|a)$ and ${\rm rc}(y|a)$ evaluated in the presence of noise. To calculate these rates we start with the effective true positive rate ${\rm tp}(y|a)$.

To this end, we first define the subgraph $G_{y}$ obtained from the RHG $G$ by keeping only links between node pairs separated by distances at most $y$. Then, it is easy to realize that the true positive rate ${\rm tp}(y)$ is proportional to the expected degree $\overline{k}_{y}$ of the $G_{y}$:

[TABLE]

$\overline{k}_{y}$ can be calculated using the hidden variable formalism

[TABLE]

To account for noise we next define noisy subgraph $G_{y}(a)$ as follows. First, noise is added to node coordinates of the original RHG as prescribed by Eq. (30) and hyperbolic distances between nodes are recalculated using the updated coordinates. Second, $G_{y}(a)$ is formed from RHG by keeping connections at recalculated distances up to $y$. It is then easy to see that the thought true positive rate is given by

[TABLE]

where $\overline{k}_{y}(a)$ is the average degree of noisy subgraph $G_{y}(a)$.

After a series of tedious calculations, which we detail in the Subsection E.2, we obtain the leading order behavior of $\overline{k}_{y}(a)$:

[TABLE]

where $\alpha\in\left(\frac{1}{2},1\right)$ is the radial node density parameter in Eq. (5) corresponding to degree distribution exponent $\gamma=2\alpha+1$. Similar to the noiseless case, $g(y)$ grows as $\exp\left(\frac{y}{2}\right)$ for $y\leq R$ and saturates to a constant value, corresponding to $\overline{k}_{y}(a)=\overline{k}$ as $y\to 2R$, Fig. 14a.

Using Eq.(69) one can rewrite the distance-dependent precision function as

[TABLE]

where $n(y|a)$ is the node pair distribution in the hyperbolic disk with coordinate noise.

Due to the uniform initial angular distribution $\rho(\theta)$, the node pair distribution is independent of noise, $n(y|a)=n(y)$, Fig. 14b. Further, in the case of sufficiently large noise amplitude $a$, ${\rm tp}(y|a)\ll\binom{N}{2}\int_{0}^{y}n(y^{\prime}|a){\rm d}y^{\prime}$ and

[TABLE]

As a result, in the case $y\leq R$, ${\rm pr}(y|a)\sim a^{1-2\alpha}$, see Fig. 14c.

Since the distance-dependent recall function is proportional to the true positive rate,

[TABLE]

The resulting AUPR score scales as

[TABLE]

where

[TABLE]

see Fig. 14d.

This result suggests that the impact of coordinate uncertainty on link prediction is higher in RHG with larger $\gamma=2\alpha+1$ values. Intuitively, this is the case since networks with larger $\gamma$ values have larger fractions of small degree nodes. Small degree nodes in the RHG are characterized by large radial coordinates, and the hyperbolic distance between the points with large radial coordinates is most affected by angular coordinate uncertainties.

E.2 The average degree of the noisy subgraph

Here we derive the leading term behavior of the average degree of the noisy subgraph $G_{y}(a)$ as a function of noise amplitude $a$

As shown in the subsection above, the number of true positives ${\rm tp}(y|a)$ is related to the average degree of noisy subgraph $G_{y}(a)$. To define $G_{y}(a)$ we add uniform noise of amplitude $a$ to original angular coordinates of the RHG and calculate noisy hyperbolic distances $\hat{x}_{ij}$ between all node pairs using noisy coordinates. $G_{y}(a)$ is the RHG subgraph formed by node pairs with noisy hyperbolic distances $\hat{x}_{ij}<y$. The average degree of $G_{y}(a)$ is given by

[TABLE]

Here $\rho(r)$ is given by Eq. (5), and $\rho(\theta|\hat{\theta})$ is the conditional probability of the true angle $\theta$, given inferred angle $\hat{\theta}$. In case of the uniform noise, $\rho(\theta|\hat{\theta})$ is also a uniform distribution centered at $\hat{\theta}$:

[TABLE]

while

[TABLE]

Throughout the calculation of $k_{y}(a)$ we will rely on the number of assumptions. We are primarily interested in RHGs with $2<\gamma<3$, which correspond to $\frac{1}{2}<\alpha<1$. To identify leading terms we will also recall on the scaling of $R$ with the system size, $N\sim e^{\frac{R}{2}}$.

Since hyperbolic distance $x$ in Eq. (3) depends on $\theta_{1}$ and $\theta_{2}$ only through their difference,

[TABLE]

and angles distributed uniformly on $[-\pi,\pi]$, $\rho(\hat{\theta}_{1,2})=\frac{1}{2\pi}$, we can simplify Eq. (87) as

[TABLE]

where

[TABLE]

and $\Theta[x]$ is the Heaviside theta function. Similar to the calculation of $\overline{k}$ in the RHGs, Ref. Krioukov et al. (2010), we can rewrite Eq. (92) as

[TABLE]

where $\overline{k}_{y}(r|a)$ is the average degree of node with radial coordinate $r$ in noisy subgraph $G_{y}(a)$:

[TABLE]

and angles $\phi\equiv\Delta\theta_{12}$ and $\hat{\phi}\equiv\Delta\hat{\theta}_{12}$ are introduced to ease the notation.

To evaluate $\overline{k}_{y}(r_{1}|a)$ we note that the integration region in Eq. (95) is given by intersection of two hyperbolic disks. The first one is of radius $R$ and is centered at the coordinate system origin, $(0,0)$. The second disk is of radius $y$ and is centered at $\left(r_{1},0\right)$.

We perform the integration for the two regimes of $y\in\left[\frac{R}{2},R\right]$ and $[R,2R]$ separately. We do not perform the integration for the $y\in\left[0,\frac{R}{2}\right]$ regime since the number of true positives here is much smaller than that in the other two regimes. This is the case since $n(y)$ grows exponentially with $y$, $n(y)\sim e^{\frac{y}{2}}$. Consequently, the number of possible true positives in the $y\in\left[0,\frac{R}{2}\right]$ regime is much smaller than that in the $y\in\left[\frac{R}{2},R\right]$ regime.

E.2.1 $y\in\left[\frac{R}{2},R\right]$

To evaluate $\overline{k}_{y}(r_{1}|a)$ we perform the integration over $r_{1}$ and $r_{2}$ values over the domain shown in Fig. 15d. Based on this domain, it is convenient to split the integration over $r_{1}$ into three regions, $0\leq r_{1}\leq R-y$, $R-y\leq r_{1}\leq y$, and $y\leq r_{1}\leq R$. However, due to specifics of the approximation techniques, it is more convenient to split the integration not into three but into five regions — (i) $0\leq r_{1}\leq\frac{R-y}{2}-\ln\frac{a}{2}$, (ii) $\frac{R-y}{2}-\ln\frac{a}{2}\leq r_{1}\leq R-y$, (iii) $R-y\leq r_{1}\leq y$, (iv) $y\leq r_{1}\leq\frac{R+y}{2}-\ln\frac{a}{2}$, and (v) $\frac{R+y}{2}-\ln\frac{a}{2}\leq R$ — which we depict for convenience in Fig. 15d with vertical dashed lines. We evaluate the contributions to $\overline{k}_{y}(r_{1}|a)$ from each of these five regions below.

Region I: $0\leq r_{1}\leq\frac{R-y}{2}-\ln\frac{a}{2}$. In this region the disk $y$ is fully contained within the disk $R$. Further, since $y>R/2$, disk $y$ is guaranteed to include the coordinate system origin for all $r_{1}\in[0,R-y]$ values, Fig. 15a. In this case the integral in $\overline{k}_{y}(r_{1}|a)$ can be evaluated as

[TABLE]

where $\phi_{y}$ is the angle given by the intersection of the disk with radius $r_{2}$ centered at $r=0$ and that of radius $y$, centered at $r=r_{1}$. To estimate $\phi_{y}$ we consider the triangle formed by the origin $(0,0)$, disk $y$ centered at $(r_{1},0)$, and the intersection of $r_{2}$ with $y$. The triangle has sides equal to $r_{1}$, $r_{2}$, and $y$ with $\phi_{y}$ being the angle between $r_{1}$ and $r_{2}$. Thus, $\phi_{y}$ is given by the hyperbolic law of cosines:

[TABLE]

In the case of sufficiently large $r_{1}$, $r_{2}$, and $y$ values we can approximate $\cos\phi_{y}$ as

[TABLE]

Since $\hat{\phi}$ in the first integral sweeps the entire $2\pi$ angle, $\mathfrak{I}_{1}$ is given by

[TABLE]

Then, since $x\left(r_{1},r_{2},\hat{\phi}\right)\leq r_{1}+r_{2}\leq R$, $p\left[x\left(r_{1},r_{2},\hat{\phi}\right)\right]\approx 1$, leading to

[TABLE]

The evaluation of $\mathfrak{I}_{2}$ is more involved and requires further approximations. We notice that $\phi_{y}\ll 1$ since $r_{2}\in[y-r_{1},y+r_{1}]$, which can be further approximated as

[TABLE]

Then, for sufficiently large noise amplitudes $a\gg\phi_{y}$, we can approximate the integral $\int_{\hat{\phi}-a}^{\hat{\phi}+a}{\rm d}\phi$ as $2\int_{0}^{a}{\rm d}\phi$, resulting in

[TABLE]

Since $r_{1}<\frac{R-y}{2}-\ln{\frac{a}{2}}$, $r_{2}<y+r_{1}$, and $y<R$, it follows that $x\left(r_{1},r_{2},\phi\right)<r_{1}+r_{2}+2\ln\frac{a}{2}<R$, and, as a result, ${\rm exp}\left(\frac{x\left(r_{1},r_{2},\phi\right)-R}{2T}\right)\ll 1$, resulting in

[TABLE]

Since $\gamma>2$ case ($\alpha>\frac{1}{2}$), $\mathfrak{I}_{2}\gg\mathfrak{I}_{1}$, and

[TABLE]

Region II: $\frac{R-y}{2}-\ln\frac{a}{2}\leq r_{1}\leq R-y$.

Similar to region I, the hyperbolic disk $y$ fully lies within disk $R$, Fig. 15a. Thus, $\overline{k}_{y}(r_{1}|a)$ is given by the same expression:

[TABLE]

The calculation of $\mathfrak{I}_{3}$ is identical to that of $\mathfrak{I}_{1}$, resulting in

[TABLE]

Different from region I is the calculation of $\mathfrak{I}_{4}$. Indeed, in the case $r_{1}\geq\frac{R-y}{2}+\ln\frac{a}{2}$, and $r_{2}\in[y-r_{1},y+r_{1}]$ hyperbolic distance $x\left(r_{1},r_{2},\phi\right)$ is no longer guaranteed to be smaller than $R$, and $p\left[x\left(r_{1},r_{2},\phi\right)\right]$ can no longer be approximated by unity. We first split $\mathfrak{I}_{4}$ into two parts and calculate them separately:

[TABLE]

where

[TABLE]

By approximating the hyperbolic law of cosines in Eq. (3) as, $x\left(r_{1},r_{2},\phi\right)\approx r_{1}+r_{2}+2\ln\frac{\phi}{2}$ and making use of Eq. (103) we obtain for $\mathfrak{I}_{41}$

[TABLE]

where $I(z;T)\equiv\int_{0}^{z}\frac{{\rm d}x}{1+x^{\frac{1}{T}}}$ is the same function as in Eq. (56).

Recall that for small $z\ll 1$ function $I(z;T)\approx z$, while, for $z\gg 1$, $I(z;T)\approx I(T)=\frac{\pi}{T\sin\left(\frac{\pi}{T}\right)}$. With these approximations in mind we split the integration in $\mathfrak{I}_{41}$ into two subregions :

[TABLE]

In the first subregion, $r_{2}\in\left[y-r_{1},R-r_{1}-2\ln\frac{a}{2}\right]$, and $\frac{a}{2}e^{\frac{r_{1}+r_{2}-R}{2}}\leq 1$, which allows us to approximate $I\left(\frac{a}{2}e^{\frac{r_{1}+r_{2}-R}{2}};T\right)\approx\frac{a}{2}e^{\frac{r_{1}+r_{2}-R}{2}}$. In the second subregion, $r_{2}\in[R-r_{1}-2\ln\frac{a}{2},y+r_{1}]$, $\frac{a}{2}e^{\frac{r_{1}+r_{2}-R}{2}}\geq 1$, and $I\left(\frac{a}{2}e^{\frac{r_{1}+r_{2}-R}{2}};T\right)\approx I\left(T\right)$. Using these approximations we obtain, to the leading order,

[TABLE]

Following the same approximation steps,

[TABLE]

where

[TABLE]

in the case $T<1/2$.

Taken together, $\mathfrak{I}_{4,1}$ and $\mathfrak{I}_{4,2}$ result in

[TABLE]

Finally, since $y<R$, we conclude that $\mathfrak{I}_{3}\ll\mathfrak{I}_{4}$, resulting in

[TABLE]

for $\frac{R-y}{2}-\ln\frac{a}{2}\leq r_{1}\leq R-y$.

Region III: $R-y\leq r_{1}\leq y$.

In this region disk $y$ is partially contained within the disk $R$. Since $r_{1}\leq y$, disk $y$ still contains the coordinate system origin, Fig. 15b. Similar to regions I and II, we split the calculation of $\overline{k}_{y}(r_{1}|a)$ into two parts:

[TABLE]

where $\phi_{y}\ll 1$ is the intersection angle of disk $r_{2}$ with that of $y$, see Fig. 15b, and is given by Eq. (103).

We first note that the integral in $\mathfrak{I}_{5}$ is identical to those in $\mathfrak{I}_{3}$ and $\mathfrak{I}_{1}$:

[TABLE]

The integration in $\mathfrak{I}_{6}$ is very similar to that in $\mathfrak{I}_{4}$ with the only difference in the upper integration bound of $r_{2}\leq R$. The evaluation of $\mathfrak{I}_{6}$ is, therefore, straightforward and requires the same approximation steps as in $\mathfrak{I}_{4}$. A quicker estimate can be obtained by noting that the upper bound for $r_{2}$ in $\mathfrak{I}_{4}$ does not contribute to the leading term. The reason is that $\mathfrak{I}_{42}$ is dominated by $r_{2}$ in the vicinity of the $r_{2}=R-r_{1}-2\ln\frac{a}{2}$ point.

Since $R>R-r_{1}-2\ln\frac{a}{2}>y-r_{1}$

[TABLE]

with integrands identical to those of $\mathfrak{I}_{41}$ and $\mathfrak{I}_{42}$. Since the integrand in $\mathfrak{I}_{42}$ is dominated by smaller $r_{2}$ values we conclude that

[TABLE]

Finally, $\mathfrak{I}_{6}$ dominates $\mathfrak{I}_{5}$ for $\alpha>\frac{1}{2}$, resulting in

[TABLE]

for $R-y\leq r_{1}\leq y$.

Region IV: $y\leq r_{1}\leq\frac{R+y}{2}-\ln\frac{a}{2}$.

In this region, hyperbolic disk $y$ is partially contained within $R$ and does not include the origin, Fig. 15c. Therefore, in this region

[TABLE]

Using the arguments similar to that of region III, we obtain

[TABLE]

for $y\leq r_{1}\leq\frac{R+y}{2}-\ln\frac{a}{2}$.

Region V: $\frac{R+y}{2}-\ln\frac{a}{2}\leq r_{1}\leq R$.

Similar to the situation in region IV, hyperbolic disk $y$ intersects disk $R$ and does not include the coordinate system origin. Different from region IV is the $r_{2}=R-r_{1}-2\ln\frac{a}{2}$ point that lies outside the $r_{2}$ integration region and we can no longer relate $\overline{k}_{y}(r_{1}|a)$ to those in other regions.

To evaluate

[TABLE]

we recall that $\phi_{y}\ll 1$, and for sufficiently large $a\gg\phi_{y}$ we obtain

[TABLE]

After straightforward approximations we obtain

[TABLE]

for $\frac{R+y}{2}-\ln\frac{a}{2}\leq r_{1}\leq R$

Merged together, Eqs. (106),(120), (127), (129), and (134) provide the solution for $\overline{k}_{y}(r_{1}|a)$:

[TABLE]

Using Eq. (135) together with Eq. (94) we finally obtain

[TABLE]

E.2.2 $y\in[R,2R]$

In the regime $y\geq R$ hyperbolic disk $y$ always contains the origin, Fig. 16. To evaluate $\overline{k}_{y}(r_{1}|a)$ in this regime we need to distinguish two cases, (VI) $0\leq r_{1}\leq y-R$, and (VII) $y-R\leq r_{1}\leq R$.

Region VI: $0\leq r_{1}\leq y-R$.

In this regime hyperbolic disk $R$ is fully contained within hyperbolic disk $y$, Fig. 16a, and $\overline{k}_{y}(r_{1}|a)=\overline{k}(r_{1})$, where $\overline{k}(r_{1})$ is the average degree of a node at $r_{1}$ in the RHG. Indeed, radial coordinates of all points are within disk $R$, and all distances from point $(r_{1},0)$ to any point within disk $R$ are guaranteed to be smaller than $y$, $x(r_{1},0,r_{2},\theta)<y$ for any $\theta\in[0,2\pi]$. Therefore in this regime

[TABLE]

Since the integral over $\hat{\phi}$ sweeps the entire circle, $\hat{\theta}\in[0,2\pi]$, synthetic noise does not affect the integration:

[TABLE]

resulting in

[TABLE]

in the case $0\leq r_{1}\leq y-R$.

Region VII: $R-y\leq r_{1}\leq R$. In this regime hyperbolic disk $R$ is partially contained within $y$ and the calculation of $\overline{k}_{y}(r_{1}|a)$ splits into two integrals:

[TABLE]

where $\phi_{y}$ is the angle of intersection of disks $R$ and $y$, Fig. 16b.

We note that the integration region for $\mathfrak{I}_{9}$ is identical to that of $\mathfrak{I}_{1}$. Different from the case of $\mathfrak{I}_{1}$ is the condition that $y>R$. In this case $x\left(r_{1},r_{2},\phi\right)$ is no longer guaranteed to be less than $R$, and $p\left[x\left(r_{1},r_{2},\phi\right)\right]$ cannot be approximated by $1$. We start evaluating $\mathfrak{I}_{9}$ by performing the integration over $\phi$, which leads to

[TABLE]

where $I(z;T)$ is given by Eq. (56). Recall that $I(z;T)\approx z$ if $z\ll 1$ and $I(z;T)\approx I(T)$ in case $x\gg 1$. Thus, to evaluate $\mathfrak{I}_{9}$ we split the integration over $r_{2}$ in to two integrals, $\int_{0}^{y-r_{1}}=\int_{0}^{R-r_{1}-2\ln\frac{\pi}{2}}+\int_{R-r_{1}-2\ln\frac{\pi}{2}}^{y-r_{1}}$. In the first integral $\frac{\pi}{2}e^{\frac{r_{1}+r_{2}-R}{2}}<1$ and we approximate $I\left(\frac{\pi}{2}e^{\frac{r_{1}+r_{2}-R}{2}};T\right)\approx\frac{\pi}{2}e^{\frac{r_{1}+r_{2}-R}{2}}$, while in the second integral $\frac{\pi}{2}e^{\frac{r_{1}+r_{2}-R}{2}}>1$ and $I\left(\frac{\pi}{2}e^{\frac{r_{1}+r_{2}-R}{2}};T\right)\approx I(T)$. The remaining integration steps in $\mathfrak{I}_{9}$ are straightforward, resulting in

[TABLE]

Finally, since $y>R$ and $\alpha>\frac{1}{2}$, we get

[TABLE]

In order to calculate $\mathfrak{I}_{10}$ we first need to estimate the cutoff angle $\phi_{y}$, which is given by the intersection of disks $R$ and $y$, and is given by Eq. (100). $\phi_{y}$ takes values from $\phi_{y}\approx 2e^{\frac{y-2R}{2}}$ at $r_{1}=r_{2}=R$ to $\phi_{y}=\pi$ at $r_{2}=y-r_{1}$. Thus, we can no longer use the $\phi_{y}\ll a$ approximation, as in $\mathfrak{I}_{2}$.

To proceed further we note that the integration domain in $\mathfrak{I}_{10}$ is given by the area above the $r_{2}=y-r_{1}$ line, Fig. 16c. We recall that the integration in the case $y<R$ is dominated by points in the vicinity of the $r_{1}+r_{2}=R-2\ln\frac{a}{2}$ line [see red dashed line in Fig. 15c]. Let us assume that this is also the case in the $y\geq R$ regime [see red dashed line in Fig. 16c]. We next note that in the vicinity of the $r_{1}+r_{2}=R-2\ln\frac{a}{2}$ line $\cos\phi_{y}\approx 1-2e^{y-R-2\ln\frac{a}{2}}$. For sufficiently small noise amplitude, such that $y<R-2\ln\frac{a}{2}$, the cutoff angle $\phi_{y}\ll 1$ and can be approximated by Eq. (103), and we can employ the same approximation techniques as in $\mathfrak{I}_{2}$.

Our strategy now is to split the integration domain of $\mathfrak{I}_{10}$ into two parts by the curve $r_{2}=\tilde{R}(r_{1})$ such that (i) this curve is below the $r_{1}+r_{2}=R-2\ln\frac{a}{2}$ line, and (ii) above this curve, $r_{2}>\tilde{R}(r_{1})$, the cutoff angle $\phi_{y}\ll 1$. One possibility for such a curve is the $\tilde{R}(r_{1})=A-r_{1}$ line, where $A=\frac{y+R}{2}-\ln\frac{a}{2}$ [see green dashed curve in Fig. 16c].

Then region VII splits into two subregions, VIIA and VIIB, corresponding to $r_{1}\in\left[y-R,2\ln\frac{2}{a}\right]$ and $r_{1}\in\left[2\ln\frac{2}{a},R\right]$, respectively, see Fig. 16c. We expect that the contribution to $k_{y}(a)$ from VIIA to be much smaller than that from VIIB since the latter contains the $r_{1}+r_{2}=R-2\ln\frac{a}{2}$ line and the former does not. Therefore, we will estimate the upper bound for $k_{y}(r_{1}|a)$ in VIIA by replacing $\phi_{y}$ with $\pi$. In subregion VIIB we split the integration over $r_{2}$ into two intervals, $r_{2}\in\left[0,\tilde{R}(r_{1})\right]$ and $r_{2}\in\left[\tilde{R}(r_{1}),R\right]$.

Subregion VIIA: $y-R\leq r_{1}\leq 2\ln\frac{2}{a}$.

Here the integral splits into

[TABLE]

Following our strategy, we evaluate the upper bound for $\mathfrak{I}_{12}$ by replacing the integration limit of ${\phi_{y}}$ with $\pi$:

[TABLE]

Then,

[TABLE]

for $y-R\leq r_{1}\leq 2\ln\frac{2}{a}$.

Subregion VIIB: $2\ln\frac{2}{a}\leq r_{1}\leq R$.

Here we distinguish three intervals:

[TABLE]

where $\tilde{R}(r_{1})=\frac{y+R}{2}-\ln\frac{a}{2}-r_{1}$.

We evaluate the upper bound for $\mathfrak{I}_{14}$ by replacing the $\phi_{y}$ cutoff with $\pi$:

[TABLE]

leading to

[TABLE]

After the same calculation steps as in $\mathfrak{I}_{9}$ we obtain

[TABLE]

To evaluate $\mathfrak{I}_{15}$ we use the $\phi_{y}\ll 1$ assumption, which enables us to use Eq. (103). This approximation holds since $r_{2}>\tilde{R}(r_{1})$. Then, by following the same simplification steps as in $\mathfrak{I}_{4}$ we obtain

[TABLE]

Following the same evaluation steps as in $\mathfrak{I}_{4}$ we confirm that both $\mathfrak{I}_{151}$ and $\mathfrak{I}_{152}$ are dominated by points in the vicinity of $r_{1}+r_{2}=R-2\ln\frac{a}{2}$, resulting in

[TABLE]

By comparing Eqs. (161) and Eqs. (157) we establish that $\mathfrak{I}_{15}\gg\mathfrak{I}_{13}+\mathfrak{I}_{14}$ since $\tilde{R(r_{1})}<R$ and $y>R$, confirming our hypothesis and resulting in

[TABLE]

in case $2\ln\frac{2}{a}\leq r_{1}\leq R$.

Taken together, our results for regions VI and VII read

[TABLE]

Using Eq. (163) together with Eq. (94) we finally obtain

[TABLE]

Finally, we conclude that $k^{2}_{y}(a)\gg k^{1}_{y}(a)$ since $y>R$, which allows us to establish

[TABLE]

for $y>R$. Equation (136) together with Eq. (167) establish the baseline for calculation of ${\rm AUPR}(a)$ in Subsection E.1.

Appendix F HYPERLINK Embedder

The original hyperbolic geometry inference algorithm was developed in Ref. Boguñá et al. (2010) and is based on MLE. While the algorithm is rather slow with the overall computational complexity of $\mathcal{O}\left(N^{3}\right)$, it has been shown to accurately infer node coordinates in $\mathbb{H}^{2}$ leading to a number of promising applications ranging from interdomain Internet routing Boguñá et al. (2010) to understanding the growth of large-scale networks Papadopoulos et al. (2012).

In recent years hyperbolic geometry inference has become an active area of research and a collection of alternative inference methods has been developed by different research teams based on the MLE Bläsius et al. (2016); Wang et al. (2016), Laplacian eigenmaps Alanis-Lobato et al. (2016a, b); Muscoloni et al. (2017) and isomap Muscoloni et al. (2017). Even though most of these methods are characterized by relatively small computational complexity, $\mathcal{O}\left(N\right)$ - $\mathcal{O}\left(N^{2}\right)$, their inference accuracy has not been well explored.

At the same time, our initial experiments indicate that even small node coordinate uncertainties drastically reduce link prediction accuracy (Fig. 5). Therefore, to optimize link prediction results one needs to maximize the accuracy of node coordinate inference. To this end, we developed an enhanced MLE-based geometry inference algorithm, which we outline below.

F.1 General MLE formulation of hyperbolic geometry inference

Given the real network of interest with randomly removed links, we aim to find the set of node coordinates $\{\mathbf{x}_{i}\}\equiv\{(r_{i},\theta_{i})\}$, $i=1,2,\ldots,N$, in the hyperbolic disk $\mathbb{H}^{2}$ maximizing the probability $\mathcal{L}\left(\{\mathbf{x}_{i}\}|a_{ij},\mathcal{P},q\right)$ that node coordinates take particular values in the case the network is generated as the RHG with a subsequent random link removal process. Here $a_{ij}$ is the network’s observed adjacency matrix, and $\mathcal{P}$ is the set of parameters of the RHG, $\mathcal{P}=\{\alpha,T,R\}$.

By the Bayes rule the thought probability is given by

[TABLE]

where $\mathcal{L}\left(a_{ij}|\{\mathbf{x}_{i}\},\mathcal{P},q\right)$ is the likelihood that network $a_{ij}$ is generated as RHG with subsequent random link removal, ${\rm Prob}(\mathbf{x}_{i})$ is the prior probability of node coordinates generated by the RHG, and $\mathcal{L}\left(a_{ij}|\mathcal{P},q\right)$ is the probability that the network has been generated as the RHG with random link removal.

In the following we assume the uniform prior probability

[TABLE]

where $\rho(r_{i})$ are given by Eq. (5). Since node pairs are connected independently, the likelihood is given by

[TABLE]

where $\tilde{p}\left(x_{ij}\right)$ is the effective connection probability in the RHG generation process with subsequent random link removal:

[TABLE]

where $p(x)$ is the RHG connection probability function prescribed by Eq. (6).

The MLE inference aims to find node coordinates $\hat{\mathbf{x}}_{i}$ maximizing the likelihood $\mathcal{L}\left(\{\mathbf{x}_{i}\}|a_{ij},\mathcal{P},q\right)$, or equivalently, its logarithm

[TABLE]

where constant $K$ absorbs all terms independent of $\{\mathbf{x}_{i}\}$.

Our hyperbolic geometry inference procedure consists of three components: (1) finite-size effects and model parameter inference, (2) MLE-based inference of radial node coordinates, (3) MLE-based inference of angular node coordinates.

F.2 Finite-size effects and model parameter inference

The RHG has four parameters: the number of nodes $N$, hyperbolic disk radius $R$, node density parameter $\alpha$ and temperature $T$.

To infer $\alpha$ we first estimate the degree distribution exponent $\gamma$ through the inspection of the network degree distribution $P(k)$. Node density $\alpha$ is related to $\gamma$ through Eq. (51):

[TABLE]

The estimation of $N$ and in $R$ is less straightforward due to finite-size effects. First, in a real network one normally can only observe nodes with nonzero degrees. In contrast, the RHG may generate nodes of zero degree, which are accounted for in the calculation of the network’s average degree, $\overline{k}$, Eq. (46).

Second, due to finite-size effects, there is a cutoff value for the smallest node radius, $R_{0}$, affecting $\langle e^{-r/2}\rangle$ and, as a result, the observable $\overline{k}(r)$ and $\overline{k}$, Eqs. (46) and (45). Specifically, with the radius cutoff $R_{0}$

[TABLE]

where $\lambda(\alpha,x)$ is the finite-size correction coefficient:

[TABLE]

In the thermodynamic limit $\lambda(\alpha,\left(R-R_{0}\right))\to 1$ as

[TABLE]

However, in networks with $\alpha$ close to $1/2$ ($\gamma$ close 2) the rate of $\lambda$ convergence is slow and one needs to account for nonzero $R_{0}$.

Third, one needs to account for missing links that affect all observable properties of the RHG. In the particular case links are missing uniformly with probability $1-q$, the connection probability function $p(x)$ gets attenuated by the factor of $q$, Eq. (171), affecting all observable network properties.

Taken together, zero degree nodes, minimum radius cutoff, and missing links affect observable network properties as follows:

[TABLE]

where $\tilde{k}_{\rm max}$ is the maximum degree observed in the network and $P(0)$ is the fraction of zero degree nodes in the network. The latter can be estimated by averaging the conditional degree distribution $P(k=0|r)$ in Eq. (49) over possible $r$ values:

[TABLE]

where $\Gamma\left[s,x\right]$ is the upper incomplete gamma function.

Equations. (177), (178), (179), (180) and (181) allow one to infer the RHG parameters $R_{0}$, $R$, $N$ as well as resulting $\overline{k}$, and $P(0)$ by measuring observables $\tilde{N}$, $\tilde{k}$, and $\tilde{k}_{\rm max}$.

The caveat here is that parameter estimation presumes the knowledge of the missing link probability $1-q$. While this information is available in our synthetic experiments, it may not be available in real networks. In case the fraction of missing links is small, one can assume that $q=1$. The most general case of substantially incomplete networks where $q\ll 1$ is beyond the scope of this paper and will be studied elsewhere.

Finally, the temperature parameter $T$ needs to be estimated numerically by finding the solution of

[TABLE]

where $c_{0}$ is the average clustering coefficient of the network of interest and $\overline{c}(T)$ is the average clustering coefficient of the RHG generated with temperature $T$. We utilize this approach to infer $T$ of real networks in Sec. IV.2, while in experiments with RHGs we use actual $T$ values.

F.3 MLE-based inference of radial node coordinates

To infer radial node coordinates we extremize the logarithm of the likelihood function,

[TABLE]

obtaining

[TABLE]

In the case of sufficiently large $r$ values ${\rm coth}\left(\alpha r_{\ell}\right)\approx 1$. Further, one can approximate $x_{\ell j}$ as

[TABLE]

resulting in $\frac{\partial x_{\ell j}}{\partial r_{\ell}}\approx 1$. Taken together, these approximations allow us to simplify Eq.(184) as

[TABLE]

for $1-q\ll 1$. Note that the first summation in Eq. (186) is the degree of node $\ell$, $\sum_{j}a_{\ell j}=k_{\ell}$, while the second summation is the expected degree of the node with $r_{\ell}$, $\tilde{k}\left(r_{\ell}\right)=q\sum_{j}p\left(x_{\ell,j}\right)$. As a result, the value of $\hat{r}_{\ell}$ extremizing the likelihood is given by

[TABLE]

where $\tilde{k}\left(r\right)$ is the observable expected degree of the node with radial coordinate $r$. Since the latter is given by

[TABLE]

one can estimate $\hat{r}_{\ell}$ as

[TABLE]

F.4 MLE inference of angular node coordinates

To infer angular node coordinates one needs to maximize the likelihood ${\rm ln}~{}\mathcal{L}\left(\{\mathbf{x}_{i}\}|a_{ij},\mathcal{P},q\right)$ in Eq. (172) with respect to angular coordinates $\{\theta_{i}\}$, given the MLE values for radial coordinates $\{\hat{r}_{i}\}$. Since the maximization of ${\rm ln}~{}\mathcal{L}\left(\{\mathbf{x}_{i}\}|a_{ij},\mathcal{P},q\right)$ with respect to $\{\theta_{i}\}$ cannot be performed analytically, we have to rely on numerical approximations. To this end, we developed an MLE-based algorithm optimized for the linked prediction problem.

Conceptually, our algorithm is similar to the one developed in Ref. Boguñá et al. (2010) but has several important differences.

Following the exposition of Ref. Boguñá et al. (2010), we make two observations based on the link independence in RHG. First, angular coordinates of any node subset $\mathbb{S}$ can be inferred independently (albeit, with lower accuracy) based only on the partial information contained in the graph $G_{\mathbb{S}}$ formed by these nodes. In other words, the inference of angular coordinates in $\mathbb{S}$ is possible by maximizing the $\mathbb{S}$-specific log likelihood:

[TABLE]

Second, any log likelihood $\mathcal{L}\left[G_{\mathbb{S}}\right]$ can be represented as a sum of local contributions $\mathcal{L}\left[G_{\mathbb{S}}\right]_{i}$:

[TABLE]

where

[TABLE]

Since the log-likelihood profile ${\rm ln}~{}\mathcal{L}\left(\{\mathbf{x}_{i}\}|a_{ij},\mathcal{P},q\right)$ is nonconvex with abundant local maxima, we do not intend to find its global maximum by optimizing all angles at once. Instead, we proceed in a nested fashion by organizing network nodes into logarithmically sized layers with nodes of larger degree belonging to inner layers. To this end, we define the set $\mathbb{C}$ of all nodes with degrees $k>1$. We then rank all nodes in $\mathbb{C}$ in the decreasing order of their degree value, and split the resulting node list into $m$ layers with logarithmically growing sizes $s_{i}$, $i=0,..,m-1$:

[TABLE]

where $N(k>1)$ is the number of nodes with degree $k>1$, and $s_{0}\ll N$. Unless otherwise noted, we set $s_{0}=20$. Finally, all $k=1$ nodes are assigned to the outer layer $s_{m}$.

Complementary to layers $\{s_{i}\}$, we also define self-enclosed cores $\{cr_{i}\}$, $i=0,..,m$, such that core $cr_{i}$ contains all layer with indices $j\leq i$, $cr_{i}=\prod_{j=0}^{i}\bigcup s_{j}$, as well as the sequence of nested subgraphs $\{G_{i}\}$, $i=0,..,m$, spanned by the nodes in corresponding cores, see Fig. 17.

We start by inferring node angular coordinates $i\in cr_{0}$ by maximizing $G_{0}$-specific likelihood ${\rm ln}\mathcal{L}\left[G_{0}\right]$. We then utilize the inferred angles $\{\theta_{i}\}\in cr_{0}$ as initial approximation to maximize ${\rm ln}\mathcal{L}\left[G_{1}\right]$. We continue the angular coordinate inference procedure in the nested fashion to find angular values maximizing ${\rm ln}\mathcal{L}\left[G_{m}\right]$:

[TABLE]

We maximize each log likelihood ${\rm ln}\mathcal{L}\left[G_{\ell}\right]$ iteratively by visiting $G_{\ell}$ nodes in rounds. At each round every node $i$ in $G_{\ell}$ is visited once and placed at $\hat{\theta}_{i}$ maximizing its local log likelihood $\mathcal{L}\left[G_{\mathbb{S}}\right]_{i}$ with respect to the current angular values of other nodes in $G_{\ell}$. The procedure is continued until we arrive at the stable angular configuration:

[TABLE]

where $0<\epsilon\ll 1$ is the precision parameter and $\Delta\hat{\theta}_{i}$ is the angular difference between angular positions of node $i$ in two consecutive rounds. In our experiments we set $\epsilon=10^{-4}$ radians.

The required total number of all-node visit rounds is typically small, of the order of the network average degree. In certain circumstances, e.g., in the case of the global ${\rm ln}~{}\mathcal{L}\left[G_{\ell}\right]$ maximum close to the second largest maximum, the procedure may require a large number of rounds to converge. To avoid these scenarios we limit the maximum number of rounds to $10$ per $G_{\ell}$.

Our experiments indicate that the resulting hyperlink link prediction accuracy is highly sensitive to the correct placement of highest degree nodes. Thus, to further improve angular inference of the most connected nodes, we split the procedure into two parts, $\ell=0,1,..,\lfloor m/2\rfloor$ and $\lfloor m/2\rfloor+1,..,m$, respectively. The first part is repeated independently for ${\rm max\_iter=20}$ times, starting from different initial angle values. For each repetition the resulting ${\rm ln}~{}\mathcal{L}\left[G_{\lfloor m/2\rfloor}\right]$ value is computed. The second part is carried out only once using $\{\theta_{i}\}$ values corresponding to the iteration with largest ${\rm ln}~{}\mathcal{L}\left[G_{\lfloor m/2\rfloor}\right]$ value. Since $\ell=0,1,..,\lfloor m/2\rfloor$ cores are significantly smaller than $\ell=\lfloor m/2\rfloor+1,..,m/2$ cores, the first part is carried out much faster than the second, despite the large number of repetitions.

After each round $\ell$ we perturb the angular coordinates $\hat{\theta}_{i}$, $i\in{\rm cr_{\ell}}$, by adding random noise:

[TABLE]

with amplitude $a(\ell)$, which we decrease linearly as $a(\ell)=\frac{\pi}{4}\left(1-\frac{\ell}{m}\right)+a_{0}$. These coordinate perturbations allow us to avoid getting trapped in local maxima of the log-likelihood function and to arrive to the optimal angles $\{\theta_{i}\}$ faster. We also stress the importance of the nonzero residual noise amplitude of $a_{0}$. In the final $\ell=m$ stage residual noise allows us to effectively ”repel” $k=1$ nodes connected to the same node. Without residual noise at the $\ell=m$ step, all $k=1$ nodes connected to the same node are likely to be placed very close to each other and their common neighbor. As a result, pairs of these $k=1$ nodes will be ranked as the most likely candidates for link prediction, and will adversely affect the hyperlink accuracy. Our experiments indicate that the hyperlink accuracy is not sensitive to specific $a_{0}$ values, as long as $a_{0}\in\left[10^{-6},10^{-3}\right]$. In all our experiments we set $a_{0}=10^{-4}$ radians.

The final part of the embedder algorithm is the series of $20$ coordinate perturbations, following local coordinate inferences in the entire network $G$. This last step often helps to further improve coordinate inference accuracy and, consequently, the accuracy of link prediction. The angular inference procedure is summarized in Alg. 1.

Having sketched the angular inference procedure, we now focus on the individual node placement subroutine. We determine $\hat{\theta}_{i}$ for each node by maximizing the corresponding local log likelihood ${\rm ln}~{}\mathcal{L}\left[G_{\ell}\right]_{i}$. To this end, we split the angular space $\left[-\pi,\pi\right]$ evenly into $\mathcal{O}(N_{\ell})$ regions, where $N_{\ell}$ is the number of nodes in $G_{\ell}$. By placing node $i$ into each of these regions we then identify $\hat{\theta}_{i}$ maximizing its local likelihood. Since ${\rm ln}~{}\mathcal{L}\left[G_{\ell}\right]_{i}$ calculation takes $\mathcal{O}(N_{\ell})$ steps for each $\theta_{i}$ value, it takes $\mathcal{O}(N_{\ell}^{3})$ steps to execute each round $\ell$. As a result, the overall running-time complexity for $m$ layers, $\mathcal{O}(mN^{3})$, is prohibitive for large networks.

To reduce the running time complexity to $\mathcal{O}(m\langle k\rangle N^{2})$, where $\langle k\rangle$ is the average degree of the entire network, we utilize the following approximation, first offered in Ref. Boguñá et al. (2010). If the number of nodes in $G_{\ell}$ is larger than or equal to $500$, for each node we first obtain the rough estimate of $\hat{\theta}_{i}$ by taking into account only its neighboring nodes in $G_{\ell}$. To this end we find the nearly optimal placement $\tilde{\theta}_{i}$ by maximizing

[TABLE]

Since the summation in Eq. (199) goes only through node $i$ neighbors, it now takes $\mathcal{O}(k_{i}N)$ steps to find $\tilde{\theta}_{i}$. Having obtained the initial approximation, we then look for the optimal angle $\hat{\theta}_{i}$ in the neighborhood of $\tilde{\theta}_{i}$ maximizing the full local likelihood ${\rm ln}~{}\mathcal{L}\left[G_{\mathbb{S}}\right]_{i}$, which takes $\mathcal{O}(LN)$ steps, where $L$ is the neighborhood centered at $\tilde{\theta}_{i}$. Specifically, we search for $\hat{\theta}_{i}$ within $L=300\frac{N_{\ell}}{N}$ regions on both sides of $\tilde{\theta}_{i}$, which takes $\mathcal{O}\left(\frac{N_{\ell}^{2}}{N}\right)$ steps, leading to the overall running time complexity of $\mathcal{O}(m\langle k\rangle N^{2})$ steps. The individual node placement subroutine is summarized in Alg. 2.

The outline of the HYPERLINK embedder above is its simplified description omitting a number of important details and presenting some of them slightly differently. The full detailed description of the algorithm exactly as used in this paper is included in its Bitbucket repository cod .

To validate the hyperbolic geometry inference algorithm we compare inferred coordinates in the RHG to its true coordinates. Parameters of the RHG are taken to be $N=5,000$, $\langle k\rangle=10$, $T=0.5$, and $\gamma=2.5$. As seen from Figs. 18 a-c, the accuracy of the angular coordinate inference does not decline significantly for small degree nodes. This is the case, mainly, due to the nested inference with inference cores $cr_{i}$ covering all network nodes, in contrast to the original algorithm of Ref. Boguñá et al. (2010), where cores only cover the most connected nodes.

As seen from Fig. 18 d, Eq. (189) allows for accurate inference of small radial coordinates. At the same time, radial coordinates inference is less accurate for large radial coordinates. To explain this observation we recall that the key assumption in Eq. (189) is that the node degree in the RHG is fully determined by its radial coordinate. In other words, we assume that possible node degree values are narrowly distributed around its expected value, which is given by Eq. (189). This is indeed the case since node degrees are distributed according to the Poisson distribution, Eq. (49). The coefficient of variation of the Poisson distribution, however, is large for small mean values. This leads to significant variation in node degree values in the case of nodes with large radial coordinates, making Eq. (189) inaccurate.

The hyperlink embedder allows for accurate node coordinate inference even in substantially incomplete networks in contrast to other mapping methods, e.g., hypermap Papadopoulos et al. (2015a) and the algorithm by Bläsius et al. Bläsius et al. (2016), which become less accurate in the case of large $T$ values, Fig. 19.

As evidenced by Fig. 18 and, indirectly, by our link prediction results in Secs. III and IV, our hyperbolic inference algorithm is sufficiently accurate for the prediction of missing links on both synthetic and real networks. At the same time, the algorithm does have limitations. First, it is designed to map networks with links removed uniformly at random. The link presence rate $q$ is the required parameter of the algorithm. In cases when the fraction of missing links is unknown, $q$ needs to be estimated and this may lead to less accurate mapping. The second limitation is the algorithm’s running time complexity of $\mathcal{O}\left(N^{2}\right)$ restricting its utility to networks of smaller size. Finally, the third limitation is the analytic estimation of radial coordinates, which is not accurate for small degree nodes. Addressing these limitations is the subject of future work that is expected to further improve the accuracy and the utility of link prediction with hyperbolic geometry.