On a Bernoulli Autoregression Framework
for Link Discovery and Prediction

Predicting and recommending links in a future network is a problem of significant interest in many communities with problems of interest ranging from recommending friends Aiello et al., (2012), patent partners Wu et al., (2013) to applications in disease propagation Myers and Leskovec, (2010) and detection of abnormal communication Huang and Lin, (2009). From a higher perspective link prediction is not markedly different from the problem of matrix completion when the network is represented as a adjacency matrix. However, in many cases past link information is predictive of the future. In one such application, Richard et al., (2010) proposed an approach for time dependent graphs where past adjacency matrices and a forecast of a network feature sequence is used to estimate links at a future time. This is a dynamic extension of the static link prediction problem Liben-Nowell and Kleinberg, (2007) with origins in using network properties of a static snapshot to predict links. As a simple example consider the problem of predicting communication between two users over a network, past information is clearly predictive.

In this paper we discuss one variant of the dynamic link prediction problem and in contrast to Richard et al., (2010) present a highly scalable model with an efficient estimation approach. The underlying idea behind our approach is to exploit a sequence of auxiliary networks with feature data to inform links in the network sequence of interest. This is a fairly common scenario in a large corporation with multiple products and services, it can be immensely useful to use network information gathered from a product A to inform usage (via link prediction) on a product B. In the online retail setting¹¹1Example dataset: https://www.kaggle.com/retailrocket/ecommerce-dataset such as Amazon a variety of user-item interactions are possible including “view”, “add-to-cart” and “purchase”. One could generate networks at a user-item level for each type of event, and leverage the “view” network or the “add-to-cart” network to infer connections on the “purchase” network, i.e., when a user eventually buys an item. In such scenarios we have two network sequences, a main sequence that is expected to be sparse and a relatively dense auxiliary sequence with additional feature information that is predictive of a link in the main network. The introduction of the auxiliary sequences serves two purposes, one to provide candidates of links to predict from and a feature set that is possibly predictive of these links. However, the model we propose is general and holds even for a single network sequence with feature information. In the standard matrix completion framework, the observed entries are directly modeled with low rank constraints on the matrix. Instead, we model link probabilities using a Bernoulli parametrization, incorporate an auto-regressive structure to allow for dependence on past data and also include an additional dependence on logit transformed feature data. The equivalence of a low rank constraint in this framework is a relatively under-parametrized model as we shall see in Section 3.

Existing work on link prediction has primarily focused on static networks with methods ranging from matrix factorization Koren, (2008), node neighborhood measures Sarkar et al., (2011); Liben-Nowell and Kleinberg, (2007) to network diffusion Myers and Leskovec, (2010). For an extensive survey see Wang et al., (2015). In contrast (dynamic) link prediction in the presence of a sequence of networks is relatively unexplored. Vu et al., (2011) propose a longitudinal network evolution approach that models edges at time $t$ as a multivariate counting process where the intensity function incorporates edge level features (network derived) using either the multiplicative Cox or additive Aalen form. Our setting is different in that we have two sequences of networks that we can exploit and the design of our estimation method enables us to discover new links which is not possible here. Moreover, the estimation methods require computation of large weight matrices and it is not clear if this can scale to networks with millions of nodes. Additionally, many real world networks tend to have low rank structures due to similar local node characteristics and explicit introduction of these constraints is a desirable feature. Closer in structure to our model is the nonparametric approach described in Sarkar et al., (2012), the edge probabilities are assumed to be Bernoulli distributed with a link function connecting the edge/node probabilities to the node/edge level features. Besides the difference in the problem setting, Sarkar et al., (2012) use a nonparametric estimator for the link function that uses a kernel weighted local neighborhood approach (in time and feature space) to estimate the link function. This requires a search over edge features that are similar to the edge being modeled and can be quite expensive, even with the use of fast search techniques it is unlikely to scale to extremely large graphs.

