\section{Introduction}
\indent The Graph Neural Network (GNN) has been one of the heated research subjects in the deep learning community during the last several years.
This is because GNN models can handle graph-structured data via aggregating neighbor information such that it performs among the best graph models in many graph applications such as node classification and link prediction \cite{Kipf2017}. However, its performance has significantly been restricted via its depth due to the over-smoothing problem (i.e. the representations of different nodes in a graph tend to be similar when stacking multiple layers), and the gradient vanishing problem (i.e. the gradient signal decays with the increase of layers) \cite{chen2020graph}, even though some models such as GraphSAGE \cite{hamilton2017inductive} have been proposed to alleviate the dilemma of the over-smoothing problem.\\
\indent On the other hand, the Graph Augmented Multi-Layer Perceptron (GA-MLP) has recently been proposed to handle graph data via augmenting graph features. The GA-MLP extracts n-hop neighbour information ($n=1,2,\cdots$) and aggregate them into an MLP model. Compared with GNN, the GA-MLP model is resistant to the over-smoothing problem \cite{chen2020graph} and has comparable performance. For example, Wu et al. shows that a two-layer GA-MLP approximates the performance of the GNN models on multiple datasets \cite{wu2019simplifying}. Despite the inferior expressiveness to GNNs on some graph tasks \cite{chen2020graph}, GA-MLP can still be considered as an alternative to GNN models.\\
\indent The GA-MLP model is supposed to perform better with the increase of its depth. However, similar to GNN, the GA-MLP still suffers from the gradient-vanishing problem, which is caused by the backpropagation algorithm. Moreover, it prevents the gradients of different layers from being calculated in parallel. This is because the calculation of the gradient in one layer is dependent on its previous layers. \\
\indent To handle these challenges, recently gradient-free optimization methods such as the Alternating Direction Method of Multipliers (ADMM) have been investigated to overcome the difficulty of the backpropagation algorithm. The spirit of ADMM is to decouple a neural network into layerwise subproblems such that each of them can be solved efficiently. ADMM does not require gradient calculation and therefore can avoid the gradient vanishing problem. Existing literature has shown its great potentials.  For example, Talyor et al. and Wang et al. proposed ADMM to train MLP models \cite{taylor2016training,wang2019admm}. Expensive experiments have demonstrated that the ADMM has outperformed most comparison methods.\\
\indent In this paper, we focus on the distributed optimization of the GA-MLP model. Specifically, we apply the recently proposed parallel  ADMM algorithm to solve the GA-MLP model, which splits an MLP model into layerwise training components \cite{wang2020toward}.\\
The contributions in this paper can be summarized as follows:\\
\vspace{-5pt}
\begin{itemize}
\item We propose a novel reformulation of the GA-MLP model, which splits a neural network into independent layer partitions and allows for ADMM to achieve model parallelism.
\item We propose a pdADMM framework to train a GA-MLP model. All subproblems generated by the ADMM algorithm are discussed. The extended pdADMM-Q algorithm reduces communication cost using the quantization technique.
\item  We provide theoretical guarantee of the proposed pdADMM algorithm and the pdADMM-Q algorithm. Specifically, they converge to a critical point when the hyperparameters are sufficiently large, and their sublinear convergence rate is $o(1/k)$.
\item We conduct extensive experiments in six benchmark datasets to show the high speedup of the pdADMM, as well as its outstanding performance when compared with all state-of-the-art optimizers.
\end{itemize}
\indent The organization of this paper is shown as follows: In Section \ref{sec:related work}, we summarize recent related research work to this paper. In Section \ref{sec:methods}, we propose the pdADMM algorithm and the pdADMM-Q algorithm to train the GA-MLP model. Section \ref{sec:convergence} details the convergence results of the proposed pdADMM algorithm and the pdADMM-Q algorithm.
Extensive experiments on benchmark datasets to demonstrate the speedup, and outstanding performance of the pdADMM are shown in Section \ref{sec:experiments}, and Section \ref{sec:conclusion} concludes this work.
\section{Related Work}
\label{sec:related work}
\indent This section summarizes existing literatures related to this research.\\
\indent \textbf{Distributed Deep Learning} With the increase of public datasets and layers of neural networks, it is imperative to establish distributed deep learning systems for large-scale applications. Many systems have been established to satisfy such needs. Famous systems include  Terngrad  \cite{wen2017terngrad}, Horovod \cite{sergeev2018horovod}, SINGA \cite{ooi2015singa} Mxnet \cite{chen2015mxnet}, TicTac \cite{hashemitictac} and Poseidon \cite{zhang2017poseidon}. They apply some parallelism techniques to reduce communication overhead and computational time. The parallelism technique can be categorized into two classes: data parallelism and model parallelism.
Data parallelism focuses on distributing data across different processors, and then aggregate results from them into a server. Scaling SGD is one of the most common ways to reach data parallelism \cite{zinkevich2010parallelized}. For example, the distributed architecture, Poseidon, is achieved by scaling SGD through overlapping communication and computation over networks.  The recently proposed ADMM \cite{taylor2016training,wang2019admm} is another way to achieve data parallelism. However, data parallelism suffers from the bottleneck of a neural network: for SGD, the gradient should be transmitted through all processors;  for ADMM, the parameters in one layer are subject to those in its previous layer. As a result, this leads to heavy communication cost and time delay. Model parallelism, however, can solve this challenge because model parallelism splits a neural network into many independent partitions. In this way, each partition can be optimized independently and reduce layer dependency. For instance, Parpas and Muir proposed a parallel-in-time method from the perspective of dynamic systems \cite{parpas2019predict}; Huo et al. introduced a feature replay algorithm to achieve model parallelism \cite{huo2018training}. Zhuang et al. broke layer dependency by introducing the delayed gradient \cite{zhuang2019fully}.\\
\indent \textbf{Deep Learning on Graphs:} Graphs are ubiquitous structures and are popular in real-world applications. There is a surge of interest to apply deep learning techniques on graphs. For a comprehensive summary please refer to \cite{wu2020comprehensive}. Graph Convolution Network (GCN) is the most popular model which performs well on downstream tasks only with two layers \cite{kipf2016semi}. However, it suffers from over-smoothing problems and many approaches have been proposed to reduce its effect \cite{hamilton2017inductive,zhao2019pairnorm,chen2020graph}. GA-MLP is such an alternative we focus on in this paper.
\section{Methods}
\label{sec:methods}
\indent We propose the pdADMM framework to GA-MLP models in this section. Specifically, Section \ref{sec:problem formulation} formulates the GA-MLP model training problem and the pdADMM framework, and Section \ref{sec:subproblem} details solutions to all subproblems of the  pdADMM framework.
\subsection{Problem Formulation}
\label{sec:problem formulation}
\begin{table}
\scriptsize
\centering
 \begin{tabular}{cc}
 \hline
 Notations&Descriptions\\ \hline
 $L$& Number of layers.\\
 $W_l$& The weight matrix for the $l$-th layer.\\
 $b_l$& The intercept vector for the $l$-th layer.\\
 $z_l$& The auxiliary variable of the linear mapping for the $l$-th layer.\\
 $f_l(z_l)$& The nonlinear activation function for the $l$-th layer.\\
 $p_l$& The input for the $l$-th layer.\\
 $q_l$& The output for the $l$-th layer.\\
 $X$& The node representation of the graph.\\
 $A$& The adjacency matrix of the graph.\\
 $y$& The predefined label vector.\\
 $R(z_L,y)$& The risk function for the $l$-th layer.\\
 $n_l$& The number of neurons for the $l$-th layer.\\
\hline
  \end{tabular}
  \captionof{table}{Important Notations}
   \label{tab:notation}

\end{table}
Let $X$ be a matrix of feature representation in the graph $G=(V, E)$, where $V$ and $E$ are sets of vertices and edges, respectively, and $\Psi=\{\psi_1(A),\cdots,\psi_K(A)\}$ be a set of (usually multi-hop) linear operators that are functions of the adjacency matrix $A$.
${X}_k=\psi_k(A)\phi(X)$ is the augmentation of node features by the $k$-hop linear operator, where $\phi(\bullet)$ is a learnable function acting as a feature transformation applied to each node separately. Then the model concatenates $X=[{X}_1,\cdots,{X}_K]$, where $X$ is denoted as the input of the MLP model, which is formulated as follows \cite{wang2019admm}:
\begin{problem}
\label{prob:problem 1}
\begin{align*}
     & \min\nolimits_{W_l,b_l,z_l,p_l} R(z_L;y)+\sum\nolimits_{l=1}^L\Omega_l(W_l) \\
     &s.t.\ z_l=W_lp_{l}+b_l, \   p_{l+1}=f_l(z_l)(l=1,\cdots,L-1)
