\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}


\section{System description}\label{sec: System Description}
\section{Neural network model for fiber nonlinearity compensation}

We investigate a 10,080 km submarine optical transmission system carrying a single channel Gbaud polarization-multiplexed (PM)-16 quadrature amplitude modulation (QAM) signal. While recent demonstrations have shown the transmission capacities have approached the Shannon limit in the linear regime of optical fibers, the nonlinear impairment remains the major limiting effect in long-distance transmission systems~\cite{essiambre2010capacity}. The evolution of the signals in the fiber link is described by the Nonlinear Schr\"odinger Equation:
\begin{equation}\label{eq: nlse}
\frac{\partial{u_{x/y}(t,z)}}{\partial{z}} + \underbrace{i\frac{\beta_{2}}{2}\frac{\partial^2{u_{x/y}(t,z)}}{\partial^2{t}}}_\text{\footnotesize dispersion}=\underbrace{i\frac{8}{9}\gamma[\mid u_{x/y}(t,z)\mid^2 + \mid u_{y/x}(t,z)\mid^2 ]\mid u_{x/y}(t,z)}_\text{\footnotesize Kerr nonlinearity},
\end{equation}
where $u_{x/y}(t,z)$ is the optical field at the x and y polarizations, respectively, $\beta_2$ is the group velocity dispersion governing the linear impairment, and $\gamma$ is the Kerr nonlinear coefficient governing the nonlinear impairment in the fiber transmission systems. Although state-of-the-art DSP can overcome linear impairments, nonlinear impairments remain difficult to overcome due to the difficulty in applying computationally expensive nonlinear compensation algorithms using ASICs. The nonlinearity impairment can be treated as the the perturbation due to the nonlinear effects, \emph{i.e.}, $u_{x/y}(t,z) = u_{0, x/y}(t, z) + \Delta u_{x/y}(t,z)$, where $u_{0, x/y}(t, z)$ is the solution of linear propagation and $\Delta u_{x/y}(t,z)$ is the nonlinear perturbation. The idea of NLC is to build up a model of the nonlinear perturbation $\Delta u_{x/y}(t,z)$ with the abundance of transmission data using a NN, and remove the distortion by subtracting the perturbation from the received optical field $u_{x/y}(t,z)$.

Here we first consider how fiber nonlinearities impair the signal on the x-polarization as an example. We define a term \textit{triplet} as $T = H_{n}H_{m+n}^* H_{m} + V_{n}V_{m+n}^* H_{m}$, where $H$ and $V$ are the received symbol sequences at the x- and y- polarizations, respectively. $m$ and $n$ are symbol indices with respect to the symbols of interest $H_0$. Triplet $T$ represents the nonlinearities added on the X-polarization signals, which is caused by the nonlinear interactions among the signals at the same polarization (i.e., $H_{n}H_{m+n}^* H_{m}$) and orthogonal polarizations ($V_{n}V_{m+n}^* H_{m}$), as shown in Eq.~\ref{eq: triplet}. Triplet $T$ corresponds to the nonlinear terms on the right side of Eq.~\ref{eq: nlse}. Then, the nonlinear perturbation can be approximated by a combination of triplets $T$~\cite{tao2011multiplier}:
\begin{equation}\label{eq: triplet}
\begin{split}
\Delta u_{x/y}(0,z) = \sum P_0^{3/2}(\underbrace{H_{n}H_{m+n}^* H_{m}}_{\begin{subarray}{c}\text{\footnotesize intra-polarization}\\ \text{\footnotesize nonlinearity}\end{subarray}}
+ \underbrace{V_{n}V_{m+n}^* H_{m})}_{\begin{subarray}{c}\text{\footnotesize inter-polarization}\\ \text{\footnotesize nonlinearity}\end{subarray}}     C_{m,n}
\end{split}
\end{equation}
where $P_0$ is the launch power, $C_{m,n}$ are the nonlinear perturbation coefficients.

