

\section{Introduction}

The past few years have witnessed the complexity of social networks and the challenges of processing. Therefore, how to learn important information from networks at a lower cost remains a problem. Recently, representation learning on graphs has provided a new paradigm for modeling networks. ~\cite{surveyHamiltonYL17}. Learning latent representations for networks, network embedding, aims to map a network into a latent space. Recent studies have shown that the learned embedding can facilitate various network applications such as vertex classification and link prediction~\cite{line,deepwalk}. These learned embeddings are widely used in various online services. The vital part of item recommendation systems at the E-Commerce platform of Alibaba requires re-embedding when the new users or items arrive online~\cite{recomalibaba}. The embedding must be quickly computed for billions of users and items. There is an urgent need to compute embedding with low latency and low memory overhead with today's billion-scale social networks.

Recent researches have shown the superiority of skip-gram based model, like DeepWalk~\cite{deepwalk}, LINE~\cite{line} and node2vec~\cite{node2vec}. The embedding learned by LINE loses the multi-hop dependencies because it only captures the first- and second-order proximities. DeepWalk and node2vec cannot handle large-scale networks because of expensive computation costs. Network embedding as matrix factorization (NetMF)~\cite{netmf} gives a theoretical foundation that DeepWalk, LINE, and node2vec can be transformed as implicit factorization of a closed-form matrix. NetMF performs better on multi-label vertex classification task than DeepWalk and LINE, which indicates that methods based on matrix factorization can learn effective embedding. However, the matrix to be factorized in NetMF is an $n\times n$ dense one, where $n$ is the number of vertices in the network. Thus it is prohibitively expensive to construct and factorize for billion-scale networks directly. Fortunately, network embedding as sparse matrix factorization (NetSMF)~\cite{netsmf} makes the aforementioned dense matrix sparse through the spectral sparsification theories. However, NetSMF is still relatively costly for time and memory due to a large number of random walks to construct a slightly sparse matrix. ProNE~\cite{prone} proposes to factorize a sparse matrix at first and present a spectral propagation scheme to enhance embedding quality. LightNE~\cite{lightne} combines the strengths of NetSMF and ProNE to propose a cost-effective and scalable embedding system that beats several state-of-art systems such as GraphVite~\cite{graphvite}, PyTorch-BigGraph~\cite{lerer2019pbg}, and ProNE~\cite{prone}. However, the sparsifier construction process of LightNE still faces a huge memory overhead. It limits LightNE to learn effective embedding from networks with hundreds of billions of edges. To obtain high quality, scalable, memory-efficient, and cost-effective network embedding method that can learn promising embedding from networks with billions of nodes and hundreds of billions of edges, we design NetMF+.

\textbf{Our Contribution.} To tackle the challenge of billion-scale networks on a ubiquitous machine equipped with only shared-memory and multi-core CPUs. We build a network embedding system NetMF+ by leveraging the following techniques: 

Firstly, we review the NetMF algorithm, whose first step is eigen-decomposition, followed by forming a dense matrix through eigenvalues and eigenvectors. Then, it generates embedding by performing singular value decomposition (SVD) on the dense matrix. To be precise, we propose a \textbf{fast randomized eigen-decomposition algorithm} that is fast, scalable, and accurate. Experiments show that it is up to 93X faster than the standard \texttt{eigsh} routine in scipy (which is a wrapper to ARPACK~\cite{arpack1998}) without sacrificing accuracy for YouTube dataset, which contains 1.1M vertices and 2.9M edges. We use the fast randomized eigen-decomposition algorithm for the first step of NetMF instead. Then we keep the eigenvectors and eigenvalues and avoid constructing the dense matrix in memory.

Secondly, we combine two randomized matrix factorization techniques: \textbf{sparse-sign random projection} and \textbf{single-pass SVD algorithm}. Instead of performing conventional SVD on dense matrix, we present a sparse-sign randomized single-pass SVD algorithm to complete generating an approximate low-rank result without constructing an $n\times n$ dense matrix in memory. The spectral propagation strategy proposed in ProNE can be used to enhance the quality of embedding.

Thirdly, Considering the CSR format of a network with hundreds of billions of edges will be a huge overhead for memory, we design \textbf{a sparse matrix matrix multiplication operation} based on state-of-the-art shared-memory graph processing techniques. It's used to compress graph and support parallel graph operation. While ensuring accuracy and efficiency, memory-efficient is realized. The main computation of NetMF+ is basic linear algebra. Thus we can use Intel Math Kernel Library (MKL) for these computations. These techniques enable NetMF+ to be the fastest and memory-efficient network embedding system compared to other state-of-the-art methods. In particular, we can embed one of the largest publicly available graphs, the Web Data Commons hyperlink 2012 graph, with over 200 billion edges by NetMF+.
In summary, the main advantages of our NetMF+ are:
\begin{enumerate}[label=(\textbf{\arabic*}),topsep=1pt,itemsep=0pt,parsep=0pt,leftmargin=20pt]
	\item \textbf{Accurate:}
	NetMF+ achieves the highest accuracy in vertex classification task and link prediction task under the same time budget.
	
	\item \textbf{Efficiency:}
	NetMF+ shows the fastest speed compared to state-of-art methods in most cases, according to experiments.
	
	\item \textbf{Memory-efficient:}
	NetMF+ has the minimal memory cost compared to other network embedding systems designed for a single machine.
	
	\item \textbf{Scalable:}
	NetMF+ can embed networks with 200 billion edges within 1.5 hours.
	
\end{enumerate}
\section{Related Work}
In this section, we review the related work of network embedding and randomized matrix factorization.

\textbf{Network Embedding.} Network embedding has been comprehensively studied over the past decade. In general, recent work about network embedding can be divided into three categories. The first category is based on skip-gram methods, inspired by word2vec~\cite{word2vec}, including DeepWalk~\cite{deepwalk}, Line~\cite{line}, node2vec~\cite{node2vec}. These methods rely on stochastic gradient descent to optimize a logistic loss. The second category is based on deep learning methods, including GCN~\cite{gcn}, GAT~\cite{gat}, GIN~\cite{gin}, GraphSAGE~\cite{surveyHamiltonYL17} and PinSAGE~\cite{pinsage}. They are beyond the scope of this paper because we focus on graphs with no additional information and pre-train general network embeddings in an unsupervised manner. In contrast, these methods pre-train in a supervised manner with vertex attributes. The third category method is based on matrix factorization, using SVD to generate the best low-rank approximation~\cite{eckart1936approximation}. Examples of methods that belong to this category include GraRep~\cite{grarep}, HOPE~\cite{hope}, NetMF~\cite{netmf}, NetSMF~\cite{netsmf}, ProNE~\cite{prone}, LightNE~\cite{lightne}. Several embedding systems for large graphs have been developed. A CPU-GPU hybrid network embedding system, GraphVite~\cite{graphvite}, is developed based on DeepWalk~\cite{deepwalk} and LINE~\cite{line}. In Graphvite, CPU is used to perform graph operation, and GPU is used to compute linear algebra. Nevertheless, the GPU memory is a disadvantage when processing billion scale networks, limiting widespread use. Based on DeepWalk~\cite{deepwalk} and LINE~\cite{line}, a distributed memory system PyTorch-BigGraph~\cite{lerer2019pbg} has been proposed. It achieves load balancing by graph partition and synchronization through a shared parameter server. In this work, we propose NetMF+, which leverages the merit of NetMF and addresses its limitation in speed and memory overhead.