In the matrix completion line of work Richard et al., (2010) propose a graph feature tracking (GFT) framework for dynamic networks where the entire network at a future step is estimated by solving an appropriately regularized optimization problem. The optimization reflects a trade-off between minimizing the error between the estimated and the last adjacency matrix, a low rank constraint and a term involving the features. This term attempts to minimize the gap between the mapped features for the estimated matrix and the unknown adjacency matrix at the future time. The main intuition being that the evolving features are indicative of links in the sequence of networks. This combination of optimization terms allows the method to uncover links that are not previously observed. See Section 4 for more details of GFT considered in the simulation study. In this work we also introduce a feature regularization term to allow for link discovery but unlike Richard et al., (2010) our approach is highly scalable. While their method involves an expensive SVD operation at every iteration we propose a stochastic gradient based algorithm with each iteration involving simple operations. A later work Richard et al., (2014) extended the GFT framework to allow for features that evolve as a vector auto-regressive process and as such is not computationally feasible for very large networks.

Perhaps most importantly our problem setup is more general and applies to arbitrary binary sequences with mapped features. In Section 5.2 we apply our approach on a hyperlink network from Reddit and show its improved performance over a baseline approach. The Reddit hyperlink network represents a single network sequence and the applicability of our approach demonstrates its potential wide usage on modeling time-dependent binary sequences. The rest of the work is organized as follows. In Section 3 we propose the BAR model with its estimation problem in Section 3.1. We compare the BAR model with baselines (GFT, logistic) in a simulation study in Section 4. On two real datasets in Section 5, one from actual products A and B ²²2To maintain anonymity we don’t reveal the actual product names. and the other from Reddit, we apply the proposed model to show its superior performance in link prediction and discovery. Finally, we conclude the work with discussions in Section 6.

Let the adjacency matrices $\mathbf{\bm{A}}^{(t)}\in\{0,1\}^{n}$ and $\mathbf{\bm{B}}^{(t)}\in\{0,1\}^{N}$ be the sequence of main and auxiliary networks for $t\in\{1,\dots,T\}$. Corresponding to the edges in the auxiliary sequence we also have a sequence of feature tensors $\mathbf{\bm{F}}^{(t)}\in\Re^{N\times N\times d}$ such that, for all $(i,j)\in\mathbf{\bm{B}}^{(t)}$ we have $\mathbf{\bm{F}}^{(t)}_{ij}\in\Re^{d}$. We assume that an edge in a main network sequence exists at least once in the auxiliary sequence and that the number of nodes ($N$ and $n$) is fixed over time. This is not an issue in practice since we can always set $N$ to be the largest values over $t$ and assume all $0$ rows at other times.

There are two important aspects of our model definition. First, we believe that the feature sequence and the history of connections are predictive of a link in the future. As an example, in the context of an email network this could imply that an increased exchange between two nodes is possibly indicative of a future calendar event. Second, the probability of link formation evolves with time and in that sense depends on the entire sequence of edge level features. For any link $\mathbf{\bm{A}}^{(t)}_{ij}$ we have

where $\lambda\in[0,1]$ controls the trade-off between the dependence on the past and the current feature vector for a link $(i,j)$. The model expresses our belief that given the features and link history, a connection at time $t$ is independent (but non-homogeneous) of all other links and that this probability evolves as an auto-regressive model with weights on probability of connection at time $t-1$ and a logit transformed probability induced by the features. We refer to the proposed model as Bernoulli Auto-Regressive (BAR). Note that the recursion of $\mathbf{\bm{Q}}^{(t)}_{ij}$ in Model (1) allows us to express the probability of connection at a given time as a weighted sum of past probabilities. This weighted sum along with a non-linear dependence on the features provides flexibility in modeling arbitrary sequences with observation level features. This is more general than modeling the probabilities as a logit transform of lagged feature vectors. In the choice of $\mathbf{\bm{P}}^{(t)}_{ij}$ we can also use different link functions to extend our model to sequences of discrete outcomes of $K$ categories where $K>2$ (e.g., by using the softmax function), providing that $\mathbf{\bm{P}}^{(t)}_{ij}$ is valued between 0 and 1 so that $\mathbf{\bm{Q}}^{(t)}_{ij}$ is a valid Bernoulli parameter.

For certain applications BAR may be under-parametrized. We can extend the model by having different coefficient vectors in the link probability for different sets of connections. These sets can be determined by using clustering on the features as a pre-processing step. More precisely we can have for $\mathbf{\bm{P}}^{(t)}_{ij}$