\begin{figure}
\includegraphics[width=1\linewidth]{figure_new/system_and_chip_2}
\spacing{1.}
\vspace*{-10mm}
\caption{\label{fig:system} (a)~(left) Neural network (NN) model for fiber nonlinearity compensation (NLC). The NN has an input layer with nonlinear triplets, 2 hidden layers with of 2 and 8 neurons in each layer, and an output layer with 2 output nodes corresponding to the real and imaginary parts of the estimated nonlinearity on one polarization, and (right) PNN implementation of the second hidden layer of the NN model based on weight-and-broadcast architecture\cite{tait2014broadcast,tait2017neuromorphic}. Red (blue) arcs are n-doped (p-doped) regions. (b) A false-color confocal micrograph of the PNN device. (c) Chip packaging and optical coupling setup. (d) A detailed circuit diagram of a photonic neuron.}
\end{figure}

Fig.~\ref{fig:system}(a) shows the NN model for NLC. The NN is a fully connected feed-forward network consisting of an input layer with 892 triplets, two hidden layers with 2 and 8 neurons, respectively, and an output layer. The model is trained using the signal on either x or y polarization only, but the trained parameters can be applied to both polarizations. Here, we choose polarization x as an example. The inputs of the NN are the triplets $T$ for the x-polarization signal~\cite{zhang2019field}.
The two outputs of the NN are the real and imaginary components of the estimated nonlinearity $\hat{H}_{k, NL}$ at the k-th symbol. The recovered symbol $\hat{H}_{k}$ is obtained by subtracting $\hat{H}_{k,~NL}$ from the received symbol $H_{k}$, i.e. $\hat{H}_{k}=H_{k}-\hat{H}_{k,~NL}$. The NN is trained by searching for the tensor parameters that minimize the mean squared error (MSE) between the $\hat{H}_{k}$ and the transmitted symbol $A_{x,k}$. After training, the NN can execute the inference task of estimating the optical nonlinearity of every incoming symbol. The NN is much simpler than those used for many other applications (e.g. image processing). Nevertheless, high-speed communications signals require fast operating speed, posing many challenges to today's digital electronic hardware. For example, to process 32 Gbaud PM-16QAM signals, NN hardware must have a throughput of approximately 232.4 TOPS. A typical NN accelerator, Tensor Processing Unit (TPU), offers a peak operation throughput of 92 TOPS\cite{jouppi2017datacenter}. As a result, the need for evaluating NN models for optical communication signals necessitates the development of new hardware.

\section{Photonic hardware implementation}
The NN model are implemented using a PNN with WDM-based broadcast-and-weight architecture, as shown in the yellow region of Fig.~\ref{fig:system}(a). The information of each neuron is encoded on a unique wavelength of light and multiplexed on a single optical waveguide.
In the feed-forward network, the neuron output in the previous layer is first broadcast to all neurons in the next layer by a power splitter, then weighted in parallel with an array of microring resonators (MRR) in an MRR weight bank configuration, and finally summed by a balanced photodetector.  Weights are applied by the MRR banks by tuning the partial transmission of signals. The weighted signals are detected, generating photocurrent to drive the optical modulators in the second layer. This process involves a nonlinear conversion of optical power into electrical current and back into the signal pathway, and hence serves as the nonlinear activation function (i.e. neuron) in the PNN\cite{tait2019silicon}.

The signal pathways in the PNN can be fabricated using standard silicon photonic platforms, which provides critical aspects of feasibility and economies of scaling to the full NN model for NLC. Several notable studies have demonstrated the integration of thousands of photonic components on a single silicon photonic chip\cite{sun2013large,poulton20208192}.
For the first time, we show that the PNN can provide accurate and fast weighting, summation, and biased nonlinearity for high-speed signals, and we experimentally implement the second hidden layer of the NN model using a 4$\times$2 PNN, as shown in Fig.~\ref{fig:system}(b). In particular, integrating fast photonic neurons on chip is a critical move toward achieving high-speed NNs for real-time applications such as NLC. In the Discussion section, we discuss that the PNN can be robustly scaled to incorporate the full NN model.