\textbf{Randomized Matrix Factorization.} Randomized matrix computation has gained significant increases in popularity as the data is becoming larger and larger. Randomized SVD can be an alternative to conventional SVD methods. Because randomized SVD involves the same or fewer floating-point operations and is more efficient for truly large high-dimensional data sets, by exploiting modern computing architectures~\cite{halko2011finding,randnla}. frPCA~\cite{frpca} is developed for large-scale sparse matrix with better performance compared with conventional PCA algorithms. Randomized SVD has shown its superiority in recent network embedding methods such as NetSMF, ProNE, and LightNE. Over the past few years, randomized matrix approximation for streaming data has been extensively studied~\cite{tropp2017practical}. Several single-pass SVD algorithms~\cite{yu2017ijcai,tropp2017practical,boutsidis2016optimal,tropp2019streaming} with different performances are introduced, which satisfy the streaming data scenario and only need to visit the target matrix once. The efficiency and memory limitation of NetMF will be solved by randomized matrix factorization.
\section{Preliminaries}
The Matlab conventions are used for specifying row/column indices of a matrix and some operations on sparse matrix.
The list of notations used in this paper can be found in Table~\ref{tab:notation}. The problem of the network embedding is as follows: Given an undirected network $G = (V,E,\vect{A})$ with $V$ as the vertex set of $n$ vertices, $E$ as the edge set of $m$ edges, and $\vect{A}$ as the adjacency matrix, the goal is to learn an embedding matrix $\vect{E}\in\mathbb{R}^{n\times d}$ that row $i$ captures the structural property of $i$-th vertex. The embedding matrix can be fed into downstream applications such as multi-label vertex classification and link prediction.
\begin{table}[tb]
 
  \setlength{\abovecaptionskip}{0.05 cm}
  \setlength{\belowcaptionskip}{-0.25 cm}
  \centering
  \caption{Notation used throughout this paper.}
  \label{tab:notation}\small{
  \begin{tabular}[htbp]{@{}l@{~}|r|@{~}l|@{~}r@{}}
    \toprule
    \textbf{Notation} & \textbf{Description} & \textbf{Notation} & \textbf{Description}                    \\  \midrule
    $G$               & input network        & $b$               & \#negative samples                      \\
    $V$               & vertex set, $|V|=n$  & $T$               & context window size                     \\
    $E$               & edge set, $|E|=m$    & $\vect{E}$             & $n\times d$ embedding matrix            \\
    $\vect{A}$        & adjacency matrix     & $d$               & embedding dimension           \\
    $\mathrm{vol}(G)$ & volume of $G$         & $\vect{U}_k$         & truncated eigenvectors            \\
    $\vect{D}$        & degree matrix        & $\vect{\Lambda}_k$ & truncated eigenvalues \\
    $\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}$  & modified Laplacian& $k$ & rank parameter     \\
    $\mathcal{L}$ & normalized  Laplacian& $q$ & power parameter\\
    $\vect{M}$ & $\sum_{r=1}^{T}(\vect{D}^{-1}\vect{A})^{r}\vect{D}^{-1}$ & $s_i$ & oversampling parameter\\
    $\vect{M}^{\prime}$ & NetMF matrix & $z$ & column sparsity\\
    \bottomrule
  \end{tabular}
  }
\end{table}


\textbf{NetMF.} One of the important network embedding methods based on matrix factorization is NetMF~\cite{netmf}. Formally, when it comes to an unweighted, undirected graph $G$, $\vect{A}$ is the adjacent matrix. The element in $\vect{A}$ corresponding to each edge is $1$, and the rest is $0$. $\vect{D}$ is the diagonal matrix, where the $i$-th diagonal entry contains the degree of the $i$-th vertex. The main contribution of NetMF is to prove that DeepWalk is approximately factorizing the matrix
\begin{equation}
	\label{eq1}
	\small
		\vect{M}^{\prime} = \mathrm{trunc\_log}^{\circ}\left(\frac{\mathrm{vol}(G)}{bT}\vect{M}\right) 
		=  \mathrm{trunc\_log}^{\circ}\left(\frac{\mathrm{vol}(G)}{bT}\sum_{r=1}^{T}(\vect{D}^{-1}\vect{A})^{r}\vect{D}^{-1}\right),
\end{equation}
where $T$ is the length of random walks, $b$ is the number of negative samples, and $\mathrm{vol}(G) = 2m$ is the total number of edges in $G$. $\mathrm{trunc\_log}^{\circ}$ is the element-wise truncated logarithm function ($\mathrm{trunc\_log}^{\circ}(x)=\mathrm{max}({0,\mathrm{log}(x)})$) and $\mathrm{log1p}(\mathrm{log1p}(x)=\mathrm{log}(1+x))$ is another optional function that has the similar effect as $\mathrm{trunc\_log}^{\circ}$~\cite{netsmf}. NetMF also shows that LINE approximately factorizes a matrix in the same form but with $T=1$. The bottleneck of factorizing the matrix in Eq.~(\ref{eq1}) is that $\vect{M}$ tends to be a dense one as the increase of $r$, and constructing the matrix is cost-prohibitive for the enormous amount of memory requirement. In addition, computing the matrix power in Eq.~(\ref{eq1}) involves dense matrix multiplication, which costs $O(n^3)$ time.
To reduce time cost in Eq.~(\ref{eq1}) for large $T$ (e.g., the default setting $T=10$ in DeepWalk), NetMF proposes to perform truncated eigen-decomposition on $\vect{D}^{-\frac{1}{2}}\vect{A}\vect{D}^{-\frac{1}{2}}$ and approximates Eq.~(\ref{eq1}) to
\begin{equation}
	\label{eq2}
	\small
	\begin{split}
		\tilde{\vect{M}}^{\prime} &
		\approx\mathrm{trunc\_log}^{\circ}\left(\frac{\mathrm{vol}(G)}{bT}\tilde{\vect{M}}\right) \\
		&=\mathrm{trunc\_log}^{\circ}\left(\frac{\mathrm{vol}(G)}{bT}\vect{D}^{-\frac{1}{2}}\vect{U}_k\left(\sum_{r=1}^{T}\vect{\Lambda}_k^{r}\right)\vect{U}_k^{\top}\vect{D}^{-\frac{1}{2}}\right).
	\end{split}
\end{equation}
However, the approximated matrix $\tilde{\vect{M}}$ is still dense, making it impossible to construct and factorize. It is worth noting that the truncated logarithm is very important to embedding quality and cannot be omitted, otherwise there exists a shortcut to the factorization without constructing the dense matrix $\tilde{\vect{M}}$, similar to NPR~\cite{nrp} and RandNE~\cite{randne}.

\textbf{ProNE.} ProNE~\cite{prone} proposes to perform truncated singular value decomposition~(SVD) on a sparse matrix which each entry defined to be $\vect{M}_{uv}=\mathrm{log}\left(\frac{\vect{A}_{uv}}{\vect{D}_u}\frac{\sum_{j}(\sum_{i}\vect{A}_{ij}/\vect{D}_i)^{\alpha}}{b(\sum_{i}^{}\vect{A}_{iv}/\vect{D}_i)^{\alpha}}\right)$, where $b = 1$ and $\alpha=0.75$. Therefore, it produces an initial embedding matrix $\vect{E}$. Consequently, ProNE proposes to improve the embedding performance by spectral propagation strategy, i.e., $\sum_{r=0}^{p}c_r\mathcal{L}^{r}\vect{E}$, where $c_r$ is coefficients of Chebyshev polynomials, $p$ is spectral propagation steps (the default setting is 10) and  $\mathcal{L} = \vect{I}-\vect{D}^{-1}\vect{A}$ is normalized graph Laplacian matrix ($\vect{I}$ is the identity matrix). Because of the universality of spectral propagation scheme, we can further use it to enhance our embedding.

\textbf{The Basic Randomized SVD.} Previous work has witnessed the advantages of randomized matrix factorization methods for solving low-rank matrix approximation~\cite{randnla,halko2011finding,musco2015nips} and network embedding problems~\cite{lightne,netsmf,prone}. The basic randomized SVD algorithm can be described as Algorithm~\ref{alg:basic_rsvd}.
\begin{algorithm}
	\caption{The basic randomized SVD}
	\small
	\label{alg:basic_rsvd}
	\SetAlgoLined
	\KwIn{Input matrix $\vect{X}\in\mathbb{R}^{m\times n}$, eigen-decomposition rank $k$, power iteration parameter $q$, oversampling parameter $s$}
	\KwOut{singular vectors $\vect{U}\in\mathbb{R}^{m\times k},\vect{V}\in\mathbb{R}^{n\times k}$, singular values $\vect{\Sigma}\in\mathbb{R}^{k\times k}$ } $\vect{\Omega}=\mathrm{randn}(n,~k+s),\vect{Y}=\vect{X}\vect{\Omega}$\;
	$\vect{Q} = \mathrm{orth}(\vect{Y})$\;
	\For{$i=1,2,...,q$}{
		$\vect{G} = \mathrm{orth}(\vect{X}^{\top}\vect{Q})$\;
		$\vect{Q} = \mathrm{orth}(\vect{X}\vect{G})$\;
	}
	$\vect{B}=\vect{Q}^{\top}\vect{X}$\;
	$[\hat{\vect{U}},\hat{\vect{\Sigma}},\hat{\vect{V}}] = \mathrm{svd}(\vect{B})$\;
	$\vect{U}=\vect{Q}\hat{\vect{U}}(:,1:k),\vect{\Sigma}=\hat{\vect{\Sigma}}(1:k,1:k),\vect{V}=\hat{\vect{V}}(:,1:k)$\;
	\KwRet $\vect{U},\vect{\Sigma},\vect{V}$\;
