Learning Based Proximity Matrix Factorization for Node Embedding

Recap from Section 1 that, the personalized PageRank (PPR) $\boldsymbol{\pi}_{u}(v)$ of node $v$ with respect to $u$ is the probability that an $\alpha$-discounted random walk from $u$ stops at node $v$, where an $\alpha$-discounted random walk has $\alpha$ probability to stop at the current node and $(1-\alpha)$ probability to randomly jump to one of its out-neighbors. Define the transition matrix $\boldsymbol{P}=\boldsymbol{D}^{-1}\boldsymbol{A}$ where $\boldsymbol{D}$ is the diagonal matrix such that $\boldsymbol{D}(i,i)$ is the out-degree (resp. degree) of node $i$ if $G$ is directed (resp. undirected), and $\boldsymbol{A}$ is the adjacency matrix. Then, the PPR proximity matrix can be expressed as:

In our trainable proximity measure, instead of fixing the stopping probability at each step to be $\alpha$, we allow the stopping probability $\alpha_{l}$ of the random walk at the $l$-th hop to be trainable. Currently, we assume that the stopping probability $\alpha_{l}$ at the $l$-th step is given and will show how to train $\alpha_{l}$ in Section 3.2. A random walk that follows such learned stopping probability at each step is denoted as a supervised random walk and the proximity derived from the supervised random walk is denoted as supervised PPR. The supervised PPR proximity matrix can be defined as:

To explain, the probability that a supervised random walk stops at exactly the $l$-th hop is $\alpha_{l}\cdot\prod_{k=0}^{l-1}(1-\alpha_{k})$. Thus, by summing up the probability to stop at each hop, we derive the supervised PPR score.

The exact supervised PPR score then can be computed with an iterative manner as shown in Algorithm 1 when $L\rightarrow\infty$. However, we observe that when $L$ is sufficiently large, the probability that the supervised random walk stops with a hop number larger than $L$ is close to zero. Hence, we discard all supervised random walks that stop with a hop larger than $L$. The following theorem shows the quality of calculated supervised PPR scores with Algorithm 1.

All the proofs of our theorems can be found in the appendix. According to our observation, the probability that the supervised random walk stops at a hop greater than $L$ is almost zero when we set $L=15$. Hence, $L$ is set to $15$ in our experiment.

Algorithm 1 makes a tight connection between the supervised PPR proximity matrix and our training parameters $\alpha_{l}$ ($0\leq l\leq L$), which allows the gradients to flow from the derived supervised PPR proximity matrix $\boldsymbol{S}_{L}$ to the training parameters via backpropagation. Let $n$ and $m$ denote the number of nodes and the number of edges, respectively. The time complexity of the algorithm can be bounded by $O(m\cdot n\cdot L)$ since $\boldsymbol{P}$ is a sparse matrix with $m$ entries.

Algorithm 2 shows the pseudo-code of the training process of Lemane. Firstly, it initializes the stopping probability $\alpha_{k}$ for $0\leq k\leq L$ such that the probability a random walk stops at the $l$-th hop follows some standard distribution, e.g., uniform, geometric, or Poisson (Line 1). Next, it uses the initial settings of these stopping probabilities to derive the supervised PPR score by invoking Algorithm 1 (Line 3). Given the supervised PPR scores, it eliminates all the proximity scores that are too small. A threshold $\delta$ is included and all scores smaller than $\delta$ are set to zero (Lines 4-5). Then, a non-linear operation, $\log({\frac{\boldsymbol{S}}{\delta}})$, is applied to the proximity matrix $\boldsymbol{S}$. The proximity scores are divided by $\delta$ to guarantee that each entry after taking the $\log$ will be non-negative (Line 6). Thereafter, a differentiable SVD, e.g., (Indyk et al., 2019), is applied on the new matrix $\boldsymbol{M}$ with input parameter $d$ (Line 7). The SVD obtains two $n\times d$ matrices $\boldsymbol{U}$ and $\boldsymbol{V}$, and a $d\times d$ diagonal matrix $\boldsymbol{\Sigma}$ such that $\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V^{T}}\approx\boldsymbol{M}$. Given the three matrices, $\boldsymbol{U}\sqrt{\Sigma}$ is returned as the first embedding matrix $\boldsymbol{X}$ and $\boldsymbol{V}\sqrt{\Sigma}$ is returned as the second embedding matrix $\boldsymbol{Y}$.