The PNN device comprises two arrays of MRR weight banks connected to two photonic neurons. The second hidden layer consists of 8 neurons, and each neuron is connected to the two neurons in the first hidden layer by a 2$\times$2 weight matrix. To emulate the second hidden layer with the 4$\times$2 PNN circuit, we feed the two first-layer outputs to the PNN four times with four subsets of weights and biases.
The photonic neuron circuit diagram is shown in Fig.~\ref{fig:system}(d). The MRR weight bank is implemented with in-ring N-doped photoconductive heaters\cite{tait2018feedback} and provides the key functionality to configure connection strengths (\emph{i.e.} weights or synapses). The weight is determined by the MRR transmission, which can be thermally tuned by adjusting the electrical current applied to the N-doped heaters\cite{tait2018feedback}. The N-doped heaters allow continuous, multi-channel control of the MRR weight bank with accuracies up to 8 bits\cite{huang2020demonstration,zhang2021microring}, which quantifies the precision of weights implementation affected by the noise in the setup. The bit resolution obtained in photonic matrix multiplier is comparable that used in DSP ASICs. The MRR weight bank has two complementary optical outputs, each of which is detected by a germanium-on-silicon photodetector.
The two photodetectors form a balanced photodetector~\cite{hai201316}, where the output photocurrent represents the subtraction operation between the two MRR weight bank outputs, resulting in a complementary –1 to +1 continuous weight range. The photocurrent, combined with a forward bias current $I_{bias}$, modulates the transmission of the MRR modulator (\emph{i.e.} photonic neuron) via free-carrier injection to the p-n junction, and hence modulates the optical power of a continuous-wave laser (labeled as "neuron pump"). The on-chip capacitor and resistor provide a network matching circuit for efficient optical-electrical-optical (OEO) conversion (see Supplementary information \textit{Photonic neuron simulation} for detailed circuit analysis). The MRR modulator exhibits nonlinear electrical-to-optical transfer functions, which produce the activation function in the NN. Furthermore, the neuron biases can be configured by adjusting the electrical current applied to the MRR modulator.

\begin{figure}
\includegraphics[width=1\linewidth]{figure_new/training_final}
\spacing{1.08}
\caption{\label{fig:training} (a) Flowchart of training process. (b) Activation function of photonic neuron. (c) The Q factor converges as the NN is trained, 32,016 symbols generated from the transmission system to monitor the training process; red trace: Q factor before nonlinearity compensation. (d) Transmission performance comparison. The neural network model for NLC is trained with the activation function of photonic neuron. The 3-step and 6-step DBP mean that fiber nonlinearity is compensated based on the conventional DSP approach digital backpropagation by equally dividing 10,080 km transmission link into 3 and 6 sections, respectively.}
\end{figure}

\section{Training with photonic neural network}
The proposed PNN provides trainable parameters including both weights and biases. Trainable biases are critical for the NLC application, because they can optimize the transfer function shapes to best approximate the higher-order nonlinearity features from the symbols. To ensure that the trained NN parameters can be accurately mapped on the photonic device, we include the activation function of photonic neurons in the training phase. As shown in Fig.~\ref{fig:training}(a), we begin by characterizing the on-chip activation function by sending a training data sequence into the weight bank and recording the corresponding neuron output waveform. The transfer function is fitted to a Lorentzian function, as shown in Fig.~\ref{fig:training}(b), since the E/O response of MRR modulators is typically a Lorentzian function. The training data have the same data rate as the test data, in order to capture the bandwidth and the noise in the device that would affect the activation function shape. Following that, we apply the characterized activation function in the NN model and train the neural network with the Adam learning algorithm\cite{kingma2014adam}.