\end{algorithm}
The basic randomized SVD can be used to compute approximate truncated SVD of an input matrix $\vect{X}$. In Alg.~\ref{alg:basic_rsvd}, $\vect{\Omega}$ is a Gaussian i.i.d matrix. The randomized low-rank approximation mainly relies on the random projection to identify the subspace capturing the dominant actions of target matrix $\vect{X}$. It can be realized by multiplying $\vect{X}$ with a random matrix on its right side to obtain the subspace's orthogonal basis matrix $\vect{Q}$. With the sketch matrix $\vect{Q}$, we can have approximation $\vect{X}\approx\vect{QB}=\vect{Q}\vect{Q}^{\top}\vect{X}$. The approximate truncated SVD of $\vect{X}$ can be obtained by performing SVD on matrix $\vect{B}$. In Alg.~\ref{alg:basic_rsvd}, steps 3\textasciitilde6 is the power iteration scheme that can improve the accuracy of the approximation. The "$\mathrm{orth}$" in Alg.~\ref{alg:basic_rsvd} is the orthonormalization operation that can be used to alleviate the round-off error in floating-point computation. QR factorization usually is a choice for orthonormalization operation. 
\section{The NetMF+ Based on Randomized Matrix Factorization}
In this section, we develop a network embedding algorithm NetMF+ based on NetMF and fast random matrix factorization, which is approximate to Eq.~(\ref{eq1}). First, we propose a fast randomized eigen-decomposition algorithm to generate $\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top}$ which is approximated to modified Laplacian matrix. Second, we propose to compute embedding by sparse-sign randomized single-pass SVD algorithm, which avoids constructing the dense matrix $\tilde{\vect{M}}^{\prime}$ in memory. Finally, spectral propagation is used to improve the embedding quality, and Graph Based Benchmark Suite(GBBS)~\cite{gbbs} is adopted as the framework to reduce the peak memory cost.
\subsection{Fast Randomized Eigen-Decomposition}
We first introduce the fast randomized eigen-decomposition algorithm. Based on Alg.~\ref{alg:basic_rsvd}, the approximation $\vect{X}\approx\vect{QB}$ is not suitable for solving eigen-decomposition. According to~\cite{halko2011finding}, the symmetric approximation formula should be $\vect{X}\approx\vect{QS}\vect{Q}^{\top}$ and truncated eigen-decomposition result of $\vect{X}$ can be derived by performing eigen-decomposition on $\vect{S}$. Combining the techniques of power iteration scheme, the symmetric approximation formula and acceleration strategy~\cite{frpca}, fast randomized eigen-decomposition can be described as Algorithm~\ref{alg:frevd}.
\begin{algorithm}[b]
	\caption{Fast randomized eigen-decomposition }
	\label{alg:frevd}
	\small
	\SetAlgoLined\DontPrintSemicolon
	\SetKwFunction{proc}{$\mathrm{freigs}$}
	\SetKwProg{myproc}{Procedure}{}{}
	\myproc{\proc{$\vect{X},k,s_1,q$}}{
		$\vect{\Omega}=\mathrm{randn}(\mathrm{size}(\vect{X},~2),~k+s_1),\vect{Y}=\vect{X}\vect{\Omega}$\;
		$[\vect{Q},\sim,\sim] = \mathrm{eigSVD}(\vect{Y})$\;
		\For{$i=1,2,...,q$}{
			$[\vect{Q},\sim,\sim] = \mathrm{eigSVD}(\vect{X}\vect{X}\vect{Q})$\;
		}
		$\vect{S}=\vect{Q}^{\top}\vect{X}\vect{Q}$\;
		$[\hat{\vect{U}},\vect{\Lambda}] = \mathrm{eig}(\vect{S})$\;
		$\vect{U}=\vect{Q}\hat{\vect{U}}$\;
		\KwRet $\vect{U}(:,1:k),\vect{\Lambda}(1:k,1:k)$\;
	}
\end{algorithm}
In Alg.~\ref{alg:frevd}, the "$\mathrm{eigSVD}$" is used as the orthonormalization operation. According to~\cite{frpca}, compared with "$\mathrm{qr}$", "$\mathrm{eigSVD}$" is much faster especially when $n\gg(k+s_1)$. But there exists a shortcoming that numerical issues may arise if matrix input to "$\mathrm{eigSVD}$" does not have full column rank. Fortunately, $\vect{D}^{-\frac{1}{2}}\vect{A}\vect{D}^{-\frac{1}{2}}$ is a symmetric matrix which eigenvalues are in the range of $[-1,1]$ according to~\cite{netmf}. Therefore, the orthonormalization operation in step 3 and step 5 is stable. Setting $l=k+s_1$, the step 2\textasciitilde3 in Alg.~\ref{alg:frevd} cost $O(ml+nl^2)$, where $n$ is the number of vertices in $\vect{A}$ and $m$ is the number of edges. According to~\cite{frpca}, step 4\textasciitilde6 cost $O(q(ml+nl^2))$, step 7 cost $O(ml+nl^2)$, step 8 cost $O(l^3)$ and step 9 cost $O(nl^2)$. Thus the total time cost is $O(q(ml+nl^2))$. Considering that $l\ll n$, $s_1$ and $q$ is a small number, the time cost of Alg.~\ref{alg:frevd} is $O(mk+nk^2)$, which is consistent with \texttt{eigsh}~\cite{arpack1998}. However, benefitting from the high parallelism of randomized matrix algorithm, Alg.~\ref{alg:frevd} show up to 93X speedup ratio when compared to \texttt{eigsh} with almost same result in the experiment.

For a large scale matrix $\vect{D}^{-\frac{1}{2}}\vect{A}\vect{D}^{-\frac{1}{2}}$, the eigenvalues decay extremely slow and \texttt{eigsh} is hard to converge. There also exists the dilemma that the approximate error is unsatisfactory for randomized algorithm to deal with such a matrix~\cite{halko2011finding}. To alleviate the accuracy problem of matrix approximation, we propose to perform Alg.~\ref{alg:frevd} on a modified Laplacian matrix $\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}$ instead, where $\alpha\in(0,0.5]$. Therefore, we have
$\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top}\approx\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}$.
It means $\vect{D}^{-\frac{1}{2}}\vect{A}\vect{D}^{-\frac{1}{2}}$ is computed as $\vect{D}^{-\frac{1}{2}+\alpha}\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top}\vect{D}^{-\frac{1}{2}+\alpha}$  approximately. The approximation error has an upper bound described by the following theorem, whose proof can be found in Appendix.
\begin{theorem}
	Let $\Vert\cdot\Vert_F$ be the matrix Frobenius norm. Then 
\begin{equation*}
\small
\label{eq:thm1}
	\left\Vert\vect{D}^{-\frac{1}{2}}\vect{A}\vect{D}^{-\frac{1}{2}}-\vect{D}^{-\frac{1}{2}+\alpha}\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top}\vect{D}^{-\frac{1}{2}+\alpha}\right\Vert_F\leq (1+\varepsilon)\sqrt{\sum_{j=k+1}^{n}|\lambda_j|^2} ~,
\end{equation*}
with high probability. Here $|\lambda_j|$ is the $j$-th largest absolute value of eigenvalue of $\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}$.
\end{theorem}
Next, the matrix $\vect{M}$ will be approximated to $\hat{\vect{M}}$:
\begin{equation}
	\label{eq3}
	\small
	\hat{\vect{M}}=\vect{D}^{-1+\alpha}\vect{U}_k\vect{\Lambda}_k\left(\sum_{r=1}^{T}\vect{K}^{r-1}\right)\vect{U}_k^{\top}\vect{D}^{-1+\alpha},
\end{equation}
where $\vect{K}=\vect{U}_k^{\top}\vect{D}^{-1+2\alpha}\vect{U}_k\vect{\Lambda}_k$ is a $k\times k$ matrix, thus the computation of $\sum_{r=1}^{T}\vect{K}^{r-1}$ is cheap. For convenience, we denote $\vect{F} = \vect{D}^{-1+\alpha}\vect{U}_k$ and $\vect{C} = \vect{\Lambda}_k\left(\sum_{r=1}^{T}\vect{K}^{r-1}\right)$. Then Eq.~(\ref{eq1}) can be approximated by
\begin{equation}
	\label{eq4}
	\small
	\hat{\vect{M}}^{\prime}=\mathrm{trunc\_log}^{\circ}\left(\hat{\vect{M}}\right)
	=\mathrm{trunc\_log}^{\circ}\left(\frac{\mathrm{vol}(G)}{bT}\vect{FC}\vect{F}^{\top}\right).
\end{equation}

\subsection{Sparse-Sign Randomized Single-Pass SVD}
In this subsection, we will introduce how to generate embedding $\vect{E}\in\mathbb{R}^{n\times d}$ without constructing the dense matrix $\hat{\vect{M}}^{\prime}$ when we have $\vect{F}\in\mathbb{R}^{n\times k}$ and $\vect{C}\in\mathbb{R}^{k\times k}$. We first introduce the definition of the sparse-sign matrix. The sparse-sign matrix with similar performance to the Gaussian matrix is another type of randomized dimension reduction map~\cite{tropp2019streaming}. Assuming that sparse-sign matrix $\vect{\Psi}\in \mathbb{R}^{n\times h}$, where $h=d+s_2$ and $s_2$ is oversampling parameter, we fix a column sparsity parameter $z$ in the range $2\leq z\leq n$. We independently generate the columns of the matrix at random. For each column, we draw $z$ i.i.d random signs and place them in $z$ uniformly random coordinates. According to~\cite{tropp2019streaming}, $z=\mathrm{min}(n,8)$ will usually be a good choice. Algorithm~\ref{alg:ssmat}  describes how to generate the sparse-sign matrix.
\begin{algorithm}[t]
	\caption{Generate sparse-sign matrix}
	\label{alg:ssmat}
	\small
	\SetAlgoLined\DontPrintSemicolon
	\SetKwFunction{proc}{$\mathrm{gen\_sparse\_sign\_matrix}$}
	\SetKwProg{myproc}{Procedure}{}{}
	\myproc{\proc{$n,h,z$}}{
		$\vect{\Psi}=\mathrm{zeros}(n,~h)$\;
		$\vect{p} = \mathrm{zeros}(z*h,~1)$\;
		\For{$j=1,2,...,h$}{
			$\vect{p}((z(j-1)+1):zj)= \mathrm{randperm}(n,z)$\;
			$\vect{\Psi}(\vect{p}((z(j-1)+1):zj),j) = \mathrm{sign}(\mathrm{randn}(z,1))$\;
		}
		$\vect{p} = \mathrm{unique}(\vect{p})$\;
		\KwRet $\vect{\Psi}, \vect{p}$\;
	}