\end{align*}
\end{problem}
where $p_1={X}$ is the input of the deep neural network, and $y$ is a predefined label vector. $p_l$ is the input for the $l$-th layer, also the output for the $(l-1)$-th layer. $R(z_L;y)$ is a risk function for the $L$-th layer, which is convex, continuous and proper; $\Omega_l(W_l)$ is a convex regularization term such as the $\ell_2$ penalty. $z_l=W_lp_l+b_l$  and $p_{l+1}=f_l(z_l)$ are linear and nonlinear mappings for the $l$-th layer, respectively.\\
\indent Problem \ref{prob:problem 1} has been addressed by deep learning Alternating Direction Method of Multipliers (dlADMM) \cite{wang2019admm}. However, parameters in one layer are dependent on its neighboring layers, and hence can not achieve parallelism. For example, the update of $p_{l+1}$ on the $l+1$-th layer needs to wait before $z_l$ on the $l$-th layer is updated.
\label{sec:pdADMM}
In order to address layer dependency, we relax Problem \ref{prob:problem 1} to Problem \ref{prob:problem 2} as follows:
\begin{problem}
\label{prob:problem 2}
\begin{align*}
    &\min\nolimits_{\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q}} F(\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q})=R(z_L;y)+\sum\nolimits_{l=1}^L\Omega(W_l)\\&+(\nu/2)(\sum\nolimits_{l=1}^{L}\Vert z_l-W_lp_l-b_l\Vert^2_2+\sum\nolimits_{l=1}^{L-1}\Vert q_l-f_l(z_l)\Vert^2_2)\\&s.t. \ p_{l+1}=q_l
\end{align*}
\end{problem}
where   $\textbf{p}=\{p_l\}_{l=1}^{L}$, $\textbf{W}=\{W_l\}_{l=1}^{L}$, $\textbf{b}=\{b_l\}_{l=1}^{L}$, $\textbf{z}=\{z_l\}_{l=1}^{L}$, $\textbf{q}=\{q_l\}_{l=1}^{L-1}$, and $\nu>0$ is a tuning parameter. As $\nu\rightarrow \infty$, Problem \ref{prob:problem 2} approaches Problem  \ref{prob:problem 1}. We reduce layer dependency by splitting the output of the $l$-th layer and the input of the ($l+1$)-th layer into two variables $p_{l+1}$ and $q_{l}$, respectively.

\indent  The high-level overview of the pdADMM algorithm is shown in Figure \ref{fig:pdADMM framework}.  Specifically,  the whole GA-MLP is split into multiple layers, each of which can be optimized by an independent client. Therefore, the layerwise training can be implemented in parallel. Moreover, the gradient vanishing problem can be avoided in this way. This is because the accumulated gradient calculated by the backpropagation algorithm is split into layerwise components. \\
\indent Now we follow the ADMM routine to solve Problem \ref{prob:problem 2}. The augmented Lagrangian function is formulated mathematically as follows:
\begin{align*}
    &L_\rho(\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q},\textbf{u})\\&\!=\!F(\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q})\!+\!\sum\nolimits_{l\!=\!1}^{L\!-\!1}(u_l^T(p_{l\!+\!1}\!-\!q_l)\!+\!(\rho/2)\Vert p_{l+1}-q_l\Vert^2_2)\\&=\!R(z_L;y)\!+\!\phi(p_1,W_1,b_1,z_1)\!+\!\sum\nolimits_{l\!=\!2}^{L}\!\phi(p_l,W_l,b_l,z_l,q_{l-1},u_{l\!-\!1})\!\\&+\!(\nu/2)\sum\nolimits_{l\!=\!1}^{L-1}\Vert q_l\!-\!f_l(z_l)\Vert^2_2
\end{align*}
where $\phi(p_1,W_1,b_1,z_1)=(\nu/2)\Vert z_1-W_1p_1-b_1\Vert^2_2$, $\phi(p_l,W_l,b_l,z_l,q_{l-1},u_{l-1})=(\nu/2)\Vert z_l-W_lp_l-b_l\Vert^2_2+u^T_{l-1}(p_l-q_{l-1})+(\rho/2)\Vert p_l-q_{l-1}\Vert^2_2$,  $u_l(l=1,\cdots,L-1)$ are dual variables, $\rho>0$ is a parameter, and $\textbf{u}=\{u_l\}_{l=1}^{L-1}$.
The detail of the pdADMM is shown in Algorithm \ref{algo:distributed ADMM}. Specifically, Lines 5-9 update primal variables $\textbf{p}$, $\textbf{W}$, $\textbf{b}$, $\textbf{z}$ and $\textbf{q}$, respectively, while Line 11 updates the dual variable $\textbf{u}$. The discussion on how to solve subproblems generated by the pdADMM is detailed in the next section.
\begin{figure}
   \centering
    \includegraphics[width=\linewidth]{framework.png}
    \caption{The overall pdADMM optimization framework: it splits a neural network into layerwise components.}
    \label{fig:pdADMM framework}
\end{figure}
\begin{algorithm}
\small
\caption{The pdADMM Algorithm to Solve Problem \ref{prob:problem 2}}
\begin{algorithmic}
\label{algo:distributed ADMM}
\REQUIRE $y$, $p_1=X$, $\rho$, $\nu$.
\ENSURE $\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q}$.
\STATE Initialize $k=0$.
\WHILE{$\textbf{p}^k,\textbf{W}^{k},\textbf{b}^{k},\textbf{z}^{k},\textbf{q}^{k}$ not converged}
\STATE Update $p_l^{k+1}$ of different $l$ by Equation \eqref{eq:update p} in parallel.
\STATE Update $W_l^{k+1}$ of different $l$ by Equation \eqref{eq:update W} in parallel.
\STATE Update $b_l^{k+1}$ of different $l$ by Equation \eqref{eq:update b} in parallel.
\STATE Update $z_l^{k+1}$ of different $l$ by Equations \eqref{eq:update z} and \eqref{eq:update zl} in parallel.
\STATE Update $q_l^{k+1}$ of different $l$ by Equation \eqref{eq:update q} in parallel.
\STATE $r^k_l\leftarrow p^{k+1}_{l+1}-q^{k+1}_l(l=1,\cdots,L)$ in parallel \# Compute residuals.
\STATE Update $u_l^{k+1}$ of different $l$ by Equation \eqref{eq:update u} in parallel.
\STATE $k\leftarrow k+1$.
\ENDWHILE
\STATE Output $\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q}$.
\end{algorithmic}
\end{algorithm}
\subsection{Solutions to All Subproblems}
\label{sec:subproblem}
\indent In this section, we discuss how to solve all subproblems generated by pdADMM in detail.\\
\textbf{1. Update $\textbf{p}^{k+1}$}\\
\indent The variable $\textbf{p}^{k+1}$ is updated as follows:
\begin{align*}
    p^{k+1}_l&\leftarrow \arg\min\nolimits_{p_l} L_\rho(\textbf{p},\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{q}^k,\textbf{u}^k)\\&=\phi(p_l,W^k_l,b^k_l,z^k_l,q^k_{l-1},u^k_{l-1})
\end{align*}
Because $W_l$ and $p_l$ are coupled in $\phi$, solving $p_l$ should require the time-consuming operation of matrix inversion of $W_l$. To handle this, we apply similar quadratic approximation techniques as used in dlADMM \cite{wang2019admm} as follows:
\begin{align}
p_l^{k+1}&\leftarrow\arg\min_{p_l} U_l(p_l;\tau_l^{k+1})
    \label{eq:update p}
\end{align}
where $
    U_l(p_l;\tau^{k+1}_l)=\phi(p^k_l,W^k_l,b^k_l,z^k_l,q^k_{l-1},u^k_{l-1})+\nabla_{p^k_l}\phi^T(p^k_l,W^k_l,b^k_l,z^k_l,q^k_{l-1},u^k_{l-1})(p_l-p^k_l)+({\tau}^{k+1}_l/2)\Vert p_{l}-p^k_{l}\Vert^2_2
$, and $\tau^{k+1}_l>0$ is a parameter. $\tau_l^{k+1}$ should satisfy $\phi(p^{k+1}_l,W^k_l,b^k_l,z^k_l,q^k_{l-1},u^k_{l-1})\leq U_l(p^{k+1}_l;\tau^{k+1}_l)$.  The solution to Equation \eqref{eq:update p} is:
${p}^{k+1}_{l}\leftarrow p^k_{l}-\nabla_{p^k_{l}}\phi(p^k_l,W^k_l,b^k_l,z^k_l,q^k_{l-1},u^k_{l-1})/{\tau}^{k+1}_{l}$.\\
\textbf{2. Update $\textbf{W}^{k+1}$}\\
\indent The variable $\textbf{W}^{k+1}$ is updated as follows:
\begin{align*}
    W^{k+1}_l&\leftarrow \arg\min\nolimits_{W_l} L_\rho(\textbf{p}^{k+1},\textbf{W},\textbf{b}^k,\textbf{z}^{k},\textbf{q}^{k},\textbf{u}^{k})\\& =\arg\min\nolimits_{W_l} \begin{cases} \phi(p^{k+1}_1,W_1,b^k_1,z^k_1)
    &l=1\\\phi(p^{k\!+\!1}_l,W_l,b^k_l,z^k_l,q^k_{l\!-\!1},u^k_{l\!-\!1})&1<\!l\!\leq\! L
    \end{cases}