The NN is trained with 32,106 training symbols in a single channel 32 Gbaud PM-16-QAM signal transmitted over a 10,080 km pure silica core fiber link. The weights are constrained to a range of -1 to 1 during the training, according to the operation range of the MRR weight bank. We monitor the Q-factor of a cross-validation (CV) data set to ensure that the NN is successfully trained to recognize the fiber nonlinear distortion. The CV data set contains 32,106 symbols and is independent with the training data set. As shown in Fig.~\ref{fig:training}(c), the Q-factor of the CV data increases gradually in the training process and finally converges to 8.1 dB after 21,600 steps. To enable efficient training, here we train the NN with a launch power higher than the optimal value. After training, a test data set with varying launch powers is sent to the optical link and processed by the trained NN. The nonlinear perturbation obtained from the NN is weighted by the difference in the launch power between the training and test set, and is then subtracted from the received signal. The Q factor of the signal after NLC is calculated and plotted in Fig.~\ref{fig:training}(d). The optical power at the Q factor peak is the optimal power that should be launched in the transmission fiber. Using the simulated PNN, the Q factor of the transmission link is increased by 0.66 dB as compared to that without NLC. Compared with 3-step DBP, the PNN is capable of achieving comparable Q factor performance.

\begin{figure}
\includegraphics[width=1\linewidth]{figure_new/experimental_setup}
\spacing{1.08}
\caption{\label{fig:setup} Experimental Setup. MZ-Mod: Mach-Zehnder modulator; AWG: arbitrary waveform generator, CPL: coupler; EDFA: erbium-doped fiber amplifier, PC: polarization controller; SMU: source meter unites; PNN: photonic neural network; PD: photodetector.
}
\end{figure}

\section{Photonic neural network implementation}
Following training, we implement the trained NN model to the PNN device, and evaluate whether the PNN can accurately compensate fiber nonlinearity impairment for high-speed signals. (See Supplementary Table.~3 for the weights and biases values). The experimental setup is shown in Fig.~\ref{fig:setup}. To emulate the first hidden layer outputs, two external cavity lasers are modulated with two waveforms generated from a two-channel arbitrary waveform generator (AWG). Each channel operates at a sampling rate of 375 MSample/s and generates a waveform of 32,106 test symbols at a symbol rate of 46.875 Mbaud/s. The test signal is independent of the signals in the training phase. The wavelengths of two lasers are initially aligned with the resonances of two MRRs in the MRR weight banks. The two signals are then combined and coupled into the PNN chip via a grating coupler. Two neuron pumps are derived from another two external cavity lasers with wavelengths close to the resonances of the MRR modulators. The modulator’s p-n junction is forward biased, allowing photocurrents from the balanced photodetector to modulate the neuron pumps through free-carrier injection. Under the carrier injection mode, the measured neuron bandwidth is 150 MHz. In Supplementary information Fig.~2, we show that, with the same photonic circuits and optimized MRR modulator geometry, the neuron bandwidth can be practically improved to 10 GHz.

\begin{figure}
\includegraphics[width=1\linewidth]{figure_new/results_final}
\spacing{1.08}
\caption{\label{fig:results_2}(a) waveforms at the second hidden layer output. Green dashed line: photonic neuron waveform; red dot: photonic neuron waveform downsampled to the data rate of 16-QAM signals; blue triangle: simulated waveform after downsampling. Ideally, the red points should overlap with the blue. And biased activation function characterized from the measured waveforms (blue dot) versus ideal activation function. (b) constellations of the 16-QAM signal on the X-polarization. Red: without nonlinearity compensation (NLC); blue: NLC using PNN; green: NLC using simulated NN. Red cross: constellation diagram of 16QAM without noise.}

\end{figure}