Subsequently, the two embedding matrices are fed into the loss function for the link prediction or node classification (Line 9). The loss functions will be discussed shortly. Then, the stopping probabilities are updated according to the loss function by backpropagation. The training process terminates until the loss function converges (Line 2). Notice that, due to the extremely long computational chain, it is infeasible to write down the explicit form of the gradients. However, like modern deep neural networks, we can use the autograd feature in PyTorch to numerically compute the gradients with respect to the training parameters. Also, the builtin SVD in PyTorch supports to compute the gradient and hence we directly adopt the PyTorch implementation. In what follows, we elaborate on the details of the loss functions.

where $\boldsymbol{x}_{u}$ is the $u$-th row of $\boldsymbol{X}$, $\boldsymbol{y}_{v}$ is the $v$-th row of $\boldsymbol{Y}$, and $d_{out}(u)$ is the out-degree of $u$. Equation 2 indicates that our learned embedding by the supervised random walk should preserve the out-degree information as much as possible.

The second part of the loss function for link prediction is the average cross-entropy over all existing edges:

where $\sigma(x)=1/(1+\exp(-x))$ is the sigmoid function and $\boldsymbol{A}_{u,v}=1$ if edge $(u,v)$ exists in the input graph and $\boldsymbol{A}_{u,v}=0$ otherwise. The final loss function for link prediction is:

where $\beta$ and $\gamma$ are two balancing hyperparameters.

where $\boldsymbol{z}_{u}$ is the concatenated representation of node $u$, $n_{c}$ is the total number of class labels in the graph, and $L_{k}$ is the unnormalized Laplacian matrix of $G_{k}$. The goal of $\mathcal{L^{\prime}}_{1}$ is to minimize pair-wise distances between node pairs with the same class label and maximize pair-wise distances between negative samples.

Next, we employ an activation function to normalize the output embedding to a probability distribution over predicted class labels: $p_{uk}={\textrm{softmax}}(\boldsymbol{Z}\cdot\boldsymbol{W}+\boldsymbol{b})_{uk}$, where $\boldsymbol{W}$ is a $2d\times n_{c}$ fixed mapping matrix generated from uniform distribution, $\boldsymbol{b}$ is the bias term, $p_{uk}$ denotes the probability that node $u$ has class label $k$, and ${\textrm{softmax}}(\boldsymbol{x})_{uk}=\exp(\boldsymbol{x}_{uk})/(\sum_{c=1}^{n_{c}}\exp(\boldsymbol{x}_{uc}))$. Let $\mathcal{Y}$ denote the $n\times n_{c}$ label matrix, where $\mathcal{Y}_{u,k}=1$ if node $u$ has class label $k$ and $\mathcal{Y}_{u,k}=0$ otherwise. Then, the second part of the loss function for node classification is the average cross-entropy over all nodes:

The final loss function for node classification is defined according to $\mathcal{L}^{\prime}_{1}$ and $\mathcal{L}^{\prime}_{2}$ as follows:

where $\beta^{\prime}$ and $\gamma^{\prime}$ are two balancing hyperparameters.

The above-mentioned training process requires calculations on a dense matrix of $O(n^{2})$ size and requires $O(L\cdot n\cdot m)$ running cost, which makes it non-scalable to large graphs. However, such cost seems unavoidable at first glance since we need to obtain the proximity matrix to do backpropagation. After a careful analysis, we make the following two observations to help us avoid the high running costs. Our first key observation is that the parameters that we need to train are only the stopping probabilities at each hop, which are node-independent. That is to say, if we can find a subgraph of the input graph $G$ such that the learned stopping probabilities on the subgraph are identical to that on $G$ for the same task, we can simply learn the parameters on the subgraph with smaller size and apply the learned parameters to the original graph directly, reducing the computational costs. Another observation is that a subgraph of $G$ with similar connectivity should share similar learning stopping probabilities as the input graph on the same task.