\end{align*}
Similar to updating $p_l$, the following subproblem should be solved instead:
\begin{align}
    W_l^{k+1}&\leftarrow\arg\min\nolimits_{W_l} V_l(W_l;\theta_l^{k+1}) \label{eq:update W}
\end{align}
where 
\begin{align*}
&V_1(W_1;\theta^{k+1}_1)= \phi(p^{k+1}_1,W^k_1,b^k_1,z^k_1)\\&+\!\nabla_{W^k_1}\phi^T(p^{k+1}_1,W^k_1,b^k_1,z^k_1)(W_1\!-\!W^k_1)   \\&+({\theta}^{k+1}_l/2)\Vert W_{1}-W^k_{1}\Vert^2_2\\&
 V_l(W_l;\theta^{k+1}_l)=
   \phi(p^{k+1}_l,W^k_l,b^k_l,z^k_l,q^k_{l-1},u^k_{l-1})\\&+\nabla_{W^k_l}\phi^T(p^{k+1}_l,W^k_l,b^k_l,z^k_l,q^k_{l-1},u^k_{l-1})(W_l-W^k_l)\\&+({\theta}^{k+1}_l/2)\Vert W_{l}-W^k_{l}\Vert^2_2
\end{align*}
and $\theta^{k+1}_l$ is a parameter, which should satisfy $\phi(p^{k+1}_1,W^{k+1}_1,b^k_1,z^k_1)\leq V(W^{k+1}_1;\theta^{k+1}_1)$ and $\phi(p^{k+1}_l,W^{k+1}_l,b^k_l,z^k_l,q^k_{l-1},u^k_{l-1})\leq V(W^{k+1}_l;\theta^{k+1}_l)(1<l<L)$. The solution to Equation \eqref{eq:update W} is shown as follows:
\begin{align*}
    {W}^{k\!+\!1}_{l}\!\leftarrow\! W^k_{l}\!-\!\begin{cases}\nabla_{W^k_1}\phi(p^{k+1}_1\!,\!W^k_1\!,\!b^k_1\!,\!z^k_1)/\theta^{k+1}_l& l=1\\\nabla_{W^k_l}\phi(p^{k\!+\!1}_l\!,\!W^k_l\!,\!b^k_l\!,\!z^k_l\!,\!q^k_{l\!-\!1}\!,\!u^k_{l\!-\!1})/\theta^{k\!+\!1}_l& 1\!<l\!\leq\! L
    \end{cases}
\end{align*}
\textbf{3. Update $\textbf{b}^{k+1}$}\\
\indent The variable $\textbf{b}^{k+1}$ is updated as follows:
\begin{align*}
    b^{k+1}_l&\leftarrow \arg\min\nolimits_{b_l} L_\rho(\textbf{p}^{k+1},\textbf{W}^{k+1},\textbf{b},\textbf{z}^{k},\textbf{q}^{k},\textbf{u}^{k})\\& =\arg\min\nolimits_{b_l} \begin{cases} \phi(p^{k+1}_1,W^{k+1}_1,b_1,z^k_1)
    &l\!=\!1\\\phi(p^{k+1}_l,W^{k\!+\!1}_l,b_l,z^k_l,q^k_{l\!-\!1},u^k_{l\!-\!1})&1<\!l\!\leq L
    \end{cases}.
\end{align*}
Similarly, we solve the following subproblems instead:
\begin{align*}
&b^{k+1}_1\leftarrow \arg\min\nolimits_{b_1} \phi(p^{k+1}_1,W^{k+1}_1,b^k_1,z^k_1)\\&+\nabla_{b^k_1}\phi^T(p^{k+1}_1,W^{k+1}_1,b^k_1,z^k_1)(b_l-b^k_l)+(\nu/2)\Vert b_{l}-b^k_{l}\Vert^2_2
\end{align*}
\begin{align}
&\nonumber b^{k+1}_l\leftarrow \arg\min\nolimits_{b_l} \phi(p^{k+1}_l,W^{k+1}_l,b^k_l,z^k_l,q^k_{l-1},u^k_{l-1})\\&\nonumber+\nabla_{b^k_l}\phi^T(p^{k+1}_l,W^{k+1}_l,b^k_l,z^k_l,q^k_{l-1},u^k_{l-1})(b_l-b^k_l)\\&+(\nu/2)\Vert b_{l}-b^k_{l}\Vert^2_2 (1< l\leq L) \label{eq:update b}
\end{align}
\indent The solution to Equation \eqref{eq:update b} is:
\begin{align*}
    {b}^{k\!+\!1}_{l}\leftarrow b^k_{l}\!-\!\begin{cases}\nabla_{b^k_1}\phi(p^{k\!+\!1}_1,W^{k\!+\!1}_1,b^k_1,z^k_1)/\nu& l=1\\\nabla_{b^k_l}\phi(p^{k\!+\!1}_l,W^{k\!+\!1}_l,b^k_l,z^k_l,q^k_{l\!-\!1},u^k_{l\!-\!1})/\nu& 1\!<\!l\!\leq L
    \end{cases}
\end{align*}
\textbf{4. Update $\textbf{z}^{k+1}$}\\
\indent The variable
$\textbf{z}^{k+1}$ is updated as follows:
\begin{align}
    \nonumber& z^{k\!+\!1}_l\!\leftarrow\! \arg\min\nolimits_{z_l} (\nu/2)\Vert z_l\!-\!W^{k\!+\!1}_l\!p^{k\!+\!1}_l\!-\!b^{k\!+\!1}_l\Vert^2_2\\&+\!(\nu/2)\Vert q^k_l\!-
    \!f_l(z_l)\Vert^2_2\!+\!(\nu/2)\Vert z_l\!-\!z^k_l\Vert^2_2(l\!<\! L)\label{eq:update z}\\
    &z^{k\!+\!1}_L\!\leftarrow\! \arg\min\nolimits_{z_l}\! R(z_L;y)\!+\!(\nu/2)\Vert z_L\!-\!W^{k\!+\!1}_Lp^{k\!+\!1}_L\!-\!b^{k\!+\!1}_L\Vert^2_2\! \label{eq:update zl}
\end{align}
where a quadratic term $(\nu/2)\Vert z_l-z^k_l\Vert^2_2$ is added in Equation \eqref{eq:update z} to control $z^{k+1}_l$ to close to $z^{k}_l$.
Equation \eqref{eq:update zl} is convex, which can be solved by Fast Iterative Soft Thresholding Algorithm (FISTA) \cite{beck2009fast}.\\
For Equation \eqref{eq:update z}, nonsmooth activations usually lead to  closed-form solutions \cite{wang2019admm,wang2017nonconvex}. For example, for Relu $f_l(z_l)=\max(z_l,0)$, the solution to Equation \eqref{eq:update z} is shown as follows:
\begin{align*}
    z^{k+1}_l=
    \begin{cases} \min((W^{k\!+\!1}_{l}p^{k\!+\!1}_{l}\!+\!b^{k+1}_l\!+\!z^k_l)/2,0)&z^{k\!+\!1}_l\leq0\\
    \max((W^{k\!+\!1}_{l}p^{k\!+\!1}_{l}\!+\!b^{k\!+\!1}_l\!+\!q^k_l\!+\!z^k_l)/3,0)&z^{k+1}_l\geq0\\
    \end{cases}
\end{align*}
For smooth activations such as tanh and sigmoid, a lookup-table is recommended \cite{wang2019admm}.\\
\textbf{5. Update $\textbf{q}^{k+1}$}\\
\indent The variable
$\textbf{q}^{k+1}$ is updated as follows:
\begin{align}
     \nonumber q^{k+1}_l&\leftarrow \arg\min\nolimits_{q_l} L_\rho(\textbf{p}^{k+1},\textbf{W}^{k+1},\textbf{b}^{k+1},\textbf{z}^{k+1},\textbf{q},\textbf{u}^k)\\&=\arg\min\nolimits_{q_l}\phi(p^{k+1}_{l+1},W^{k+1}_{l+1},b^{k+1}_{l+1},z^{k+1}_{l+1},q_l,u^{k}_l). \label{eq:update q}
\end{align}
Equation \eqref{eq:update q} has a closed-form solution as follows:\\
\begin{align*}
    q^{k+1}_l\leftarrow (p_{l+1}^{k+1}+u^k_l/\rho+f_l(z^{k+1}_l))/2