The trained NN parameters are grouped into four subsets, each with a 2$\times$2 weight matrix and 2 neuron biases.  Each weight matrix is mapped to a current matrix $\bm{I_h}$ and applied to the in-ring heaters of the MRR weight banks using the control method described in Ref~\cite{huang2020demonstration}. To apply neuron biases, the modulator bias current $I_{bias}$ is tuned to approximate the trained transfer function. The 2$\times$2 current matrix and 2 modulator biases are simultaneously applied to the PNN chip. Fig.~\ref{fig:results_2}(a) shows the waveform snapshots from the PNN output, which presents the outputs of the second hidden layer, together with the measured transfer function of the eight neurons. Each neuron's transfer function corresponds to the different piece of the characterized Lorentzian function under intended neuron biases. We simulate the same NN with a 32-bit GPU-assisted workstation to benchmark the accuracy of PNN implementation. Each neuron is evaluated by comparing the experimental PNN output (red dots) with the simulated NN (blue dots). The mean squared errors (MSEs) of photonic neurons are calculated over 32,016 test symbols. MSEs are used to quantify the errors caused by noises and derivations from target weights and activation functions. The MSEs vary between 0.008 to 0.014. The signal's Q factor is then calculated from the eight photonic neuron outputs (see Supplementary information \textit{Q factor calculation}). Without nonlinearity compensation, the optimal Q-factor of the signal is 9.34 dB, obtained at the launch power of -1 dBm. With the PNN, the signal's Q factor is improved by 0.60 dB, and the constellation diagram after NLC is plotted in Fig.~\ref{fig:results_2}(b.i), showing less noise in contrast to that without NLC. The Q-factor obtained from the simulated NN indicates a 0.65 dB gain in Q-factor. That is, the penalty of loading the neural network to the PNN is only 0.05 dB after taking into account all inherent noise from the physical components and equipment.

\section{Scalability for NLC}
Mainstream silicon photonic platforms offer an opportunity to implement the entire neural network for NLC on a single silicon photonic chip. Silicon photonic circuits with a few thousands of elements have been demonstrated~\cite{sun2013large,poulton20208192}. Weight pruning\cite{zhang2019field} is an effective pathway for reducing the PNN footprint and power consumption, which can be accomplished by having a few more neurons in hidden layers (see Fig.~\ref{fig:activation_function_comparison}) while reducing the number of input triplets. The current NN model has 892 input triplets. After weight pruning, the number of inputs can be reduced by three times to 300. The photonic neural network, in this case, has 952 optoelectronic components. These components require 0.38 mm\textsuperscript{2} chip space. Electrical traces and pads connecting to the optoelectronic components provide electrical biases to the photonic circuit.
Electrical traces and pads connecting to the optoelectronic components provide electrical biases to the photonic circuit. The amount of chip space used by electrical traces and pads varies significantly depending on the chip floorplan and packaging technologies. For example, the electrical traces can be very short using the flip-chip bonding technology, and the pads take up 0.76 mm\textsuperscript{2} areas (each bonding pad takes 20$\mu$m $\times$ 20$\mu$m area)\cite{carroll2016photonic}. Ref~\cite{sun2013large} have demonstrated that silicon photonic devices with comparable footprints are feasible.

A challenge of large-scale silicon photonic circuits is device non-uniformity, i.e., the dimensions and performances of the devices across the chip are not uniform and deviate from the design. In PNNs, device non-uniformity will cause variations of the relationship between the actuated MRR currents and the actual weights, the activation function shape, etc. The influences of device variations can be minimized by pre-calibrating each device or using control methods, such as the feedback MRR control~\cite{tait2018feedback,huang2020demonstration}. Device variances can be further reduced through post-fabrication technologies such as trimming~\cite{atabaki2013accurate} or in-situ training\cite{wang2019situ}. Another challenge deals with light sources. Frequency combs-based WDM sources\cite{brasch2016photonic} allow the generation of broadband and evenly spaced emission wavelengths from a single laser, and thus offer a practical approach of generating WDM light sources for the PNN. Soliton-based frequency combs built on CMOS-compatible photonic integrated circuits can generate hundreds of comb lines. Such systems can integrate with the PNN device.

\section{Claiming cascadability for NLC}