Motivated by this, we present our sub-sampling based training method for Lemane on large graphs. Obviously, a straightforward solution is to sample a number $n_{s}$ of nodes and then consider the subgraph containing these $n_{s}$ nodes. However, such a solution severely degrades the connectivity among the nodes. Simple edge sampling strategies will face a similar dilemma which hampers the connectivity among the sampled nodes.

To keep the connectivity among the sampled nodes, we apply a BFS style traversal for subgraph sampling. Our goal is still to sample a subgraph with a constant number $n_{s}$ of nodes. To sample such a subgraph, firstly a source node $u$ is randomly sampled from the input graph $G$. Thereafter, a BFS traversal is applied from the source $u$ to explore the local community of node $u$. If the number of visited nodes by the BFS from $u$ is smaller than $n_{s}$, another node $v$ is randomly sampled as the source to do BFS. The BFS sampling stops as soon as in total $n_{s}$ nodes are visited.

However, the weights trained on a single subgraph might be biased and makes the learned stopping probabilities non-generalizable to the original input graph $G$. To make the learned parameters generalizable to the input graph, we sample multiple subgraphs by the above strategy to learn the parameters. In particular, in each iteration, we sample a subgraph by the BFS strategy and then update the training parameters, i.e., the stopping probabilities, according to the loss functions defined on this subgraph. The loss functions on the subgraph are modified accordingly where the total number $n$ of nodes in Equations 2 and 6 are replaced by sample size $n_{s}$; the total number $m$ of edges in Equation 3 is also replaced by sample size $n_{s}$; the set of nodes with label $k$ in Equation 5 is changed to $V_{Sk}=V_{k}\cap V_{S}$, where $V_{S}$ is the node set of the subgraph, $V_{k}$ is the set of nodes with label $k$ in $G$; and the negative sample set in Equation 5 is changed to $\mathcal{N_{S}}\subset V_{S}\times V_{S}$.

With such a sampling technique, the time complexity of Algorithm 2 can be bounded by $O(h\cdot L\cdot n_{s}^{3})$ where $h$ is the number of training iterations by Algorithm 2 until it converges. Since $n_{s}$ is a controllable constant, we set $n_{s}$ such that the proximity matrix can be fed into the GPU memory for more efficient training.

Our main idea to compute the supervised PPR proximity matrix is to derive a sparsified proximity matrix such that the supervise PPR scores no larger than $\delta$ can be safely discarded. But still, how much cost should we take to derive a sufficiently accurate approximation proximity score? In STRAP (Yin and Wei, 2019), the authors propose a solution to derive the PPR estimations such that $|\boldsymbol{\pi}_{u}(v)-\hat{\boldsymbol{\pi}}_{u}(v)|<\delta$ with a cost of $O(\frac{m}{\delta})$. But, can we further reduce the running cost without sacrificing the embedding quality? Here, we give an affirmative answer. Our solution is inspired by the local graph clustering algorithm Local-Push (Andersen et al., 2006) which returns approximate PPRs with respect to a source $s$ in $O(\frac{1}{\delta})$ running time and guarantees that

on undirected graphs where $d_{out}(v)$ is the degree of node $v$. The Local-Push algorithm suggests that if we only want to compute the approximate PPR scores around the local graph cluster with respect to a source, the running time can be reduced. In our case, the node embedding aims to find nodes that are their representatives and the nodes in their local graph cluster stand as perfect representatives. However, the Local-Push algorithm only works for PPR, not for our supervised PPR. Thus, we present a generalized push algorithm that works for arbitrary stopping probabilities.