\end{algorithm}
After generating the sparse-sign matrix, we need to use it as the randomized matrix, which multiplies $\hat{\vect{M}}^{\prime}$ on its right side to obtain sketch matrix $\vect{Y}$. 
Considering that a column of sparse-sign matrix $\vect{\Psi}$ only has $z$ nonzeros, it will generate at most $z*h$ coordinates in the range of $[1,n]$. Therefore, we can perform row sampling according to these unique coordinates $\vect{p}$. Assuming that $[v,\sim]=\mathrm{size}(\vect{p},1)$, which $v\leq z*h$. One can immediately observe that we can have a sampling matrix $\vect{F}(\vect{p},:)^{\top}$. Therefore, the computation $\vect{Y} = \hat{\vect{M}}^{\prime}\vect{\Psi}$ is accomplished by
\begin{equation}
	\label{eq5}
	\small
	\vect{Y} = \mathrm{trunc\_log}^{\circ}\left(\frac{\mathrm{vol}(G)}{bT}\vect{F}\vect{C}\vect{F}\left(\vect{p},:\right)^{\top}\right)\vect{\Psi}.
\end{equation}
The time cost of Eq.~(\ref{eq5}) is $O(nk^2+nvk+nvh)$, and the memory cost is $O(nk+nv)$. When calculating $\vect{F}\vect{C}\vect{F}\left(\vect{p},:\right)^{\top}$ for a network with billion vertices, it will introduce $O(nv)$ memory cost, which is still a disaster. To solve this, we adopt the batch matrix matrix multiplication, selecting the fixed-size rows of $\vect{F}$ in turn to complete multiplication. Therefore, the memory cost will be reduced to $O(nk)$. After orthogonalizing $\vect{Y}\in\mathbb{R}^{n\times h}$, we have an orthogonal basis matrix $\vect{Q}\in\mathbb{R}^{n\times h}$ which contains the information of $\hat{\vect{M}}^{\prime}$. 

There is still a dilemma that we cannot generate $\vect{B}=\vect{Q}^{\top}\hat{\vect{M}}^{\prime}$, where $\vect{Q}$ is a dense matrix.
To generate the embedding of $\hat{\vect{M}}^{\prime}$, we leverage a randomized single-pass SVD algorithm~\cite{tropp2019streaming}. The randomized single-pass SVD algorithm draws multiple sketch matrices that capture the row and column dominant actions of matrix $\hat{\vect{M}}^{\prime}$ and computes SVD based on these sketch matrices. In~\cite{tropp2019streaming}, four random matrices
$\vect{\Phi}\in\mathbb{R}^{x\times n}, \vect{\Psi}\in\mathbb{R}^{l\times x},\vect{\Theta}\in\mathbb{R}^{y\times n}, \vect{\Pi}\in\mathbb{R}^{l\times y}
$ are drawn for target matrix $\vect{X}\in\mathbb{R}^{n\times l}$. Then three sketch matrices
$\vect{K}=\vect{\Phi}\vect{X}$, 
$\vect{Y}=\vect{X}\vect{\Psi}$ and
$\vect{Z}=\vect{\Theta}\vect{X}\vect{\Pi}$
are generated respectively. By the way, $\hat{\vect{M}}^{\prime}$ is symmetric, which indicates that $\vect{K}$ and $\vect{Y}$ can be the same and $\vect{Z}$ can be computed with one random matrix. Briefly, we need two random matrices $\vect{\Psi},\vect{\Pi}$ to generate two sketch matrices $\vect{Y}=\hat{\vect{M}}^{\prime}\vect{\Psi}$ and $\vect{Z}=\vect{\Pi}^{\top}\hat{\vect{M}}^{\prime}\vect{\Pi}$.  With $\vect{Q}=\mathrm{orth}(\vect{Y})$ and $\hat{\vect{M}}\approx\vect{Q}\left(\vect{Q}^{\top}\hat{\vect{M}}^{\prime}\vect{Q}\right)\vect{Q}^{\top}$, we can have	$\vect{Z}=\vect{\Pi}^{\top}\hat{\vect{M}}^{\prime}\vect{\Pi}\approx\left(\vect{\Pi }^{\top}\vect{Q}\right)\left(\vect{Q}^{\top}\hat{\vect{M}}^{\prime}\vect{Q}\right)\left(\vect{Q}^{\top}\vect{\Pi}\right)$.
Then we get the core approximation matrix
$\vect{S}=\left(\vect{\Pi}^{\top}\vect{Q}\right)^{\dagger}\vect{Z}\left(\vect{Q}^{\top}\vect{\Pi}\right)^{\dagger}\approx\vect{Q}^{\top}\hat{\vect{M}}^{\prime}\vect{Q}$ by solving least-squares problems. Finally, it will form a low rank approximation of the target matrix $\hat{\vect{M}}^{\prime}$ via
$\hat{\vect{M}'}\approx\vect{Q}\vect{S}\vect{Q}^{\top}$,
and the approximate truncated SVD of $\hat{\vect{M}}^{\prime}$ can be derived from performing SVD on $\vect{S}$.

With the sparse-sign matrix as the random matrix, the sparse-sign randomized single-pass SVD algorithm for $\hat{\vect{M}}^{\prime}$ can be described as Algorithm~\ref{alg:ssrspSVD}.
\begin{algorithm}
	\caption{Sparse-sign randomized single-pass SVD}
	\label{alg:ssrspSVD}
	\small
	\SetAlgoLined\DontPrintSemicolon
	\SetKwFunction{proc}{$\mathrm{sparse\_sign\_rand\_single\_pass\_SVD}$}
	\SetKwProg{myproc}{Procedure}{}{}
	\myproc{\proc{$\vect{F},\vect{C},d,s_2,s_3,z$}}{
		$[\vect{\Psi},\vect{p}]=\mathrm{gen\_sparse\_sign\_matrix}(\mathrm{size}(\vect{F},1),d+s_2,z)$\;
		$\vect{Y} = \mathrm{trunc\_log}^{\circ}\left(\frac{\mathrm{vol}(G)}{bT}\vect{FC}\vect{F}\left(\vect{p},:\right)^{\top}\right)\vect{\Psi}$\;
		$[\vect{Q},\sim,\sim] = \mathrm{eigSVD}(\vect{Y})$\;
		$[\vect{\Pi},\vect{p}]=\mathrm{gen\_sparse\_sign\_matrix}(\mathrm{size}(\vect{F},1),d+s_3,z)$\;
		$\vect{Z} = \vect{\Pi}^{\top} \mathrm{trunc\_log}^{\circ}\left(\frac{\mathrm{vol}(G)}{bT}\vect{F}\left(\vect{p},:\right)\vect{C}\vect{F}\left(\vect{p},:\right)^{\top}\right)\vect{\Pi}$\;
		$\vect{S}=(\vect{\Pi}^{\top}\vect{Q})^{\dagger}\vect{Z}(\vect{Q}^{\top}\vect{\Pi})^{\dagger}$\;
		$[\hat{\vect{U}},\hat{\vect{\Sigma}},\hat{\vect{V}}] = \mathrm{svd}(\vect{S})$\;
		$\vect{U}=\vect{Q}\hat{\vect{U}}(:,1:d),\vect{\Sigma}=\vect{Q}\hat{\vect{\Sigma}}(:,1:d),\vect{V}=\vect{Q}\hat{\vect{V}}(:,1:d)$\;
		\KwRet $\vect{U},\vect{\Sigma},\vect{V}$\;
	}
\end{algorithm}
In Alg.~\ref{alg:ssrspSVD}, step 4 use "$\mathrm{eigSVD}$" in Alg.~\ref{alg:frevd} to orthogonal sketch matrix $\vect{Y}$. Step 2 and step 5 in Alg.~\ref{alg:ssrspSVD} is to generate two random matrices. Considering that $s_2$ and $s_3$ are small number, the time cost of step 3 is $O(nk^2+nvk+nvd)$ and step 6 cost $O(v(k^2+d)+v^2(k+d))$. Step 4 cost $O(nd^2)$, steps 7\textasciitilde8 cost $O(d^3)$ and step 9 cost $O(nd^2)$. Therefore, the time cost of Alg.~\ref{alg:ssrspSVD} is $O(nk^2+nvk+nvd+nd^2)$.

\subsection{The NetMF+ Algorithm}
In this section, we formally describe the NetMF+ algorithm as Algorithm~\ref{alg:NetMF+}, which consists of three steps: fast randomized eigen-decomposition, sparse-sign randomized single-pass SVD and spectral propagation.
In Alg.~\ref{alg:NetMF+}, step 1 compute fast randomized eigen-decomposition for modified Laplacian matrix $\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}$ and steps 2\textasciitilde3 form matrix $\vect{F},\vect{C}$ for Eq.~(\ref{eq4}). Step 4 in Alg.~\ref{alg:NetMF+} is to compute the SVD of $\hat{\vect{M}}^{\prime}$ by $\vect{F}$ and $\vect{C}$ through sparse-sign randomized single-pass SVD. Step 6 can be used to improve the embedding quality.
\begin{algorithm}[t]
	\caption{NetMF+}
	\label{alg:NetMF+}
	\small
	\SetAlgoLined
	\KwIn{A social network $G=(V,E,\vect{A})$; Normalized parameter $\alpha$; eigen-decomposition parameter $k,s_1,q$; Sparse-sign matrix parameter $s_2,s_3,z$; Spectral propagation step $p$; The dimension of embedding $d$; }
	\KwOut{An embedding matrix $\vect{E}\in \mathbb{R}^{n\times d}$}
	$\vect{D} = \mathrm{sum}(\vect{A},1), [\vect{U}_k,\vect{\Lambda}_k] = \mathrm{freigs}(\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha},k,s_1,q)$\;
	$\vect{K}=\vect{U}_k^{\top}\vect{D}^{-1+2\alpha}\vect{U}_k\vect{\Lambda}_k$\;
	$\vect{F}=\vect{D}^{-1+\alpha}\vect{U}_k,\vect{C}=\vect{\Lambda}_k\left(\sum_{r=1}^{T}\vect{K}^{r-1}\right)$\;
	$[\vect{T},\vect{\Sigma},\sim]=\mathrm{sparse\_sign\_rand\_single\_pass\_SVD}(\vect{F},\vect{C},d,s_2,s_3,z)$\;
	$\vect{E}=\vect{T}\sqrt{\vect{\Sigma}}$\;	$\vect{E}=\sum_{r=0}^{p}c_r\mathcal{L}^{r}\vect{E}$\tcp*{spectral propagation}
	\KwRet $\vect{E}$ as network embedding\;