We discuss the feasibility of scaling the current chip to a multi-layer neural network. This issue refers to the need for cascadability, that is, the ability to have the one neuron output to excite the neurons in the next layer. Cascadability requires 1) the neuron output can evoke an equivalent response in neurons in the subsquent layer; 2) to avoid noise propagation along with multiple neural layers. The main conclusion is, in order to have one layer driving the next layer, enough power must be provided to compensate for loss and power splitting due to fan-out to the next layer. Such power can come from optical lasers or electrical gain via a high transimpedance. The electrical gain can be provided with a passive resistor or an active TIA, both of which are compatible with today’s silicon photonic platform. For example, we have implemented an on-chip resistor in our current PNN design. Alternatively, we can use an efficient modulator (i.e. small $V_{\pi}$ and junction capacitor C) to reduce the power requirement. We derive the devices and power requirement of maintaining the cascadability in the Supplementary information \textit{
Performance comparison between DSP and PNNs}.

Regarding noise propagation, we found that if the voltage swing to the modulator is large enough, the nonlinear transfer function of the modulator can suppress the noise and avoid noise propagation across the network\cite{de2019noise}. Such requirements can also be addressed with sufficient power or with an efficient modulator, as argued in~\cite{nozaki2019femtofarad,de2019noise,nahmias2019photonic}.

\section{Impact of activation function and the number of neurons}

In the NLC application, the NN is trained to recognize the fiber nonlinearity perturbation from the received symbols. An appropriate choice of activation function can benefit such regression problems by including higher-order nonlinearities among the symbols. Fiber nonlinearity perturbation manifests itself as long, heavy-tailed symbols due to nonlinear phase noise. To extract such features, especially the tail edge, an activation function with an unbounded range is preferred. Fig.~\ref{fig:activation_function_comparison}(a) shows the performance gap between the Lorentzian function of an MRR modulator neuron and the Leaky ReLU function which has no upper boundary for positive input. Generally, for different NN applications, activation functions need to be chosen to synthesize particular tasks. In contrast to a single cavity MRR, photonic devices with multiple coupled cavities can exhibit different lineshapes beyond the Lorentzian shape. Additional transfer function programmability could be achieved at the expense of additional circuits and calibration complexity using multiple-cavity modulators for some or all of the neurons\cite{huang2020chip}. Reconfigurable elements such as heaters or phase change materials allow optimizing device's transfer functions for different machine learning tasks.

In addition, as shown in Fig.~\ref{fig:activation_function_comparison}(b), the Q-factor improves as the number of neurons on the second-hidden layer increases. This can be expected because, with more neurons, the network can more precisely constitute different parts of nonlinearity perturbation over symbols. PNNs add almost no delays as the NN size increases, because the major operations in a PNN are computed at a single time step.

\begin{figure}[ht]
\center\includegraphics[width=0.5\linewidth]{figure_new/discussion_final}
\caption{\label{fig:activation_function_comparison}(a)~Leaky ReLU, optimized for fiber nonlinearity compensation as demonstrated by Zhang et al.\cite{zhang2019field}, and Lorentzian, characterized from the PNN. (b)~The comparison of Q factor improvement between Leaky ReLU (red) and Lorentzian (blue) activation functions.}
\end{figure}

\subsection{}

\section{Conclusion}

We have reported a reconfigurable photonic-electronic integrated neural network platform that incorporates photonic components on a silicon photonic platform for high-speed information processing. We demonstrated our system is capable of performing essential functionalities such as accurate weighting, summation, and biased nonlinearity for neural network computing. We applied this system for fibre nonlinearity compensation in trans-Pacific transmission data links, which require extremely high throughput and real-time signal processing. We experimentally implemented the second layer of the neural network using silicon photonic neurons, and showed in simulations that the extended system can accurately model the fibre nonlinearity, leading to a Q-factor improvement comparable with numerical simulations on a conventional computer. The system is capable of processing optical signals in the analogue domain, thereby significantly relaxing the stringent requirements of complexity and speed that are typically placed on conventional DSP circuits. Our approach could help in the development of reconfigurable photonic–electronic integrated circuits for solving a range of problems in high-speed intelligent signal processing.

\begin{methods}

\subsection{Device Fabrication}