\end{align*}
\textbf{6. Update $\textbf{u}^{k+1}$}\\
\indent The variable
$\textbf{u}^{k+1}$ is updated as follows:
\begin{align}
u^{k+1}_l&\leftarrow u^k_l+\rho(p^{k+1}_{l+1}-q^{k+1}_{l}) \label{eq:update u}
\end{align}
\indent Finally, Our proposed pdADMM can be efficient for training a deep GA-MLP model. To achieve this,  we begin from training a shallow GA-MLP model with the first few layers of the deep GA-MLP, then more layers are added for training step by step until finally all layers are involved in the training process. The pdADMM can achieve  excellent performance as well as reduce training cost by this strategy.
\subsection{pdADMM Extension for Quantization (pdADMM-Q)}
\indent In the proposed pdADMM algorithm, $p_l$ and $q_l$ are transmitted back and forth among layers (i.e. clients). However, the communication overheads of $p_l$ and $q_l$ surge for a large-scale  graph $G$ with millions of nodes. To alleviate this challenge, the quantization technique is commonly utilized to reduce communication cost by mapping continuous values into a discrete set \cite{huang2021alternating}. In other words, $p_l$ and $q_l$ are required to fit into a countable set $\Delta$, which is shown as follows:
\begin{problem}
\label{prob:problem 3}
\begin{align*}
    & \min\nolimits_{\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q}} F(\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q})=R(z_L;y)+\sum\nolimits_{l=1}^L\Omega(W_l)\\&+(\nu/2)(\sum\nolimits_{l=1}^{L}\Vert z_l-W_lp_l-b_l\Vert^2_2+\sum\nolimits_{l=1}^{L-1}\Vert q_l-f_l(z_l)\Vert^2_2)\\&s.t. \ p_{l+1}=q_l, \ p_{l}\in \Delta=\{\delta_1,\cdots,\delta_m\}
\end{align*}
\end{problem}
where $\delta_i(i=1,\cdots,m)\in \Delta$ are quantized values, which can be integers or  low-precision values. $m=\vert\Delta\vert$ is the cardinality of $\Delta$. To address Problem \ref{prob:problem 3}, we rewrite it into the following form:
\begin{align*}
    &\min\nolimits_{\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q}} R(z_L;y)+\sum\nolimits_{l=1}^L\Omega(W_l)+\sum\nolimits_{l=2}^L\mathbb{I}(p_l)\\&+(\nu/2)(\sum\nolimits_{l=1}^{L}\Vert z_l-W_lp_l-b_l\Vert^2_2+\sum\nolimits_{l=1}^{L-1}\Vert q_l-f_l(z_l)\Vert^2_2)\\&s.t. \ p_{l+1}=q_l
\end{align*}
where the indicator function $\mathbb{I}(p_l)$ is defined as follows: $\mathbb{I}(p_l)=0$ if $p_l\in\Delta$, and $\mathbb{I}(p_l)=+\infty$ if $p_l\not\in\Delta$. The Augmented Lagrangian of Problem \ref{prob:problem 3} is shown as follows:
\begin{align*}
    \beta_\rho(\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q},\textbf{u})=L_\rho(\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q},\textbf{u})+\sum\nolimits_{l=2}^L \mathbb{I}(p_l)
\end{align*}
where $L_\rho$ is the Augmented Lagrangian of Problem \ref{prob:problem 2}. The extended pdADMM-Q algorithm follows the same routine as the  pdADMM algorithm, which is shown as follows:\\
\begin{algorithm}
\small
\caption{The pdADMM-Q Algorithm to Solve Problem \ref{prob:problem 3}}
\begin{algorithmic}
\label{algo:distributed ADMM-Q}
\REQUIRE $y$, $p_1=X$, $\rho$, $\nu$.
\ENSURE $\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q}$.
\STATE Initialize $k=0$.
\WHILE{$\textbf{p}^k,\textbf{W}^{k},\textbf{b}^{k},\textbf{z}^{k},\textbf{q}^{k}$ not converged}
\STATE Update $p_l^{k+1}$ of different $l$ by Equation \eqref{eq:update p quantized} in parallel.
\STATE Update $W_l^{k+1}$ of different $l$ by Equation \eqref{eq:update W} in parallel.
\STATE Update $b_l^{k+1}$ of different $l$ by Equation \eqref{eq:update b} in parallel.
\STATE Update $z_l^{k+1}$ of different $l$ by Equations \eqref{eq:update z} and \eqref{eq:update zl} in parallel.
\STATE Update $q_l^{k+1}$ of different $l$ by Equation \eqref{eq:update q} in parallel.
\STATE $r^k_l\leftarrow p^{k+1}_{l+1}-q^{k+1}_l(l=1,\cdots,L)$ in parallel \# Compute residuals.
\STATE Update $u_l^{k+1}$ of different $l$ by Equation \eqref{eq:update u} in parallel.
\STATE $k\leftarrow k+1$.
\ENDWHILE
\STATE Output $\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q}$.
\end{algorithmic}
\end{algorithm}
\indent Obivously, the only difference between the pdADMM-Q algorithm and the pdADMM algorithm is the $p_l$-subproblem in Line 6, which is outlined in the following:\\
\begin{align}
p_l^{k+1}&\leftarrow\arg\min\nolimits_{p_l} U_l(p_l;\tau_l^{k+1})+\mathbb{I}(p_l)
    \label{eq:update p quantized}
\end{align}
where $U_l$ follows Equation \eqref{eq:update p}. The solution to Equation \eqref{eq:update p quantized} is \cite{huang2021alternating}:
${p}^{k+1}_{l}\leftarrow\arg\min\nolimits_{\delta\in \Delta}\Vert \delta-(p^k_{l}-\nabla_{p^k_{l}}\phi(p^k_l,W^k_l,b^k_l,z^k_l,q^k_{l-1},u^k_{l-1})/{\tau}^{k+1}_{l})\Vert$
\section{Convergence Analysis}
\label{sec:convergence}
\indent In this section, the theoretical convergence of the proposed pdADMM algorithm and the pdADMM-Q algorithm is provided, whose proofs are detailed in the Appendix. Firstly, the Lipschitz continuity and coercivity are defined as follows:
\begin{definition}(Lipschitz Continuity)
A function $g(x)$ is Lipschitz continuous if there exists a constant $D>0$ such that $\forall x_1,x_2$, the following holds
\begin{align*}
    \Vert g(x_1)-g(x_2)\Vert\leq D\Vert x_1-x_2\Vert. 
\end{align*}
\end{definition}
\begin{definition}(Coercivity)
A function $h(x)$ is coerce over the feasible set  $\mathscr{F}$ means that $h(x)\rightarrow \infty$ if $x\in \mathscr{F}$ and $\Vert x\Vert\rightarrow \infty$.
\end{definition}
Next, the definition of a quantized critical point is shown as follows \cite{huang2021alternating}:
\begin{definition}(quantized critical point)
The $p_l$ is a quantized critical point of of Problem \ref{prob:problem 3} if there exists $\tau>0$ such that
\begin{align*}
    p_l\in\arg\min\nolimits_{\delta\in \Delta}\Vert \delta-(p_l-\nabla_{p_l} F(\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q})/\tau)\Vert
\end{align*}
\end{definition}
The quantized critical point is an extension of the critical point in the discrete setting, and any global solution $p_l$ to Problem \ref{prob:problem 3} is a quantized critical point to Problem \ref{prob:problem 3} (Lemma 3.7 in \cite{huang2021alternating}). Then the following assumption is required for convergence analysis.
\begin{assumption}
$f_l(z_l)$ is Lipschitz continuous with coefficient $S>0$, and $F(\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q})$ is coercive. Moreover, $\partial f_l(z_l)$ is bounded, i.e. there exists $M>0$ such that $\Vert\partial f_l(z_l)\Vert \leq M$. 
\label{ass:lipschitz continuous}
\end{assumption}
Assumption \ref{ass:lipschitz continuous} is mild to satisfy: most common activation functions such as Relu and leaky Relu satisfy Assumption \ref{ass:lipschitz continuous}.
No assumption is needed on the risk function $R(z_l;y)$, which shows that the convergence condition of our proposed pdADMM is milder than that of the dlADMM, which requires $R(z_l;y)$ to be Lipschitz differentiable \cite{wang2019admm}. Due to space limit, the detailed proofs are provided in the Appendix. The technical proofs follow the similar routine as dlADMM \cite{wang2019admm}. The difference consists in the fact that the dual variable $u_l$ is controlled by $q_l$ and $z_l$ (Lemma \ref{lemma:u bound} in the Appendix), which holds under Assumption \ref{ass:lipschitz continuous}, while $u_l$ can  be controlled only by $z_l$ in the convergence proof of dlADMM. The first lemma shows that the objective keeps decreasing when $\rho$ is sufficently large.