where $C^{(t)}_{ij}$ is the cluster membership for link $(i,j)$ at time $t$. We assume only one cluster for edge membership for the rest of the work.

The decaying contribution of historical features is more obvious once we express the recursion of $\mathbf{\bm{Q}}^{(t)}_{ij}$ in Model (1) equivalently as

The longer the time has elapsed since time $s$, the smaller the weight ($\lambda^{t-s}(1-\lambda)$) is applied on the probability $\mathbf{\bm{P}}^{(s)}_{ij}$ and its underlying features from time $s$. The diminishing contribution from features over elapsed time is a valid assumption in many applications (e.g., attribution model). The choice of $\lambda$ governs the decaying rate: larger $\lambda$ makes $\mathbf{\bm{Q}}^{(t)}_{ij}$ more dependent on historical probabilities, whereas smaller $\lambda$ magnifies the dependence on current features. To illustrate the relation between $\mathbf{\bm{Q}}^{(t)}_{ij}$ and $\lambda$, we consider a simplified scenario: we let $\mathbf{\bm{Q}}^{(0)}_{ij}$ be 0.2 and $\mathbf{\bm{P}}^{(t)}_{ij}$ be 0.8 for all $t\in\{1,\ldots,15\}$. In Figure 1 we show how $\mathbf{\bm{Q}}^{(t)}_{ij}$ evolves with $t$ at various $\lambda$ values. At extremes, $\mathbf{\bm{Q}}^{(t)}_{ij}$ takes constant value of $\mathbf{\bm{Q}}^{(0)}_{ij}$ (when $\lambda=1$) or $\mathbf{\bm{P}}^{(t)}_{ij}$ (when $\lambda=0$). For $0<\lambda<1$, the larger $\lambda$ is, the more slowly the contribution of historical features to the $\mathbf{\bm{Q}}^{(t)}_{ij}$ decays and the more slowly $\mathbf{\bm{Q}}^{(t)}_{ij}$ deviates from its initial assignment (0.2) towards the probability derived from current features (0.8).

We simulate data from the generative model (1) and perform several experiments to gain more insight into its behavior and performance compared to competing models. Moreover, since the network sequences and edge features are highly likely to be statistically dependent, we employ settings that are indicative of real-world applications. To generate the first network in the auxiliary sequence we sample Erdős-Renyi graphs with different combinations of $N$ and edge generation probabilities $p$. For subsequent time steps ($T=15$) we vary

For every edge $\mathbf{\bm{B}}^{(t)}_{ij}$ we sample the features $\mathbf{\bm{F}}^{(t)}_{ij}\in\Re^{10}$ such that each element is $\text{Poisson}(\mathbf{\bm{\mu}}_{t})$ distributed, where $\mathbf{\bm{\mu}}_{t}=\mathbf{\bm{\mu}}_{t-1}+\mathbf{\bm{\epsilon}}_{t}\sim\mathcal{N_{\text{10}}}(\mathbf{\bm{0}},0.05\times\mathbf{\bm{I}}_{10})$ is a simple random walk. The coefficient vector $\mathbf{\bm{\beta}}\sim\text{Unif}[0,1]^{10}$ is normalized to have unit Euclidean norm. Finally, the main network sequence $\mathbf{\bm{A}}^{(t)}_{ij}$ is generated by Bernoulli trials with probability $\mathbf{\bm{Q}}^{(t)}_{ij}$ computed using the recursion in (1), where $\mathbf{\bm{Q}}^{(0)}_{ij}$ is initialized as the maximum likelihood estimate of edge probability.