\end{algorithm}

\textbf{Complexity Analysis.}
Now we analyze the time and space complexity of NetMF+. As for step 1, the input matrix $\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}$ is still an $n\times n$ matrix and have $2m$ nonzeros. According to section 4.1, it requires $O(q(mk+nk^2))$ time and $O(m+nk)$ space. As for steps 2\textasciitilde3, it requires $O(nk^2)$ time and $O(nk)$ space. As for step 4, the time cost is $O(n(k^2+d^2+vk+vd))$ and space cost is $O(n(k+d))$. Steps 5\textasciitilde6 require $O(pmd+nd)$ time and $O(m+nd)$ space. In total, the NetMF+ cost time $O(q(mk+nk^2)+nd^2+nvk+nvd+pmd)$ and the space cost is $O(m+nk)$.






\subsection{System Design}
Figure~\ref{img1} is the overview of the system design of NetMF+. 
\begin{figure}[h]
    \setlength{\abovecaptionskip}{0.08 cm}
	\centering
	\includegraphics[width=1\linewidth]{img1}
	\caption{The system design of NetMF+.}
	\label{img1}
\end{figure}
Each step of NetMF+ can be parallelized, and the storage overhead is mainly depends on the sparse adjacency matrix and the temporary matrix. Below we discuss the detail system design of each step. 

\textbf{Compression.} First, the system is built on the Graph Based Benchmark Suite (GBBS)~\cite{gbbs}, which is an extension of the Ligra~\cite{ligra} interface. The GBBS is easy to use and has already shown its practicality for many large scale fundamental graph problems. LightNE~\cite{lightne} introduces GBBS to network embedding problems and shows its superiority to real-world networks with hundreds of billions of edges. The main benefit of GBBS to NetMF+ is the data compression. In the beginning, a sparse adjacency matrix is always stored as CSR format, which is normally regarded as an excellent compressed graph representation~\cite{kepner2011graph}. However, CSR format is still a huge memory overhead for networks with hundreds of billions of edges. For example, storing a network with 1 billion vertices and 100 billion edges will cost 1121GB of memory which is unsatisfactory. Therefore, we need to compress it further and reduce memory cost. A compressed CSR format for graph from Ligra+~\cite{ligraplus}, which supports fast parallel graph encoding and decoding, is a good solution.

\textbf{Parallelization.} Second, the main time cost of the NetMF+ is Sparse Matrix-Matrix multiplication (SPMM) and Matrix-Matrix product (GEMM). In step 1 of NetMF+ system, it performs fast randomized eigen-decomposition on $\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}$. The random projection and the power iteration process in Alg.~\ref{alg:frevd} are, in essence, a product of an $n\times n$ sparse matrix and a $n\times (k+s_1)$ dense Gaussian random matrix, which is supported in MKL's Sparse BLAS Routines. However, MKL's Sparse BLAS Routine requires the sparse matrix in CSR format as the input, which contradicts the original intention of using GBBS. Fortunately, GBBS supports traversing all neighbors of a vertex $u$ for the compressed CSR format, and we can propose an SPMM operation with the help of GBBS. It can be computed as follows: First of all, We need to traverse n vertices in parallel, and then traverse neighbor $v$ of the vertex $u$ to obtain the element $\vect{D}(u,u)^{-\alpha}\vect{D}(v,v)^{-\alpha}$ corresponding to the sparse matrix. Finally, with the support of \texttt{cblas\_saxpy} in MKL, we multiply the element with the $v$-th row of the row-major matrix and add the result to the $u$-th row of the target matrix. The SPMM operation based on GBBS is slightly slower than MKL's SPMM operation, but ensuring memory-efficient. "eigSVD", "eig" and other operations are well supported by Intel MKL BLAS routines and Intel MKL LAPACK routines. Step 3 of NetMF+ system also involves applying repeated SPMM between an $n\times n$ sparse Laplacian matrix $\mathcal{L}$ and a $n\times d$ dense temporary matrix. MKL's Sparse BLAS Routines support the operation. Step 2 of the system involves matrix row sampling and batch GEMM operation. In detail, we need to compute $\vect{Y}(\vect{r},:)=\mathrm{trunc\_log}^{\circ}\left(\frac{\mathrm{vol}(G)}{b}\frac{1}{T}\vect{F}(\vect{r},:)\vect{C}\vect{F}(\vect{p},:)^{\top}\right)\vect{\Psi}$, where $\vect{r}$ is the batch row index. For each column of $\vect{\Psi}$, we can get $z$ nonzeros, and there are corresponding $z$ columns of $\vect{F}^{\top}$ to complete matrix multiplication based on element multiple and summation. This process is easy to implement through for loop under OpenMP~\cite{openmp}.

In conclusion, NetMF+ is implemented by C++. We use GBBS to reduce memory usage and implement a GBBS-based SPMM operation. For efficiency, we use the Intel MKL library for basic linear algebra operations, and we apply OpenMP for parallel programming especially in step 4 and step 7 of Alg.~\ref{alg:ssrspSVD}.


\section{Experiments}

\begin{table*}[t]
	\setlength{\abovecaptionskip}{0.05 cm}
	\caption{Statistics of datasets.}
	\label{tab:datasets}
	\small{
	\renewcommand\arraystretch{0.9}
	\begin{tabular}{r|rrrrr|rrrr}
		\toprule
		&\multicolumn{5}{c}{Multi-label Vertex Classification Task} & \multicolumn{4}{|c}{Link Prediction Task} \\
		\midrule
		& BlogCatalog& YouTube & Friendster-small & Friendster& OAG& Livejournal & ClueWeb & Hyperlink2014 & Hyperlink2012\\
		\midrule
		$|V|$ & 10,312& 1,138,499& 7,944,949& 65,608,376 & 67,768,244 &4,847,571 &978,408,098 &1,724,573,718 & 3,563,602,789\\
		$|E|$ & 333,983& 2,990,443& 447,219,610& 1,806,067,142 & 895,368,962 &68,993,773 &74,744,358,622 &124,141,874,032 & 225,840,663,232 \\
		\bottomrule
	\end{tabular}
	}
\end{table*}