Algorithm 3 shows the pseudo-code for our generalized push algorithm. Given a source node $s$, a vector $\boldsymbol{\hat{\pi}}_{s}$ is maintained to store the portion of supervised random walks that has stopped at each node, and is the estimated supervised PPR scores with respect to source $s$. Besides, for each hop $0\leq k\leq L$, where $L$ is the maximum length of a truncated supervised random walk introduced in Section 3.1, an additional residue vector $\boldsymbol{r}_{s}^{(k)}$ is maintained. The vector $\boldsymbol{r}_{s}^{(k)}$ indicates the portion of supervised random walks from $s$ that currently stay at the $k$-th hop but have not stopped yet. Thus, if the residue vectors are all zero, it returns the exact supervised PPR scores. Initially, the residue vectors are all zero except for $\boldsymbol{r}_{s}^{(0)}(s)=1$ (Lines 1-2), indicating that the supervised random walks initially all stay at $s$ and has not stopped yet. Then, if any entry $\boldsymbol{r}_{s}^{(k)}(v)$ in the $L$ residue vectors is above $\delta\cdot d_{out}(v)$ (Line 3), a push operation (Lines 4-7) is invoked. In particular, it first converts $\alpha_{k}\cdot\boldsymbol{r}_{s}^{(k)}(v)$ to $\boldsymbol{\hat{\pi}}_{s}(v)$ (Line 4). To explain, $\alpha_{k}$ portion of the $\boldsymbol{r}_{s}^{(k)}(v)$ random walks stop at the $k$-th hop. Next, the remaining $(1-\alpha_{k})$ portion of the $\boldsymbol{r}_{s}^{(k)}(v)$ random walks randomly jump to the out-neighbors of $v$ (Lines 5-6). Thus, for each $u$ that is an out-neighbor of $v$, the residue $\boldsymbol{r}_{s}^{(k+1)}(u)$ is incremented by $(1-\alpha_{k})\cdot\boldsymbol{r}_{s}^{(k)}(v)/d_{out}(v)$. After the push operation, the residue $\boldsymbol{r}_{s}^{(k)}(v)$ is set to zero (Line 7). The algorithm terminates when there exists no residue $\boldsymbol{r}_{s}^{(k)}(v)$ for any $k$ such that it is larger than $d_{out}(v)\cdot\delta$.

The above analysis shows that our generalized push algorithm can provide identical result quality as the Local-Push algorithm with the identical asymptotic running cost, thus returning high-quality results for the representative nodes of the source node.

Given the generalized push algorithm, we finally show how to output the embedding for the graph. Following STRAP (Yin and Wei, 2019), we compute the supervised PPR on both the input graph $G$ and the transpose graph $G^{T}$ by reversing the direction of each edge of $G$ and set $\boldsymbol{S}(u,v)$ as $\boldsymbol{\pi}_{u}(v)+\boldsymbol{\pi}_{v}^{T}(u)$, where $\boldsymbol{\pi}_{u}(v)$ is the supervised PPR of $v$ with respect to $u$ and $\boldsymbol{\pi}_{v}^{T}(u)$ is the supervised PPR of $u$ with respect to $v$ on the transpose graph $G^{T}$. Note that we do not bring this part into our training phase to reduce the computational costs since we use the same stopping probabilities for both the input graph $G$ and the transpose graph $G^{T}$. Algorithm 4 shows how to compute approximate proximity scores. For each node $u$, we compute the approximate supervised PPR on the input graph $G$ and the transpose graph $G^{T}$ (Lines 3-4). For any approximate supervised PPR score $\boldsymbol{\hat{\pi}}_{u}(v)$, it is added to $\boldsymbol{S}(u,v)$ only if it is larger than the threshold $\delta$. Similarly, each approximate supervised PPR score $\boldsymbol{\hat{\pi}}^{T}_{u}(v)$ on the transpose graph $G^{T}$ is added to $\boldsymbol{S}(v,u)$ only if $\boldsymbol{\hat{\pi}}^{T}_{u}(v)$ is larger than $\delta$ (Lines 5-9). With such a pruning strategy, the proximity matrix $\boldsymbol{S}$ is sparsified to include $O(\frac{n}{\delta})$ non-zero entries. Then, a non-linear operation $\log(\frac{\boldsymbol{S}}{\delta})$, is applied to $\boldsymbol{S}$. Notice that all the entries are divided by $\delta$ before taking the logarithm to guarantee that the values will be non-negative (Line 10). The resulting matrix $\boldsymbol{M}$ is then fed to a sparse SVD to derive final embedding matrices $\boldsymbol{X}$ and $\boldsymbol{Y}$ (Lines 11-12). We have the following theorem for the running time and decomposition quality with respect to $\boldsymbol{M}$.