\begin{lemma}[Decreasing Objective]
\label{lemma:objective decrease} For both the pdADMM algorithm and the pdADMM-Q algorithm, 
if $\rho>\max(4\nu S^2,(\sqrt{17}+1)\nu/2)$, there exist $C_1=\nu/2-2\nu^2S^2/\rho>0$ and $C_2=\rho/2-2\nu^2/\rho-\nu/2>0$ such that it holds for any $k\in \mathbb{N}$ that
\begin{align}
    \nonumber &L_\rho(\textbf{p}^{k},\!\textbf{W}^{k},\!\textbf{b}^k,\!\textbf{z}^{k},\!\textbf{q}^k,\!\textbf{u}^{k})\!-\!L_\rho(\textbf{p}^{k\!+\!1},\textbf{W}^{k\!+\!1},\textbf{b}^{k\!+\!1},\textbf{z}^{k\!+\!1},\textbf{q}^{k\!+\!1},\textbf{u}^{k\!+\!1})\\\nonumber &\geq \sum\nolimits_{l\!=\!2}^L (\tau^{k\!+\!1}_l/2)\Vert p^{k+1}_l\!-\!p^k_l\Vert^2_2\!+\!\sum\nolimits_{l\!=\!1}^{L}(\theta^{k\!+\!1}_l/2)\Vert  W^{k\!+\!1}_l\!-\!W^k_l\Vert^2_2\!\\&\nonumber\!+\!\sum\nolimits_{l\!=\!1}^{L}(\nu/2)\Vert  b^{k\!+\!1}_l\!-\!b^k_l\Vert^2_2\!+\!\sum\nolimits_{l\!=\!1}^{L\!-\!1} C_1\Vert z^{k+1}_l\!-\!z^k_l\Vert^2_2\\&\!+\!(\nu/2)\Vert z^{k+1}_L-z^k_L\Vert^2_2+\sum\nolimits_{l=1}^{L-1}C_2\Vert q^{k+1}_l-q^k_l\Vert^2_2
    \label{eq:pdadmm objective decrease}
\\\nonumber &\beta_\rho(\textbf{p}^{k},\!\textbf{W}^{k},\!\textbf{b}^k,\!\textbf{z}^{k},\!\textbf{q}^k,\!\textbf{u}^{k})\!-\!\beta_\rho(\textbf{p}^{k\!+\!1},\textbf{W}^{k\!+\!1},\textbf{b}^{k\!+\!1},\textbf{z}^{k\!+\!1},\textbf{q}^{k\!+\!1},\textbf{u}^{k\!+\!1})\\\nonumber &\geq \!\sum\nolimits_{l\!=\!1}^{L}(\theta^{k\!+\!1}_l/2)\Vert  W^{k\!+\!1}_l\!-\!W^k_l\Vert^2_2\!+\!\sum\nolimits_{l\!=\!1}^{L}(\nu/2)\Vert  b^{k\!+\!1}_l\!-\!b^k_l\Vert^2_2\!\\\nonumber&+\!\sum\nolimits_{l\!=\!1}^{L\!-\!1} C_1\Vert z^{k+1}_l\!-\!z^k_l\Vert^2_2+\!(\nu/2)\Vert z^{k+1}_L-z^k_L\Vert^2_2\\&+\sum\nolimits_{l=1}^{L-1}C_2\Vert q^{k+1}_l-q^k_l\Vert^2_2
    \label{eq:pdadmm-q objective decrease}
\end{align}
\end{lemma}
\indent The second Lemma illustrates that the objective is bounded from below when $\rho$ is large enough, and all variables are bounded.
\begin{lemma} [Bounded Objective]
\label{lemma:lower bounded}

(1). For the pdADMM algorithm, if $\rho> \nu$, then
 $L_\rho(\textbf{p}^k,\textbf{W}^{k},\textbf{b}^{k},\textbf{z}^{k},\textbf{q}^{k},\textbf{u}^{k})$
 is lower bounded. Moreover, $\textbf{p}^k,\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{q}^{k}$,and $\textbf{u}^k$ are bounded, i.e. there exist $\mathbb{N}_\textbf{p}$, $\mathbb{N}_\textbf{W}$, $\mathbb{N}_\textbf{b}$,  $\mathbb{N}_\textbf{z}$,  $\mathbb{N}_\textbf{q}$, and $\mathbb{N}_\textbf{u}>0$, such that $\Vert \textbf{p}^k\Vert\leq \mathbb{N}_\textbf{p}$, $\Vert \textbf{W}^k\Vert\leq \mathbb{N}_\textbf{W}$, $\Vert \textbf{b}^k\Vert\leq \mathbb{N}_\textbf{b}$, $\Vert \textbf{z}^k\Vert\leq \mathbb{N}_\textbf{z}$, $\Vert \textbf{q}^k\Vert\leq \mathbb{N}_\textbf{q}$, and $\Vert \textbf{u}^k\Vert\leq \mathbb{N}_\textbf{u}$. \\
 (2). For the pdADMM-Q algorithm, If $\rho> \nu$, then
 $\beta_\rho(\textbf{p}^k,\textbf{W}^{k},\textbf{b}^{k},\textbf{z}^{k},\textbf{q}^{k},\textbf{u}^{k})$
 is lower bounded. Moreover, $\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{q}^{k}$,and $\textbf{u}^k$ are bounded, i.e. there exist  $\mathbb{N}_\textbf{W}$, $\mathbb{N}_\textbf{b}$,  $\mathbb{N}_\textbf{z}$,  $\mathbb{N}_\textbf{q}$, and $\mathbb{N}_\textbf{u}>0$, such that $\Vert \textbf{W}^k\Vert\leq \mathbb{N}_\textbf{W}$, $\Vert \textbf{b}^k\Vert\leq \mathbb{N}_\textbf{b}$, $\Vert \textbf{z}^k\Vert\leq \mathbb{N}_\textbf{z}$, $\Vert \textbf{q}^k\Vert\leq \mathbb{N}_\textbf{q}$, and $\Vert \textbf{u}^k\Vert\leq \mathbb{N}_\textbf{u}$.\end{lemma}
 \indent Based on Lemmas \ref{lemma:objective decrease} and \ref{lemma:lower bounded}, the following theorem ensures that the objective  is convergent.
 \begin{theorem}[Convergent Objective] \label{theo:convergent variable}
(1). For the pdADMM algorithm, if $\rho>\max(4\nu S^2,(\sqrt{17}+1)\nu/2)$, then
$L_\rho(\textbf{p}^k,\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{q}^{k},\textbf{u}^k)$ is convergent. Moreover, $\lim_{k\rightarrow\infty}\Vert\textbf{p}^{k+1}-\textbf{p}^{k}\Vert^2_2=0$, $\lim_{k\rightarrow\infty}\Vert\textbf{W}^{k+1}-\textbf{W}^{k}\Vert^2_2=0$, $\lim_{k\rightarrow\infty}\Vert\textbf{b}^{k+1}-\textbf{b}^{k}\Vert^2_2=0$, $\lim_{k\rightarrow\infty}\Vert\textbf{z}^{k+1}-\textbf{z}^{k}\Vert^2_2=0$, $\lim_{k\rightarrow\infty}\Vert\textbf{q}^{k+1}-\textbf{q}^{k}\Vert^2_2=0$, $\lim_{k\rightarrow\infty}\Vert\textbf{u}^{k+1}-\textbf{u}^{k}\Vert^2_2=0$.\\
(2). For the pdADMM-Q algorithm, if $\rho>\max(4\nu S^2,(\sqrt{17}+1)\nu/2)$, then
$\beta_\rho(\textbf{p}^k,\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{q}^{k},\textbf{u}^k)$ is convergent. Moreover, $\lim_{k\rightarrow\infty}\Vert\textbf{W}^{k+1}-\textbf{W}^{k}\Vert^2_2=0$, $\lim_{k\rightarrow\infty}\Vert\textbf{b}^{k+1}-\textbf{b}^{k}\Vert^2_2=0$, $\lim_{k\rightarrow\infty}\Vert\textbf{z}^{k+1}-\textbf{z}^{k}\Vert^2_2=0$, $\lim_{k\rightarrow\infty}\Vert\textbf{q}^{k+1}-\textbf{q}^{k}\Vert^2_2=0$, $\lim_{k\rightarrow\infty}\Vert\textbf{u}^{k+1}-\textbf{u}^{k}\Vert^2_2=0$.
\end{theorem}
\begin{proof}
(1). From Lemmas \ref{lemma:objective decrease} and \ref{lemma:lower bounded}, we know that $L_\rho(\textbf{p}^{k},\textbf{W}^{k},\textbf{b}^k,\textbf{z}^{k},\textbf{q}^k,\textbf{u}^{k})$ is convergent because a monotone bounded sequence converges. Moreover, we take the limit on the both sides of Inequality \eqref{eq:pdadmm objective decrease} to obtain
\begin{align*}
     &0=\lim\nolimits_{k\rightarrow\infty}L_\rho(\textbf{p}^{k},\textbf{W}^{k},\textbf{b}^k,\textbf{z}^{k},\textbf{q}^k,\textbf{u}^{k})\\&-\lim\nolimits_{k\rightarrow\infty}L_\rho(\textbf{p}^{k+1},\textbf{W}^{k+1},\textbf{b}^{k+1},\textbf{z}^{k+1},\textbf{q}^{k+1},\textbf{u}^{k+1})\\ &\geq \lim\nolimits_{k\rightarrow\infty}(\sum\nolimits_{l=2}^L (\tau^{k+1}_l/2)\Vert p^{k+1}_l\!-\!p^k_l\Vert^2_2\\&+\sum\nolimits_{l\!=\!1}^{L}(\theta^{k\!+\!1}_l/2)\Vert  W^{k\!+\!1}_l\!-\!W^k_l\Vert^2_2+\sum\nolimits_{l\!=\!1}^{L}(\nu/2)\Vert  b^{k\!+\!1}_l\!-\!b^k_l\Vert^2_2\\&+\sum\nolimits_{l=1}^{L-1} C_1\Vert z^{k+1}_l-z^k_l\Vert^2_2+(\nu/2)\Vert z^{k+1}_L-z^k_L\Vert^2_2\\&+\sum\nolimits_{l=1}^{L-1}C_2\Vert q^{k+1}_l-q^k_l\Vert^2_2)\geq 0