We compare the BAR model with two competing models. The first model is graph feature tracking (GFT) from Richard et al., (2010) which incorporates a feature prediction and regularization procedure for estimating $\mathbf{\bm{A}}^{(T+1)}$. In specific, Richard et al., (2010) define $n\times k$ matrix of features $\mathbf{\bm{\Omega}}$ that can capture the dynamics of $\mathbf{\bm{A}}^{(t)}$, and project $\mathbf{\bm{A}}^{(t)}$ with a linear feature map $\mathbf{\bm{G}}^{(t)}=\mathbf{\bm{A}}^{(t)}\mathbf{\bm{\Omega}}$. With historical feature maps $\{\mathbf{\bm{G}}^{(1)},\ldots,\mathbf{\bm{G}}^{(T)}\}$ Richard et al., (2010) propose using ridge regression to estimate $\hat{\mathbf{\bm{G}}}$ to substitute $\mathbf{\bm{A}}^{(T+1)}\mathbf{\bm{\Omega}}$. In the directed graph setting, we let the features $\mathbf{\bm{\Omega}}=\mathbf{\bm{V}}_{(:,1:k)}\mathbf{\bm{\Sigma}}^{-1}_{(:,1:k)}$ where $\mathbf{\bm{V}}_{(:,1:k)}$ and $\sigma_{1:k}$ are $k$ leading singular vectors and singular values from the SVD of $\mathbf{\bm{A}}^{(T)}=\mathbf{\bm{U}}\mathbf{\bm{\Sigma}}\mathbf{\bm{V}}^{\prime}$, respectively. We solve Problem (10) to get a GFT estimator for $\mathbf{\bm{A}}^{(T+1)}$:

where $\|\mathbf{\bm{S}}\|_{*}$ is a nuclear norm for inducing low rankness in $\mathbf{\bm{S}}$. Since the GFT estimator is not guaranteed to have its values between 0 and 1, we threshold any value of $\mathbf{\bm{S}}$ outside $[0,1]$ to its closer boundary. In addition to GFT, we compare the BAR model to a logistic model with features being averaged historical features up to the current time.

where averaged features across historical time $\mathbf{\bm{F}}^{\text{avg}(t)}$ is defined as follows:

All three models under the comparison have access to historical data though they leverage the data differently.

The comparisons are made over two tasks: a standard setting where networks of different densities are considered, and a setting with varying network dependence for which adjacent network variation is gradually increased. For each scenario, we fix $N$ and $n$ to 10 thousands and $p_{\text{add}}$ to $10^{-3}p$, but vary $p$ and $p_{\text{del}}$ over a range of values. We include all simulation settings in Table 1. For the task with varying network dependence, we gradually increase $p_{\text{del}}$ but fix all other parameters. Larger $p_{\text{del}}$ increases the probability of deleting an existing edge, injecting more variability into networks over time. By varying these parameters we highlight the advantages of our approach to a more traditional matrix completion framework (e.g., GFT from Richard et al., (2010)) and also gain insight into the workings of our method. Under each parameter setting, we generate 10 data samples and fit each model onto the samples. We finely tune each model with training data from the first 14 periods and test on the last period. We report the averaged AUC ROC over 10 repetitions for each model in Table 1. On all tasks the BAR significantly outperforms the GFT and the logistic by achieving higher AUC ROC. When network becomes more varied and more sparse with increasing $p_{\text{del}}$, the challenge increases for all models; the BAR still constantly do better than the other two methods under this scenario.

We exploit the performance of the BAR model on two real experiments. The first experiment uses product B as an auxiliary network to predict and discover connections on product A (the main network). The superior performance of the BAR model demonstrates the connections between the two underlying products. Moreover, the BAR model can fully leverage such connections to infer links on the product of interest (A). Our second experiment focuses on sentiment classification with a Reddit hyperlink network. We construct temporal network sequences for hyperlinks and use the BAR model (along with the baselines) to classify sentiment for the Reddit posts. We show that the BAR model performs at least as well as the logistic baselines even when minimal time dependence is present in network sequence.

We proposed a generalized auto-regressive model for link prediction and discovery with network sequences. The proposed BAR model can exploit signals from the auxiliary network to infer connections on the main network. We also develop an optimization framework for estimating the model, and a stochastic gradient descent algorithm that is scalable to networks with millions of nodes. To demonstrate the superior performance of our model, we apply the BAR model along with baselines on both simulated data and real data, and show significant improvement in BAR’s performance over baselines. Our experiment on the product A’s network shows practical use of the BAR model at scale, while the Reddit hyperlink experiment validates the generalizability of BAR. For future line of research, an extension of the BAR model to continuous time would be useful for many real-world applications. We also plan to investigate the theoretical properties of our model.

The authors thank Sushama Murthy and Richard Johnston for supporting the work.