Following previous work (Yin and Wei, 2019), we treat $\epsilon$ as a constant and use the default setting of builtin SVD implementations. The running time of Algorithm 4 is $O(\frac{n}{\delta}+n\cdot d^{2})$.

Datasets. We test on six real datasets that are used in recent node embedding studies (Perozzi et al., 2014; Grover and Leskovec, 2016; Tsitsulin et al., 2018; Lerer et al., 2019; Yin and Wei, 2019; Yang et al., 2019, 2020). The statistics of these datasets are shown in Table 1. BlogCatalog (Tang and Liu, 2009), Slashdot (Leskovec et al., 2009), TWeibo (twe, ggle) and Orkut (Mislove et al., 2007) are four social networks in which links represent a friendship/following relationship between users. Wikipedia (wik, mark) is a co-occurrence network of words appearing in the Wikipedia dump. Wikivote (Leskovec et al., 2010) is a who-votes-on-whom network on Wikipedia. All datasets and the node labels (if any) can be downloaded from public sources (blo, orks; wik, mark; sna, SNAP; twe, ggle).

• •

  Skip-gram methods: DeepWalk (Perozzi et al., 2014), Node2vec (Grover and Leskovec, 2016), VERSE (Tsitsulin et al., 2018), and InfiniteWalk (Chanpuriya and Musco, 2020);

• •

  Matrix factorization methods: AROPE (Zhang et al., 2018b), RandNE (Zhang et al., 2018a), ProNE (Zhang et al., 2019), NetSMF (Qiu et al., 2019), STRAP (Yin and Wei, 2019), and NRP²²2There were some implementation issues in the released code of NRP on undirected graphs, which was fixed recently by the inventors. (Yang et al., 2020);

• •

  Neural network methods: GraphGAN (Wang et al., 2018) and AW (Abu-El-Haija et al., 2018);

• •

  Other methods: NetHiex (Ma et al., 2018), GraphZOOM³³3We use the default embedding method DeepWalk for GraphZoom. (Deng et al., 2020), Louvain (Bhowmick et al., 2020).

Since the backpropagation in Lemane is complicated, our objective function may easily fall into a local minimum. To tackle this issue, the stopping probabilities are initialized with different distributions. Specifically, we set $\alpha_{0},\alpha_{1},\cdots,\alpha_{L}$ such that the probability that the supervised random walk stops at the $l$-th hop follows a geometric distribution or Poisson distribution with respect to $l$ and we report the best results of Lemane. The parameters of Lemane are optimized by Stochastic Gradient Descent (SGD) optimizer.

Link prediction aims to predict which pairs of nodes are likely to form edges. Following previous work (Yang et al., 2020; Zhang et al., 2018b), we randomly hide $30\%$ of the edges for testing and train the embedding vectors on the rest of the graph. Then, the testing set is generated by including (i) the node pairs corresponding to the $30\%$ removed edges, and (ii) an equal number of node pairs that are not connected by any edge in the original graph $G$. Given a node pair $(u,v)$ in $E_{test}$, we compute a score for $(u,v)$ based on the embedding vectors, and evaluate model performance using precision score.

Following previous work (Yin and Wei, 2019; Yang et al., 2020), for DeepWalk, Node2vec, VERSE, InfiniteWalk, GraphGAN, and Louvain, we use the edge-feature approach introduced in (Ma et al., 2018): (i) randomly select $30\%$ existing edges which are not in $E_{test}$ and the same number of non-existing edges as training set $E_{train}$ on each dataset; (ii) for each node pair $(u,v)\in E_{train}\cup E_{test}$, we concatenate $d$-dimension embedding vectors of node $u$ and that of node $v$; (iii) we consider the $2d$-length vectors as the features of node pairs in $E_{train}$ and feed them into a binary logistic regression classifier; (iv) then the trained classifier is used to perform link prediction on $E_{test}$. For AROPE, RandNE, ProNE, AW, NetHiex, and GraphZOOM, the score of a node pair $(u,v)$ is the inner product of the embedding vector of node $u$ and that of node $v$; for NetSMF, STRAP, NRP, and our Lemane, the score of a node pair $(u,v)$ is the inner product of embedding $\boldsymbol{x}_{u}$ from $\boldsymbol{X}$ and $\boldsymbol{y}_{v}$ from $\boldsymbol{Y}$ (Ref. to Section 3 for the definitions of $\boldsymbol{X}$ and $\boldsymbol{Y}$).