\begin{figure*}[t]
	\setlength{\abovecaptionskip}{0.05 cm}
	\centering
	\includegraphics[width=\linewidth]{img2-eps-converted-to.pdf}
	\caption{Classification performance w.r.t. the ratio of training data. \small The x-axis represents the ratio of labeled data (\%), and the y-axis in the top and bottom rows denote the Micro-F1 and Macro-F1 scores (\%) respectively. For methods which fail to finish computation in one day or cannot handle the computation, their results are not available and thus not plotted in this figure.}
	\label{img2}
\end{figure*}



In this section, we evaluate the proposed NetMF+ method on the multi-label vertex classification task and link prediction task, which have been commonly used to evaluate previous network embedding techniques~\cite{prone,netsmf,deepwalk,lightne}. We introduce our datasets, experimental setup, and results in Section 5.1, Section 5.2, and Section 5.3, respectively. The ablation study is in Section 5.4.


\subsection{Datasets}
The statistics of datasets are listed in Table~\ref{tab:datasets}. 
We employ five datasets for the multi-label vertex classification task. BlogCatalog~\cite{blogdata} and YouTube~\cite{youtubedata} are small graphs with less than 10M edges, while the others are large graphs with more than 10M but less than 10B edges.
For the link prediction task, We have four datasets in which vertex labels are not available. Livejournal~\cite{livejournaldata} is the large graph, while the others are very large graphs with more than 10B edges. These datasets are of different scales but have been widely used in network embedding literature~\cite{lightne}. The two largest networks are hyperlinked networks extracted from the Common Crawl web corpus, containing billions of edges.

\subsection{Experimental Setup}
\begin{table}[t]
	\setlength{\abovecaptionskip}{0.05 cm}
	\setlength{\belowcaptionskip}{-0.25 cm}
	\caption{The most similar hardware configuration on Azure compared to the machine used in the original paper.}
	\label{tab:hardwareconfig}
	\centering
	\small{
	\renewcommand\arraystretch{0.9}
	\begin{tabular}{cccccc}
    \toprule
    Azure    & Method     & vCores        & RAM       & GPU     & Price(\$/h)             \\ \midrule
    NC24s v2   & GraphVite      & 24     & 448GB   & 4X P100    & 8.28     \\\midrule 
    E48 v3  & PBG      & 48      & 384GB    & 0    & 3.024      \\ \midrule
    \multirow{3}{*}{M128s} & \multirow{3}{*}{\begin{tabular}[c]{@{}l@{}}LightNE\\ ProNE\\ NetMF+\end{tabular}} & \multirow{3}{*}{128} & \multirow{3}{*}{2048GB} & \multirow{3}{*}{0} & \multirow{3}{*}{13.338} \\
                           &                                                                                   &                      &                         &                    &                         \\
                           &                                                                                   &                      &                         &                    &                         \\ \bottomrule
    \end{tabular}
    }
\end{table}
For NetMF+, all experiments are conducted on a server with two Intel\textsuperscript{\textregistered} Xeon\textsuperscript{\textregistered} E5-2699 v4 CPUs (88 virtual cores in total) and 1.5TB memory. In fact, LightNE~\cite{lightne} is essentially an enhanced version of NetSMF~\cite{netsmf} and is the most important method to compare. Thus we make a full comparison between NetMF+ and LightNE. We follow the experiment and evaluation procedures that were performed in LightNE. For PyTorch-BigGraph (PBG)~\cite{lerer2019pbg}, we compare NetMF+ with it on Livejournal which is an example dataset from the original paper~\cite{lerer2019pbg}. We compare NetMF+ with GraphVite~\cite{graphvite}, ProNE~\cite{prone} and LightNE in vertex classification task. We used the evaluation scripts and hyperparameters provided by the corresponding paper's GitHub Repos for comparison. For ProNE, we use the high-performance version released by LightNE GitHub repos. Considering that GraphVite and PBG are developed for different scenarios, we evaluate them and NetMF+ based on the cost of cloud server rental. We list the most suitable Azure Cloud server for each system in Table.~\ref{tab:hardwareconfig}.

\textbf{Vertex Classification Setting.} After generating the embedding, we randomly sample a portion of labeled vertices for training and use the remaining for testing. The task is completed by the one-vs-rest logistic regression model, which is implemented by LIBLINEAR~\cite{fan2008liblinear}. We repeat the prediction procedure five times and evaluate the average performance in terms of both Micro-F1 and Macro-F1 scores~\cite{tsoumakas2009mining}. Across all datasets, we set the embedding dimension $d$ to be 128 except on Friendster, where we set 96. We recommend that the oversampling parameter $s_1$ can usually be a small number of 10, 30, 50. We set parameter $b = 1, z = 8, k=256, s_1=30$ in most cases for NetMF+. We set $s_2 = 100$, and $s_3 = 1000$ according to~\cite{tropp2019streaming} which indicate that $s_3$ should more larger than $s_2$ for a better result. We observe that $T$, $\alpha$ and $q$ are more sensitive to the experimental results so that we use cross-validation to select them. 

\textbf{Link Prediction Setting.} We follow LightNE to set up link prediction evaluation for very large graphs --- we randomly exclude 0.00001\% edges from the training graph for evaluation. When training NetMF+ on the three very large graphs, we skip the spectral propagation step (due to memory issue) and set $d = 32$ except Hyperlink2012, where we use $d = 16$. We set $b=1, T = 5, z = 8, k=32, s_1=0$ and $s_2 = 0,s_3 = 5000$ for NetMF+. The ranking metrics on the test set are obtained by ranking positive edges among randomly sampled corrupted edges. We evaluate the link prediction task with metrics to be HITS@1, HITS@10, HITS@50 and AUC.

\subsection{Experiment Results}
\textbf{Comparison with PBG.} First, We follow exactly the same experiment setting in PBG~\cite{lerer2019pbg}. We set $T=5, \alpha=0.375, q=5, k=1024, d=1024, z = 8, s_1=50, s_2=100, s_3=500$ for NetMF+. Combining the price in Table~\ref{tab:hardwareconfig}, the results are reported in the Table~\ref{tab:ljcompare}.
\begin{table}[t]
    \setlength{\abovecaptionskip}{0.05 cm}
	\setlength{\belowcaptionskip}{-0.25 cm}
	\caption{Comparison with PBG on Livejournal.}
	\label{tab:ljcompare}
	\centering
	\small
	\renewcommand\arraystretch{0.9}
	\begin{tabular}{cccccc}
    \toprule
    Systems   & Time & Cost & MR & MRR & HITS@10 \\ \midrule
    PBG    &  7.3 hours &  \$22.1    & 4.25   & 0.87 &0.93         \\
    NetMF+ &  \textbf{9.7 mins}    & \textbf{\$2.2}     & \textbf{2.27}   & \textbf{0.90}     & \textbf{0.97}         \\ \bottomrule
    \end{tabular}
    \normalsize
\end{table}
It's obvious that NetMF+ has better performance for all metrics with less time. To be precise, NetMF+ is $45\times$ faster and $10\times$ cheaper than PBG.