\end{align*}
Because $L_\rho(\textbf{p}^{k},\textbf{W}^{k},\textbf{b}^k,\textbf{z}^{k},\textbf{q}^k,\textbf{u}^{k})$ is convergent, then $\lim_{k\rightarrow\infty}\Vert\textbf{p}^{k+1}-\textbf{p}^{k}\Vert^2_2=0$, $\lim_{k\rightarrow\infty}\Vert\textbf{W}^{k+1}-\textbf{W}^{k}\Vert^2_2=0$, $\lim_{k\rightarrow\infty}\Vert\textbf{b}^{k+1}-\textbf{b}^{k}\Vert^2_2=0$, $\lim_{k\rightarrow\infty}\Vert\textbf{z}^{k+1}-\textbf{z}^{k}\Vert^2_2=0$, and $\lim_{k\rightarrow\infty}\Vert\textbf{q}^{k+1}-\textbf{q}^{k}\Vert^2_2=0$. $\lim_{k\rightarrow\infty}\Vert\textbf{u}^{k+1}-\textbf{u}^{k}\Vert^2_2=0$ is derived from Lemma \ref{lemma:u bound} in the Appendix.\\
(2). The proof follows the same procedure as (1).
\end{proof}
 \indent The third lemma guarantees that the subgradient of the objective is upper bounded, which is stated as follows:
 \begin{lemma}[Bounded Subgradient]
\label{lemma:subgradient bound}
(1). For the pdADMM algorithm, there exists a constant $C>0$ and $g^{k+1}\in \partial L_\rho(\textbf{p}^{k+1},\textbf{W}^{k+1},\textbf{b}^{k+1},\textbf{z}^{k+1},\textbf{q}^{k+1},\textbf{u}^{k+1})$ such that
\begin{align*}
    &\Vert g^{k+1}\Vert \leq C(\Vert\textbf{p}^{k+1}-\textbf{p}^{k}\Vert+\Vert\textbf{W}^{k+1}-\textbf{W}^{k}\Vert+\Vert\textbf{b}^{k+1}-\textbf{b}^{k}\Vert\\&+\Vert\textbf{z}^{k+1}-\textbf{z}^{k}\Vert+\Vert\textbf{q}^{k+1}-\textbf{q}^{k}\Vert+\Vert\textbf{u}^{k+1}-\textbf{u}^{k}\Vert)
\end{align*}
(2). For the pdADMM-Q algorithm, there exists a constant $\overline{C}>0$, $\overline{g}_\textbf{W}^{k+1}\in \nabla_{\textbf{W}^{k+1}} \beta_\rho(\textbf{p}^{k+1},\textbf{W}^{k+1},\textbf{b}^{k+1},\textbf{z}^{k+1},\textbf{q}^{k+1},\textbf{u}^{k+1})$, $\overline{g}_\textbf{b}^{k+1}\in \nabla_{\textbf{b}^{k+1}} \beta_\rho(\textbf{p}^{k+1},\textbf{W}^{k+1},\textbf{b}^{k+1},\textbf{z}^{k+1},\textbf{q}^{k+1},\textbf{u}^{k+1})$, $\overline{g}_\textbf{z}^{k+1}\in \partial_{\textbf{z}^{k+1}} \beta_\rho(\textbf{p}^{k+1},\textbf{W}^{k+1},\textbf{b}^{k+1},\textbf{z}^{k+1},\textbf{q}^{k+1},\textbf{u}^{k+1})$, $\overline{g}_\textbf{q}^{k+1}\in \nabla_{\textbf{q}^{k+1}} \beta_\rho(\textbf{p}^{k+1},\textbf{W}^{k+1},\textbf{b}^{k+1},\textbf{z}^{k+1},\textbf{q}^{k+1},\textbf{u}^{k+1})$, $\overline{g}_\textbf{u}^{k+1}\in \nabla_{\textbf{u}^{k+1}} \beta_\rho(\textbf{p}^{k+1},\textbf{W}^{k+1},\textbf{b}^{k+1},\textbf{z}^{k+1},\textbf{q}^{k+1},\textbf{u}^{k+1})$ such that
\begin{align*}
    &\Vert \overline{g}_\textbf{W}^{k+1}\Vert \leq \overline{C}(\Vert\textbf{W}^{k+1}-\textbf{W}^{k}\Vert+\Vert\textbf{b}^{k+1}-\textbf{b}^{k}\Vert\\&+\Vert\textbf{z}^{k+1}-\textbf{z}^{k}\Vert+\Vert\textbf{q}^{k+1}-\textbf{q}^{k}\Vert+\Vert\textbf{u}^{k+1}-\textbf{u}^{k}\Vert)\\&
    \Vert \overline{g}_\textbf{b}^{k+1}\Vert \leq \overline{C}(\Vert\textbf{W}^{k+1}-\textbf{W}^{k}\Vert+\Vert\textbf{b}^{k+1}-\textbf{b}^{k}\Vert\\&+\Vert\textbf{z}^{k+1}-\textbf{z}^{k}\Vert+\Vert\textbf{q}^{k+1}-\textbf{q}^{k}\Vert+\Vert\textbf{u}^{k+1}-\textbf{u}^{k}\Vert)\\&\Vert \overline{g}_\textbf{z}^{k+1}\Vert \leq \overline{C}(\Vert\textbf{W}^{k+1}-\textbf{W}^{k}\Vert+\Vert\textbf{b}^{k+1}-\textbf{b}^{k}\Vert\\&+\Vert\textbf{z}^{k+1}-\textbf{z}^{k}\Vert+\Vert\textbf{q}^{k+1}-\textbf{q}^{k}\Vert+\Vert\textbf{u}^{k+1}-\textbf{u}^{k}\Vert)\\&
    \Vert \overline{g}_\textbf{q}^{k+1}\Vert \leq \overline{C}(\Vert\textbf{W}^{k+1}-\textbf{W}^{k}\Vert+\Vert\textbf{b}^{k+1}-\textbf{b}^{k}\Vert\\&+\Vert\textbf{z}^{k+1}-\textbf{z}^{k}\Vert+\Vert\textbf{q}^{k+1}-\textbf{q}^{k}\Vert+\Vert\textbf{u}^{k+1}-\textbf{u}^{k}\Vert)\\&\Vert \overline{g}_\textbf{u}^{k+1}\Vert \leq \overline{C}(\Vert\textbf{W}^{k+1}-\textbf{W}^{k}\Vert+\Vert\textbf{b}^{k+1}-\textbf{b}^{k}\Vert\\&+\Vert\textbf{z}^{k+1}-\textbf{z}^{k}\Vert+\Vert\textbf{q}^{k+1}-\textbf{q}^{k}\Vert+\Vert\textbf{u}^{k+1}-\textbf{u}^{k}\Vert)
\end{align*}
\end{lemma}
 \indent Now based on Theorem \ref{theo:convergent variable}, and Lemma \ref{lemma:subgradient bound}, the convergence of the pdADMM algorithm to a critical point is presented in the following theorem.
\begin{theorem} [Convergence of the pdADMM algorithm]
\label{theo: pdadmm global convergence}
If $\rho>\max(4\nu S^2,(\sqrt{17}+1)\nu/2)$, then for the variables $(\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q},\textbf{u})$ in Problem \ref{prob:problem 2}, starting from any $(\textbf{p}^{0},\textbf{W}^{0},\textbf{b}^{0},\textbf{z}^{0},\textbf{q}^{0},\textbf{u}^{0})$, $(\textbf{p}^{k},\textbf{W}^{k},\textbf{b}^{k},\textbf{z}^{k},\textbf{q}^{k},\textbf{u}^{k})$ has at least a limit point $(\textbf{p}^*,\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{q}^*,\textbf{u}^*)$, and any limit point is a critical point of Problem \ref{prob:problem 2}. That is, $0\in \partial L_\rho(\textbf{p}^*,\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{q}^*,\textbf{u}^*)$. In other words, 
\begin{align*}
    & p^*_{l+1}=q^*_l, \ \nabla_{\textbf{p}^*} L_\rho(\textbf{p}^*,\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{q}^*,\textbf{u}^*)=0,\\& \
    \nabla_{\textbf{W}^*} L_\rho(\textbf{p}^*,\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{q}^*,\textbf{u}^*)=0,\\& \nabla_{\textbf{b}^*} L_\rho(\textbf{p}^*,\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{q}^*,\textbf{u}^*)=0, \\& 0\in\partial_{\textbf{z}^*} L_\rho(\textbf{p}^*,\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{q}^*,\textbf{u}^*),\\& 
     \nabla_{\textbf{q}^*} L_\rho(\textbf{p}^*,\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{q}^*,\textbf{u}^*)=0.