Table 2 reports the performance of Lemane-F against the 15 competitors on three small datasets. Table 3 further reports the performance of Lemane-S against 6 methods which scale to large graphs. For the other 9 methods, they either cannot finish training in 24 hours or run out of memory on the large graphs. As we can observe, our Lemane shows the best performance on 4 social networks where link prediction finds extensive applications. On the Wikipedia dataset, a co-occurrence network of words appearing in the Wikipedia, our Lemane still achieves high performance and is the best method among all matrix factorization methods. Compared to two state-of-the-art matrix factorization methods STRAP and NRP, our Lemane takes the lead by more than $1\%$ on Wikipedia and Slashdot, and takes the lead by almost $3\%$ on Orkut. This demonstrates the effectiveness of our learning based method.

Node classification task aims to predict the label(s) of each node based on the embeddings. As we mentioned in Section 3, we first randomly sample $5\%$ labeled nodes for parameter training. Then the classification task is performed with the following steps: (i) following (Yang et al., 2020), for Lemane and other matrix factorization methods ⁴⁴4For matrix factorization based methods, there was some implementation issues in the evaluation code for node classification on directed graphs. We have rerun the experiment for matrix factorization based methods on directed graphs., we first normalize the embedding vector⁵⁵5We observe some significant improvement on Orkut dataset when normalization is applied. Thus, we take normalization for the embedding vectors on all datasets. $\boldsymbol{x}_{v}$ from $\boldsymbol{X}$ and embedding vector $\boldsymbol{y}_{v}$ from $\boldsymbol{Y}$ of each node $v$, and then concatenate them to get the representation of $v$; (ii) for methods without factorization operation, we use the embedding vector of node $v$ as its representation. Specifically, we randomly split the node and label sets into the training set and testing set and the training ratio varies from $10\%$ to $90\%$. To make a fair comparison, note that if the training ratio is $x\%$, then for Lemane, it will sample an additional $(x-5)\%$ and include the $5\%$ labeled nodes to train the classifiers. Following previous work (Perozzi et al., 2014; Grover and Leskovec, 2016), we employ a one-vs-rest logistic regression classifier implemented by LIBLINEAR (Fan et al., 2008) with default parameters for all methods. Micro-F1 score is used as the evaluation metric for the classification task.

For node classification, we test on four datasets Wikipedia, BlogCatalog, TWeibo, and Orkut, which include label information. Table 4 and Table 5 show the performance of our Lemane-F against all 15 methods on the two small datasets: Wikipedia and BlogCatalog, respectively. Table 6 and Table 7 show the performance of our Lemane-S against 6 methods that scale to TWeibo and Orkut.

We make the following observations. Firstly, our Lemane achieves the best Micro-F1 scores on three datasets BlogCatalog, Tweibo, and Orkut in all of the tested training ratios. Besides, compared to STRAP, which takes PPR without training the stopping probabilities, our Lemane achieves more than $1\%$ lead on Wikipedia datasets and up to $3\%$ on the Orkut dataset. Compared to the second-best matrix factorization method NRP, our Lemane further achieves about $1\%$ lead on the BlogCatalog dataset in all of the tested training ratios.

In summary, the experimental study reveals that our learning based Lemane can learn proximity measures that most suit the task in most scenarios. Compared to other matrix factorization methods, e.g., STRAP (Yin and Wei, 2019) and NRP (Yang et al., 2020), that take personalized PageRank as the proximity measure without learning, our Lemane can train the proximity measure, i.e., the supervised PPR, to gain better performance.