\textbf{Vertex Classification Results.}
We summarize the multi-label vertex classification performance in Figure~\ref{img2}. To compare different algorithms' runtime and memory cost, we also list the efficiency comparison in Table~\ref{tab:runtimememory}. In BlogCatalog, NetMF+ achieves the best performance compared to other methods. It achieves significantly better Micro-F1 and Macro-F1 than the second place method LightNE (by 3.5\% on average), while only cost 2 seconds and 17 GB. In YouTube, NetMF+ show comparable Micro-F1 and Macro-F1 to LightNE while show better result than ProNE. In Friendster-small~\cite{friendsterdata}, NetMF+ shows slightly better Micro-F1 than the other methods while achieving comparable Macro-F1 results. In Friendster~\cite{friendsterdata}, NetMF+ is comparable to LightNE and ProNE. In OAG~\cite{oagdata}, the result of LightNE is produced when choosing large edge samples $M=17Tm$, and NetMF+ achieves better performance than LightNE (Micro-F1 improved by 4.2\% on average). Overall, NetMF+ has significantly better or comparable classification results compared to other methods with a lightweight memory and time cost. Except in OAG, the speed of NetMF+ is slower than ProNE, but both the Micro-F1 and Macro-F1 result of ProNE is inferior, which proves the effectiveness of the NetMF+ algorithm in learning challenging datasets. In general, the vertex classification results illustrate the effectiveness and lightweight of the NetMF+.
\begin{table}[t]
	\setlength{\abovecaptionskip}{0.05 cm}
	\setlength{\belowcaptionskip}{-0.25 cm}
	\caption{Efficiency comparison.}
	\label{tab:runtimememory}
	\centering
	\small{
	\renewcommand\arraystretch{0.9}
	\begin{threeparttable}
    \begin{tabular}{l|l|rrrr}
    \toprule
                                         Metrics  &  Datasets   & GVT\tnote{2} & ProNE\tnote{3} & LightNE & NetMF+ \\ \midrule
    \multirow{8}{*}{Time} & BlogCatalog      &  20.1~s   &79.0~s & 152.0~s & \textbf{2.0~s}        \\
                          & YouTube          &  93.0~s   &65.0~s & 96.0~s  & \textbf{32.0~s}         \\
                         & Friendster-small &  2.8~h    &5.3~m  & 7.5~m   & \textbf{4.9~m}         \\
                          & Friendster       &  20.3~h   &19.5~m & 37.6~m  & \textbf{17.7~m}        \\
                          & OAG              &  $-$      &\textbf{22.6~m}  & 1.5~h  & 1.2~h        \\
                          & ClueWeb          &  $-$      &$\times$ & 1.3~h & \textbf{38.9~m}        \\
                          & Hyperlink2014    &  $-$      &$\times$ & 1.8~h & \textbf{1.0~h}        \\
                         & Hyperlink2012    &  $-$      &$\times$ & 5.6~h & \textbf{1.2~h}        \\ \midrule
    \multirow{8}{*}{\makecell{Mem\tnote{1}\\(GB)}}& BlogCatalog      &  $\star$  & 18  & 273   & \textbf{17 }       \\
                        & YouTube          &  $\star$  & 28  & 83    & \textbf{27 }        \\
                       & Friendster-small &  $\star$  & 84  & 541   & \textbf{56 }        \\
                          & Friendster       &  $\star$  & 326   & 559  & \textbf{236 }         \\
                          & OAG              &  $\star$  &403  & 1391  & \textbf{384 }       \\
                         & ClueWeb          &  $\star$  &$\times$ & 1493  & \textbf{612 }       \\
                         & Hyperlink2014    &  $\star$  &$\times$ & 1500  & \textbf{1076 }         \\
                         & Hyperlink2012    & $\star$   &$\times$ & 1500  & \textbf{1323 }       \\ \bottomrule
    \end{tabular}
    \begin{tablenotes}
    \footnotesize
    \item[1] "Mem" represents the peak memory cost.
    \item[2] GVT is short for GraphVite. "$-$" indicates that the corresponding algorithm fails to complete within one day. "$\star$" indicates that we use rental cost to evaluate CPU-GPU hybrid system.
    \item[3] "$\times$" indicates that the corresponding algorithm is unable to handle the computation due to excessive space and memory consumption.
    \end{tablenotes}
    \end{threeparttable}
	}
	\vspace{-0.05in}
\end{table}

Next, we compare the rental costs of GraphVite and NetMF+ when applied to large graph data. Based on Table~\ref{tab:hardwareconfig} and Table~\ref{tab:runtimememory}, we get the cost of GraphVite on Friendster-small, and Friendster is \$23.2 and \$168.1, respectively, while the cost of NetMF+ is \$1.1 and \$3.9 respectively. Therefore, NetMF+ shows up to $69\times$ faster and $43\times$ cheaper compared to GraphVite.

\textbf{Link Prediction Results.}
 In addition, we need to compare the performance of different algorithms on link prediction task for very large graphs. Due to limited memory, ProNE cannot accomplish learning an embedding. GraphVite and PBG cannot finish the job in one day. Therefore, we only report the results of LightNE and NetMF+ in Table~\ref{tab:linkpred}. The results of LightNE are reported by choosing large edge samples for a best possible result in the 1.5 TB memory machine. Not only does NetMF+ achieve better performance regarding all metrics, but it also reduces time and memory costs significantly. The runtime of NetMF+ on ClueWeb~\cite{cluewebdata} is 38.9 minutes in Table~\ref{tab:runtimememory}, and the peak memory is 612GB, which is less than half of LightNE. It is worth mentioning that NetMF+ will take 27 mins and 1498 GB to learn on ClueWeb if MKL's Sparse BLAS Routine is used in step 1 of the system. The uncompressed ClueWeb graph (in CSR format) occupies 839 GB while the total peak memory cost of NetMF+ (using GBBS-based SPMM) on ClueWeb is 612GB, which proves the effects of compression of GBBS. The experimental results on Hyperlink2014~\cite{hyperlinkdata} also demonstrate the effectiveness and efficiency of NetMF+. For the largest network, Hyperlink2012~\cite{hyperlinkdata}, LightNE can only have small edge samples, making embedding ineffective in the link prediction task compared to NetMF+. It describes the fact that NetMF+ has better performance when it uses the same or similar amount of memory as LightNE. In conclusion, the results on the three very large graphs show that NetMF+ is significantly better than LightNE with less time and memory.
\begin{table}[t]
	\setlength{\abovecaptionskip}{0.05 cm}
	\setlength{\belowcaptionskip}{-0.25 cm}
	\caption{Link prediction comparison.}
	\label{tab:linkpred}
	\small
	\renewcommand\arraystretch{0.9}
    \begin{tabular}{l|l|llll}
    \toprule
                   Datasets     & Systems & HITS@1 & HITS@10 & HITS@50 & AUC \\ \midrule
    \multirow{2}{*}{ClueWeb}       & LightNE & 0.684  &  0.753  & 0.803   & 0.903\\
                                   & NetMF+  & \textbf{0.703} &  \textbf{0.772} &  \textbf{0.868}& \textbf{0.972}    \\ \midrule
    \multirow{2}{*}{Hyperlink2014} & LightNE & 0.364 & 0.5 & 0.628  & 0.874\\
                                   & NetMF+  & \textbf{0.503} & \textbf{0.646}&\textbf{0.739}&\textbf{0.913} \\ \midrule
    \multirow{2}{*}{Hyperlink2012} & LightNE & 0.073 & 0.189 & 0.348 & 0.751\\
                                   & NetMF+  &\textbf{0.615}  & \textbf{0.716} &\textbf{0.798}  &\textbf{0.928}  \\ 
    \bottomrule 
    \end{tabular}
	\normalsize
\end{table}

\subsection{Ablation Study} We first perform ablation studies on the effects of each step of the NetMF+ system and then focus on the parameter analysis.

\textbf{Efficiency and Effectiveness of Each Step of NetMF+.} First, we need to focus on the fast randomized eigen-decomposition (the step~1 of NetMF+ system). Halko et al.~\cite{halko2011finding} has shown that it's challenging to perform eigen-decomposition on $\vect{D}^{-0.5}\vect{A}\vect{D}^{-0.5}$. Although Halko et al. stated that this problem could be improved by using power iteration scheme, it is still impractical for large real-world networks. We choose YouTube as an example, and its related eigenvalues are shown in Figure~\ref{eigenvaluesimg}. 
The eigenvalues of $\vect{D}^{-0.5}\vect{A}\vect{D}^{-0.5}$ with $q=20$ and $q=100$ are far from the correct result, which the largest eigenvalue should be 1. The \texttt{eigsh}~\cite{arpack1998} is not able to complete this job in one day. The results illustrate the requirement to use modified Laplacian matrix as input for fast randomized eigen-decomposition.
When we choose $\alpha = 0.4$, the eigenvalues of fast randomized eigen-decomposition ($q=20$) on $\vect{D}^{-0.4}\vect{A}\vect{D}^{-0.4}$ is indistinguishable from the eigenvalues computed by \texttt{eigsh} when $k<140$. The runtime of fast randomized eigen-decomposition is 20 seconds while \texttt{eigsh} costs 31 minutes, which proves the accuracy and efficiency of fast randomized eigen-decomposition. 
\begin{figure}[t]
	\centering
	\includegraphics[width=0.625\linewidth]{img3-eps-converted-to.pdf}
	\vspace{-0.2in}
	\caption{The computed eigenvalues on YouTube data with different $\alpha$.}
	\vspace{-0.05in}
	\label{eigenvaluesimg}
\end{figure}

Next, we focus on the effects of sparse-sign randomized single-pass SVD and spectral propagation.  NetMF+ (w/o spectral) and NetMF+ (w/ spectral) represent embeddings from step~2 and step~3 of the NetMF+ system, respectively, and we need to compare them. The vertex classification results are shown in Figure~\ref{comparespectral}.
\begin{figure}[t]
	\setlength{\abovecaptionskip}{0.05 cm}
	\setlength{\belowcaptionskip}{-0.1 cm}
	\centering
	\includegraphics[width=.9\linewidth]{img4-eps-converted-to.pdf}
	\caption{The embedding performance comparison.}
	\label{comparespectral}
\end{figure}
NetMF+ (w/o spectral) results are satisfactory, which proves the effectiveness of the sparse-sign randomized single-pass SVD. Combining with step~3 of the system, NetMF+ (w/ spectral) show better results, especially the significant improvement on Macro-F1, which describes the effects of spectral propagation.