\end{align*}
\end{theorem}
\begin{proof}
 From Lemma \ref{lemma:lower bounded}(1), $(\textbf{p}^k,\textbf{W}^k, \textbf{b}^k, \textbf{z}^k, \textbf{q}^k,\textbf{u}^k)$ has at least a limit point $(\textbf{p}^*,\textbf{W}^*, \textbf{b}^*, \textbf{z}^*, \textbf{q}^*,\textbf{u}^*)$ because a bounded sequence has at least a limit point. From Lemma \ref{lemma:subgradient bound} and Theorem \ref{theo:convergent variable},
$\Vert g^{k+1}\Vert \rightarrow 0$ as $k\rightarrow \infty$. According to the definition of general subgradient (Defintion 8.3 in \cite{rockafellar2009variational}), we have $0\in \partial L_\rho(\textbf{p}^*,\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{q}^*,\textbf{u}^*)$. In other words, every limit point $(\textbf{p}^*,\textbf{W}^*, \textbf{b}^*, \textbf{z}^*, \textbf{q}^*,\textbf{u}^*)$ is a critical point.
\end{proof}
\indent Theorem \ref{theo: pdadmm global convergence} shows that our proposed pdADMM algorithm converges for sufficiently large $\rho$, which is consistent with previous literature \cite{wang2019admm}. Similarly, the convergence of the proposed pdADMM-Q algorithm is shown as follows:
\begin{theorem} [Convergence of the pdADMM-Q algorithm]
\label{theo: pdadmm-q global convergence}
If $\rho>\max(4\nu S^2,(\sqrt{17}+1)\nu/2)$, then for the variables $(\textbf{p},\textbf{W},\textbf{b},\textbf{z},\textbf{q},\textbf{u})$ in Problem \ref{prob:problem 3}, starting from any $(\textbf{p}^{0},\textbf{W}^{0},\textbf{b}^{0},\textbf{z}^{0},\textbf{q}^{0},\textbf{u}^{0})$, $(\textbf{p}^{k},\textbf{W}^{k},\textbf{b}^{k},\textbf{z}^{k},\textbf{q}^{k},\textbf{u}^{k})$ has at least a limit point $(\textbf{p}^*,\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{q}^*,\textbf{u}^*)$, and any limit point $(\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{q}^*,\textbf{u}^*)$ is a critical point of Problem \ref{prob:problem 3}. Moreover, if $\tau^{k+1}_l$ is bounded, then $\textbf{p}^*$ is a quantized critical point of Problem \ref{prob:problem 3}. That is 
\begin{align*}
    & p^*_{l+1}=q^*_l, \
    \nabla_{\textbf{W}^*} \beta_\rho(\textbf{p}^*,\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{q}^*,\textbf{u}^*)=0,\\& \nabla_{\textbf{b}^*} \beta_\rho(\textbf{p}^*,\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{q}^*,\textbf{u}^*)=0, \\& 0\in\partial_{\textbf{z}^*} \beta_\rho(\textbf{p}^*,\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{q}^*,\textbf{u}^*),\\& 
    \nabla_{\textbf{q}^*} \beta_\rho(\textbf{p}^*,\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{q}^*,\textbf{u}^*)=0,\\&
    p^*_l\in\arg\min\nolimits_{\delta\in \Delta}\Vert \delta-(p^*_l-\nabla_{p^*_l} F(\textbf{p}^*,\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{q}^*)/\tau^*_l)\Vert.
\end{align*}
where $\tau^*_l$ is a limit point of $\tau^k_l$.
\end{theorem}
The only difference between Theorems \ref{theo: pdadmm global convergence} and \ref{theo: pdadmm-q global convergence} is that $\textbf{p}^*$ is a critical point in Problem \ref{prob:problem 2} and a quantized critical point in Problem \ref{prob:problem 3}.  Next, the following theorem ensures the sublinear convergence rate $o(1/k)$ of the proposed pdADMM algorithm and the pdADMM-Q algorithm.
\begin{theorem}[Convergence Rate]
(1). For the pdADMM algorithm and a sequence $(\textbf{p}^k,\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{q}^k,\textbf{u}^k)$, define $c_k=\min\nolimits_{0\leq i\leq k}(\sum\nolimits_{l=2}^L (\tau^{i+1}_l/2)\Vert p^{i+1}_l\!-\!p^i_l\Vert^2_2+\sum\nolimits_{l\!=\!1}^{L}(\theta^{i\!+\!1}_l/2)\Vert  W^{i\!+\!1}_l\!-\!W^i_l\Vert^2_2+\sum\nolimits_{l\!=\!1}^{L}(\nu/2)\Vert  b^{i\!+\!1}_l\!-\!b^i_l\Vert^2_2+\sum\nolimits_{l=1}^{L-1} C_1\Vert z^{i+1}_l-z^i_l\Vert^2_2+(\nu/2)\Vert z^{i+1}_L-z^i_L\Vert^2_2+\sum\nolimits_{l=1}^{L-1}C_2\Vert q^{i+1}_l-q^i_l\Vert^2_2)$ where $C_1=\nu/2-2\nu^2S^2/\rho>0$ and $C_2=\rho/2-2\nu^2/\rho-\nu/2>0$, then the convergence rate of $c_k$ is $o(1/k)$.\\
(2). For the pdADMM-Q algorithm and a sequence $(\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{q}^k,\textbf{u}^k)$, define $d_k=\min\nolimits_{0\leq i\leq k}(\sum\nolimits_{l\!=\!1}^{L}(\theta^{i\!+\!1}_l/2)\Vert  W^{i\!+\!1}_l\!-\!W^i_l\Vert^2_2+\sum\nolimits_{l\!=\!1}^{L}(\nu/2)\Vert  b^{i\!+\!1}_l\!-\!b^i_l\Vert^2_2+\sum\nolimits_{l=1}^{L-1} C_1\Vert z^{i+1}_l-z^i_l\Vert^2_2+(\nu/2)\Vert z^{i+1}_L-z^i_L\Vert^2_2+\sum\nolimits_{l=1}^{L-1}C_2\Vert q^{i+1}_l-q^i_l\Vert^2_2)$ where $C_1=\nu/2-2\nu^2S^2/\rho>0$ and $C_2=\rho/2-2\nu^2/\rho-\nu/2>0$, then the convergence rate of $d_k$ is $o(1/k)$.
\label{theo: convergence rate}
\end{theorem}
\section{Experiments}
\label{sec:experiments}
\indent In this section, we evaluate the performance of the pdADMM on GA-MLP models using six benchmark datasets.  Speedup and accuracy are compared with several state-of-the-art optimizers. All experiments were conducted on a 64-bit machine with Intel Xeon(R) silver 4114 Processor and 48GB RAM. 
\subsection{Datasets}
 In this experiment, six benchmark datasets were used for performance evaluation, which is shown in Table \ref{tab:dataset}, each of which is detailed below:\\
 1. Cora \cite{sen2008collective}. The Cora dataset consists of 2708 scientific publications classified into one of seven classes. The citation network consists of 5429 links. Each publication in the dataset is described by a 0/1-valued word vector indicating the absence/presence of the corresponding word from the dictionary. The dictionary consists of 1433 unique words.\\