The device sample is fabricated on a silicon-on-insulator (SOI) wafer with a silicon thickness of 220 nm and a buried oxide thickness of 2 $\mu$m. The waveguide is 500 nm wide. The weight bank consists of four MRRs coupled with two bus waveguides in an add/drop configuration. The four MRRs have radii of 8.0, 8.1, 8.2, 8.3 $\mu$m, respectively. Only the first two MRRs are used in the experiments. A slight difference is introduced in the ring radii to avoid resonance collision. The gap between the ring and bus waveguide is 200 nm, yielding a Q factor of $\sim$6000. For the purpose of weight control, in-ring N-doped photoconductive heaters are used. To implement the N-doped heater, each MRR consists of a ring waveguide etched to a 90 nm thick pedestal that hosts the phosphorus dopants. A 10 $\mu$m wide N doping section is patterned to follow the MRR, outside of which heavy N++ doping is used to make ohmic contacts. Metal vias and traces are deposited to connect the heater contacts of the MRR weight bank to electrical metal pads.

The MRR modulator consists of a ring waveguide with a radius of 8 $\mu$m coupled to two bus waveguides with gaps of 0.2 and 0.5 $\mu$m respectively. The ring waveguide is etched to a 90-nm-thick pedestal that hosts boron and phosphorus to form n+ and n-, and p+ and p- regions, respectively. The lightly p- and n-doped regions are extended to the waveguide core to form a p-n junction for high-speed modulation. The highly p+ or n+ doped regions are located at the pedestal edge to serve as p-type and n-type ohmic contact areas. To fabricate the germanium-on-silicon photodetector, a germanium layer is deposited and patterned on top of the silicon waveguide. Boron and phosphorus are implanted in the sidewalls and slabs of the Ge waveguide to form a horizontal p-i-n junction as well as p-type and n-type ohmic contact areas. Metal vias and traces are deposited to connect the electrical ports of the devices to the electrical metal pads. The cross-sections of MRR resonator embedded with n-doped heater, MRR modulator, and germanium-on-silicon photodetector are illustrated in Supplementary Fig.~1.

The on-chip capacitor is formed by two layers of metal and has a capacitance of 2 pF. The resistor is formed by a slab waveguide etched to a 90 nm thick pedestal that heavily doped with phosphorus (N++). The resistor has a resistance of 300 $\Omega$.

\end{methods}

\begin{addendum}

 \item [Data availability] All data used in this study are available from the corresponding author upon reasonable request.
 \item [Code availability] All codes used in this study are available from the corresponding author upon reasonable request.

\end{addendum}

\begin{addendum}
 \item This research is supported by the Office of Naval Research (ONR) (N00014-18-1-2297), Defense Advanced Research
Projects Agency (HR00111990049), National Science Foundation (NSF) (Grants No. ECCS 1642962 and No. DGE 1148900), and CUHK Research Direct Grant. The devices were fabricated at the IME ASTAR foundry in Singapore. Fabrication support was provided via the Natural Sciences and Engineering Research Council of Canada (NSERC) Silicon Electronic-Photonic Integrated Circuits (SiEPIC) Program and the Canadian Microelectronics Corporation (CMC).

 \item[Author Contributions]  C.H., S.F., and T.L. conceived the ideas and implement the experimental setup, designed the experiment, conducted the experiment measurement, and analyzed the results. C.H., T.L., E.B., S.B., A.J., H.P. designed the silicon photonic chip. A.T., Y.T., F.Y, B.S. provided the theoretical support. T.L., E.B., S.B. performed the chip packaging. C.H., S.F., T.L., A.T., B.S. wrote the manuscript. T.W. and P.P. supervised the research and contributed to the general concept and interpretation of the results. All authors discussed the data and contributed to the manuscript.
 \item[Competing Interests] The authors declare no competing interests.
\item[corresponding authors ] Chaoran Huang~(email: chaoranh@princeton.edu, ORCID: 0000-0001-6997-758X), Paul R. Prucnal~(email: prucnal@princeton.edu, ORCID: 0000-0002-5291-5830).
\end{addendum}

\newpage

\subsection{References}