\textbf{Parameter Analysis}
The peak memory cost is constant when we fix the oversampling parameters $s_1,s_2,s_3$, column sparsity $z$, eigen-decomposition rank $k$, and embedding dimension $d$. The parameters that strongly impact the quality of learned embedding are $\alpha$, $T$ and $q$. $\alpha$ determines the how slow decaying speed of the eigenvalues of the modified Laplacian and thus determines the accuracy of the fast randomized eigen-decomposition and the approximation error of $\vect{M}$. In general, $\alpha\in[0.35,0.5)$ and $T=5$ or $10$ will always be an appropriate choices. We need to pay attention to $q$, which determines the accuracy of the fast randomized eigen-decomposition. Thus, we make a tradeoff between the quality of the learned embedding and the overall system runtime.
We choose OAG dataset for test. When we fix $\alpha=0.45, T=10, k=256,s_1=20,s_2=300$ and $s_3=5000$ we enumerate the $q$ from $\{10,15,20,25,30\}$, the efficiency-effectiveness trade-off result of OAG is shown as Figure~\ref{oagee}. 
We also add LightNE results with different edge samples $M$ from $\{1Tm,5Tm,7Tm,10Tm,13Tm,17Tm\}$. The peak memory of NetMF+ is 384GB while that of LightNE is 553 GB, 682 GB, 776 GB, 936 GB, 1118 GB and 1391 GB, respectively. Figure~\ref{oagee} shows that NetMF+ works better than LightNE when the runtime is fixed. The experiments prove that users can adjust NetMF+ flexibly according to time/memory budgets and performance requirements.
\begin{figure}[t]
	\setlength{\abovecaptionskip}{0.05 cm}
	\setlength{\belowcaptionskip}{-0.1 cm}
	\centering
	\includegraphics[width=.9\linewidth]{img5-eps-converted-to.pdf}
	\caption{The trade-offs between efficiency and performance.}
	\label{oagee}
\end{figure}
\section{Conclusion}
In this work, we propose NetMF+ employing the techniques of fast randomized matrix factorization,  sparse-sign randomized single-pass SVD and GBBS, which produce a fast and memory-efficient network embedding system. Compared with state-of-the-art network embedding systems, NetMF+ achieves the best performance on nine benchmarks. With the help of fast randomized eigen-decomposition and a novel sparse-sign randomized single-pass SVD, NetMF+ consumes much less computational time and memory
while generating a promising embedding.
With GBBS and GBBS-based SPMM, the memory cost is further reduced. 
The main computations of NetMF+ are highly parallelizable, which is thus well supported by the MKL library and OpenMP. With these techniques, NetMF+ can learn high-quality embeddings from networks with more than 100 billion edges in about an hour with  1 TB or less memory cost. 

In the future, we will extend the method to handle dynamic network embedding problems.
\section{Appendix}
The following lemmas will be useful in our proof.
\setcounter{lemma}{0}
\renewcommand{\thelemma}{\textbf{\arabic{lemma}}}
\setcounter{theorem}{0}
\renewcommand{\thelemma}{\textbf{\arabic{lemma}}}
\setcounter{equation}{0}
\renewcommand\theequation{\arabic{equation}}

\begin{lemma}
	\label{lemma1}
	(\cite{trefethen1997numerical}) Singular values of a real symmetric matrix are the absolute values of its eigenvalues.
\end{lemma}

\begin{lemma}	
	\label{lemma2}
	(\cite{halko2011finding,musco2015nips}) For input matrix $X$, the Alg.~\ref{alg:basic_rsvd} with the PI scheme has the following guarantee: 
	\begin{equation*}
	    \Vert\vect{X}-\vect{Q}\vect{Q}^{\top}\vect{X}\Vert_F=\Vert\vect{X}-\vect{U}\vect{\Sigma}\vect{V}^{\top}\Vert_F\leq (1+\varepsilon)\Vert\vect{X}-\vect{X}_k\Vert_F,
	\end{equation*}
	with high probability. Here, $\vect{X}_k$ is the best rank-$k$ approximation of $\vect{X}$.
\end{lemma}

\begin{lemma}	
	\label{lemma3}
	(\cite{horn1994topics}) Let $\vect{B},\vect{C}$ be two $n \times n$ symmetric matrices. Then for the decreasingly ordered singular values $\sigma$ of $\vect{B},\vect{C}$ and $\vect{BC}$, $\sigma_{i+j-1}(\vect{BC})\leq \sigma_i(\vect{B})\times \sigma_j(\vect{C})$ holds for any $1\leq i,j\leq n$ and $i+j\leq n+1$.
\end{lemma}


\begin{theorem}
	Let $\Vert\cdot\Vert_F$ be the matrix Frobenius norm. Then \\
	$\Vert\vect{D}^{-\frac{1}{2}}\vect{A}\vect{D}^{-\frac{1}{2}}-\vect{D}^{-\frac{1}{2}+\alpha}\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top}\vect{D}^{-\frac{1}{2}+\alpha}\Vert_F\leq(1+\varepsilon)\sqrt{\sum_{j=k+1}^{n}|\lambda_j|^2}.$
	with a high probability. Here $|\lambda_j|$ is the $j$-th largest absolute value of eigenvalue of $\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}$. 
\end{theorem}
\begin{proof}
	Apply Lemma~\ref{lemma1} and Lemma~\ref{lemma2}, we will have the approximation for Alg.~\ref{alg:frevd} as
	\begin{equation}
	    \label{appro_eq}
	    \begin{split}
	         \Vert\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}-\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top}\Vert_F
	         &\leq (1+\varepsilon) \Vert\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}-(\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha})_k\Vert_F\\
	         &=(1+\varepsilon)\sqrt{\sum_{j=k+1}^{n}|\lambda_j|^2},
	    \end{split}
	\end{equation}
	where $|\lambda_j|$ is the $j$-th largest absolute value of eigenvalue of $\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}$.
	
	Consider the approximation in Thm.~1, we have
	\begin{equation}
		\label{proof_eq1}
		\begin{split}
			&\left\Vert\vect{D}^{-\frac{1}{2}}\vect{A}\vect{D}^{-\frac{1}{2}}-\vect{D}^{-\frac{1}{2}+\alpha}\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top}\vect{D}^{-\frac{1}{2}+\alpha}\right\Vert_F\\
			&= \left\Vert\vect{D}^{-\frac{1}{2}+\alpha}\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}\vect{D}^{-\frac{1}{2}+\alpha} - \vect{D}^{-\frac{1}{2}+\alpha}\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top}\vect{D}^{-\frac{1}{2}+\alpha} \right\Vert_F\\
			&=\sqrt{\sum_{i=1}^{n}\sigma_i\left(\vect{D}^{-\frac{1}{2}+\alpha}\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}\vect{D}^{-\frac{1}{2}+\alpha}-\vect{D}^{-\frac{1}{2}+\alpha}\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top}\vect{D}^{-\frac{1}{2}+\alpha}\right)^2}.
		\end{split}
	\end{equation}
	Here, we have $\alpha\in(0,0.5]$ and we use $d_i$ to denote $\vect{D}(i,i)$. Since every $d_i$ is positive, the decreasingly ordered singular values of the matrix $\vect{D}^{-\frac{1}{2}+\alpha}$ can be constructed by sorting $d_i^{-\frac{1}{2}+\alpha}$ in the non-increasing order. In particular, we mark $d_{min}^{-\frac{1}{2}+\alpha}$ is largest singular values where $d_{min}$ is the smallest vertex degree.
	Apply Lemma~\ref{lemma3} twice for Eq.~(\ref{proof_eq1}), we will have
    \begin{equation}
        \label{proof_eq2}
        \begin{split}
            &\Vert\vect{D}^{-\frac{1}{2}}\vect{A}\vect{D}^{-\frac{1}{2}}-\vect{D}^{-\frac{1}{2}+\alpha}\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top}\vect{D}^{-\frac{1}{2}+\alpha}\Vert_F\\
            &\leq \sqrt{\sum_{i=1}^{n}\sigma_{1}(\vect{D}^{-\frac{1}{2}+\alpha})^2\sigma_{i}(\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}-\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top})^2 \sigma_{1}(\vect{D}^{-\frac{1}{2}+\alpha})^2}\\
            &=\sigma_{1}(\vect{D}^{-\frac{1}{2}+\alpha})^2\sqrt{\sigma_{i}(\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}-\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top})^2}\\
            &=d_{min}^{-1+2\alpha}\sqrt{\sigma_{i}(\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}-\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top})^2}\\
            &\leq\sqrt{\sigma_{i}(\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}-\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top})^2}.
        \end{split}
    \end{equation}
	Combining Eq.~(\ref{appro_eq}) with Eq.~(\ref{proof_eq2}), we will have
	\begin{equation}
        \label{proof_eq3}
        \begin{split}
            &\Vert\vect{D}^{-\frac{1}{2}}\vect{A}\vect{D}^{-\frac{1}{2}}-\vect{D}^{-\frac{1}{2}+\alpha}\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top}\vect{D}^{-\frac{1}{2}+\alpha}\Vert_F\\
            &\leq\sqrt{\sigma_{i}(\vect{D}^{-\alpha}\vect{A}\vect{D}^{-\alpha}-\vect{U}_k\vect{\Lambda}_k\vect{U}_k^{\top})^2}\\
            &\leq(1+\varepsilon)\sqrt{\sum_{j=k+1}^{n}|\lambda_j|^2}.
        \end{split}
    \end{equation}
	
\end{proof}