2. PubMed \cite{sen2008collective}. PubMed comprises 30M+ citations for biomedical literature that have been collected from sources such as MEDLINE, life science journals, and published online e-books. It also includes links to text content from PubMed Central and other publishers' websites.\\
3. Citeseer \cite{sen2008collective}. The Citeseer dataset was collected from the Tagged.com social network website. It contains 5.6 million users and 858 million links between them. Each user has 4 features and is manually labeled as ``spammer" or ``not spammer". Each link represents an action between two users and includes a timestamp and a type. The network contains 7 anonymized types of links. The original task on the dataset is to identify (i.e., classify) the spammer users based on their relational and non-relational features. \\
4. Amazon Computers and Amazon Photo \cite{mcauley2015image}. Amazon Computers and Amazon Photo are segments of the Amazon co-purchase graph, where nodes represent goods, edges indicate that two goods are frequently bought
together, node features are bag-of-words encoded product reviews, and class labels are given by the
product category.\\
5. Coauthor CS \cite{shchur2018pitfalls}. Coauthor CS is a co-authorship graph based on the Microsoft Academic Graph
from the KDD Cup 2016 challenge 3. Here, nodes are authors, that are connected by an edge if they
co-authored a paper; node features represent paper keywords for each author’s papers, and class
labels indicate the most active fields of study for each author.
\begin{table}[]
    \small
    \centering
    \begin{tabular}{p{1.2cm}|p{0.8cm}|p{1.1cm}|p{0.9cm}|p{0.8cm}|p{0.8cm}}
    \hline\hline
         Dataset&\tabincell{c}{ Node\\\#}&\tabincell{c}{Training\\ Sample\\\#}&\tabincell{c}{Test\\ Sample\\\#}&\tabincell{c}{ Class\\\#}&
\tabincell{c}{Feature\\\#}\\\hline Cora&2708&140&1000&7&1433\\\hline PubMed&19717&60&1000&3&500\\\hline Citeseer&2110&120&1000&6&3703\\\hline 
\tabincell{c}{Amazon\\  Computers}&13381&200&1000&10&767\\\hline 
\tabincell{c}{Amazon\\ Photo}&7487&160&1000&8&745\\\hline 
\tabincell{c}{Coauthor\\ CS} &54493&100&1000&5&8415\\\hline\hline
    \end{tabular}
    \caption{Six benchmark datasets}
    \label{tab:dataset}
\end{table}


\subsection{Comparison Methods}
\indent SGD and its variants are state-of-the-art methods and hence were served as comparison methods. For SGD-based methods, the full batch dataset is used for training models. All parameters were chosen by maximizing the accuracy of training datasets. The baselines are described as follows: \\
\indent 1. Gradient Descent (GD) \cite{bottou2010large}. The GD and its variants are the most popular deep learning optimizers, whose convergence has been studied extensively in the literature. The learning rate of GD was set to $10^{-2}$.\\
\indent 2. Adaptive learning rate method (Adadelta) \cite{zeiler2012adadelta}. The Adadelta is proposed to overcome the sensitivity to hyperparameter selection. The learning rate of Adadelta was set to $10^{-4}$.\\
\indent 3. Adaptive gradient algorithm (Adagrad) \cite{duchi2011adaptive}. Adagrad is an improved version of SGD: rather than fixing the learning rate during training, it adapts the learning rate to the hyperparameter. The learning rate of Adagrad was set to $10^{-3}$.\\
\indent 4. Adaptive momentum estimation (Adam) \cite{kingma2014adam}. Adam is the most popular optimization method for deep learning models. It estimates the first and second momentum to correct the biased gradient and thus makes convergence fast. The learning rate of Adam was set to $10^{-3}$.


\subsection{Speedup}
  \begin{figure}[]
  \centering
  \centerline{\includegraphics[width=0.9\linewidth]{speedup.pdf}}
  \caption{The relationship between speedup and: (a) the number of layers; (b) the number of neurons in six datasets: the speedup increases linearly with the number of layers and the number of neurons.}
  \label{fig:speedup}
\end{figure}

\indent Firstly, we investigate the speedup of the pdADMM  concerning the number of layers and the number of neurons in the large-scale deep GA-MLP models. The activation function was set to the Rectified linear unit (Relu). The loss function was the cross-entropy loss. The running time per epoch was the average of 10 epochs. $\rho$ and $\nu$ were both set to $10^{-3}$.\\
 \indent For the relationship between speedup and the number of neurons, the pdADMM in the GA-MLP models with 10 layers was tested. The number of neurons in each layer ranged from 3,000 to 4,200. The speedup is shown in Figure \ref{fig:speedup}(a) in six datasets.  The speedup increases linearly with the number of neurons. Moreover, the speedup in the Coauthor CS dataset is more significant than others: it can achieve roughly 4 even for 3000 neurons, while the best speedup of others is about 3.5 for 4,200 neurons.\\
  \indent Now the relationship between speedup and the number of layers is investigated. We set up a GA-MLP model with different numbers of hidden layers, which ranged from 8 to 17. The number of neurons in each layer was fixed to 2,000.  Figure \ref{fig:speedup}(b) shows that the speedup scales linearly with the number of layers. Again, the speedup in the Coauthor CS dataset is the most significant: it reaches a speedup between 5 and 6, whereas the speedup of PubMed is only 2. The other four datasets achieve a speedup between 3 and 4.  We conclude by Figure \ref{fig:speedup} that the number of neurons has a more significant effect on speedup than the number of layers: for example, in the Cora dataset, the speedup increases approximately by $50\%$ when the growth of the number of neurons is about $40\%$; However, while the number of layers doubles, the speedup grows only about $30\%$.
 \subsection{Performance}
  \begin{figure*}
    \begin{minipage}{0.3\linewidth}
    \includegraphics[width=\linewidth]{cora_train.pdf}
    \centerline{(a). Cora}
    \end{minipage}
    \hfill
    \begin{minipage}{0.3\linewidth}
    \includegraphics[width=\columnwidth]{pubmed_train.pdf}
    \centerline{(b). PubMed}
    \end{minipage}
    \hfill
    \begin{minipage}{0.3\linewidth}
    \includegraphics[width=\linewidth]{citeseer_train.pdf}
        \centerline{(c). Citeseer}
    \end{minipage}
    \vfill
    \begin{minipage}{0.3\linewidth}
    \includegraphics[width=\linewidth]{amazon_computers_train.pdf}
    \centerline{(d). Amazon Computers}
    \end{minipage}
    \hfill
    \begin{minipage}{0.3\linewidth}
    \includegraphics[width=\linewidth]{amazon_photo_train.pdf}
        \centerline{(e). Amazon Photo}
    \end{minipage}
    \hfill
    \begin{minipage}{0.3\linewidth}
    \includegraphics[width=\linewidth]{coauthor_cs_train.pdf}
            \centerline{(f). Coauthor CS}
    \end{minipage}
    \caption{Training Accuracy of all methods: all methods achieve nearly $100\%$ accuracy in every dataset.}
    \label{fig:train performance}
\end{figure*}
\begin{figure*}
    \begin{minipage}{0.3\linewidth}
    \includegraphics[width=\linewidth]{cora_test.pdf}
    \centerline{(a). Cora}
    \end{minipage}
    \hfill
    \begin{minipage}{0.3\linewidth}
    \includegraphics[width=\columnwidth]{pubmed_test.pdf}
    \centerline{(b). PubMed}
    \end{minipage}
    \hfill
    \begin{minipage}{0.3\linewidth}
    \includegraphics[width=\linewidth]{citeseer_test.pdf}
        \centerline{(c). Citeseer}
    \end{minipage}
    \vfill
    \begin{minipage}{0.3\linewidth}
    \includegraphics[width=\linewidth]{amazon_computers_test.pdf}
    \centerline{(d). Amazon Computers}
    \end{minipage}
    \hfill
    \begin{minipage}{0.3\linewidth}
    \includegraphics[width=\linewidth]{amazon_photo_test.pdf}
        \centerline{(e). Amazon Photo}
    \end{minipage}
    \hfill
    \begin{minipage}{0.3\linewidth}
    \includegraphics[width=\linewidth]{coauthor_cs_test.pdf}
            \centerline{(f). Coauthor CS}
    \end{minipage}
    \caption{Test Accuracy of all methods: the pdADMM outperforms other comparison methods in every dataset.}
    \label{fig:test performance}
\end{figure*}
\indent Finally, the performance of the pdADMM is analyzed against comparison methods. All methods were tested on a GA-MLP model with $10$ layers. The number of neurons was fixed to $2000$. Our pdADMM was firstly trained in a two-layer GA-MLP model, and then the layer number of the GA-MLP model increased to five, and finally, all layers were involved in training. Figures \ref{fig:train performance} and \ref{fig:test performance} show the training accuracy and the test accuracy of the pdADMM against all comparison methods in six datasets, respectively. X-axis and Y-axis represent epoch and accuracy, respectively. Overall, all methods perform well in the training sets: they reach nearly $100\%$ accruacy in every datasets. However, the pdADMM outperforms others marginally in the test set. For example, the pdADMM accomplishes $80\%$ test accuracy in the cora datasets, while the performance of other methods is all below $80\%$. In the Coauthor CS dataset, the pdADMM reaches $85\%$ test accuracy, about $5\%$ better than the best comparison method Adam.
\section{Conclusion}
\label{sec:conclusion}
The Graph Augmented Multi-Layer Perceptron (GA-MLP) is an attractive model due to resistance to the over-smoothing problem. In this paper, we propose the pdADMM algorithm and the pdADMM-Q algorithm to achieve model parallelism. Their theoretical convergence to a critical point is guaranteed with a sublinear convergence $o(1/k)$. Extensive experiments verify that the pdADMM algorithm not only accelerates the training of deep GA-MLP models but also stands out among all the existing state-of-the-art comparison methods in six benchmark datasets.
\bibliographystyle{plain}
\small
