\section{Introduction}
\vspace*{-6pt}
Alternating direction method of multipliers (ADMM) \cite{gabay1976dual,boyd2011distributed} is a powerful optimization tool
for the composite or constrained problems in machine learning. In general, it
considers the following optimization problem:
\vspace*{-4pt}
\begin{equation} \label{eq:1}
\min_{x,y} f(x) + g(y), \quad \mbox{s.t.} \ Ax + By =c, \nonumber
\vspace*{-4pt}
\end{equation}
where $f(x):\mathbb{R}^d\rightarrow \mathbb{R}$ and $g(y):\mathbb{R}^p\rightarrow \mathbb{R}$ are convex functions.
For example, in machine learning, $f(x)$ can be used for the empirical loss, $g(y)$ for the structure regularizer,
and the constraint for encoding the structure pattern of model parameters.
Due to the flexibility in splitting the objective function into loss $f(x)$ and regularizer $g(y)$,
the ADMM can relatively easily solve some complicated structure problems in machine learning,
such as the graph-guided fused lasso \cite{kim2009multivariate} and the overlapping group lasso,
which are too complicated for the other popular optimization methods
such as proximal gradient methods \cite{nesterov2005smooth,beck2009fast}.
Thus, the ADMM has been extensively studied in recent years \cite{boyd2011distributed,nishihara2015general,xu2017admm}.

\vspace*{-2pt}
The above deterministic ADMM generally needs to compute the gradients of empirical loss function
on all examples at each iteration,
which makes it unsuitable for solving big data problems.
Thus, the online and stochastic versions of ADMM \cite{wang2012online,suzuki2013dual,ouyang2013stochastic}
are developed. However, due to large variance of stochastic gradients,
these stochastic methods suffer from a slow convergence rate.
Recently, some fast stochastic ADMM methods \cite{zhong2014fast,suzuki2014stochastic,zheng2016fast} have been proposed
by using the variance reduced (VR) techniques.

\begin{table*}
\vspace*{-0.4cm}
  \centering
  \caption{  IFO complexity comparison of the \textbf{non-convex ADMM} methods for finding
  an $\epsilon$-approximate stationary point of the problem \eqref{eq:2},
  i.e., $\mathbb{E}\|\nabla \mathcal {L}(x,y_{[m]},z)\|^2\leq \epsilon$. $n$ denotes the sample size.}
  \label{tab:1}
\begin{tabular}{c|c|c|c}
  \hline
 
  \textbf{Problem} & \textbf{Algorithm} & \textbf{Reference}  & \textbf{IFO}  \\ \hline
  \multirow{4}*{Finite-sum} & ADMM & \citet{jiang2019structured}    &  $\mathcal{O}(n\epsilon^{-1})$ \\ \cline{2-4}
  &SVRG-ADMM & \citet{huang2016stochastic,zheng2016stochastic}    &  $\mathcal{O}(n+n^{\frac{2}{3}}\epsilon^{-1})$ \\ \cline{2-4}
  &SAGA-ADMM & \citet{huang2016stochastic}    & $\mathcal{O}(n+n^{\frac{2}{3}}\epsilon^{-1})$ \\ \cline{2-4}
  &SPIDER-ADMM  & Ours  & {\color{red}{$\mathcal{O}( n+n^{\frac{1}{2}}\epsilon^{-1})$ }}\\ \hline
  \multirow{2}*{Online} & SADMM  &\citet{huang2018mini}   & $\mathcal{O}(\epsilon^{-2})$ \\ \cline{2-4}
  &Online SPIDER-ADMM  & Ours  & {\color{red}{$\mathcal{O}(\epsilon^{-\frac{3}{2}})$}}  \\
  \hline
\end{tabular}
\vspace*{-0.5cm}
\end{table*}

\vspace*{-2pt}
So far, the above discussed ADMM methods build on the convexity of objective functions.
In fact, ADMM is also highly successful in solving various nonconvex problems
such as tensor decomposition \cite{kolda2009tensor} and training neural networks \cite{Taylor2016Training}.
Thus, some works \cite{li2015global,wang2015convergence,wang2015global,hong2016convergence,jiang2019structured}
have devoted to studying the non-convex ADMM methods.
More recently, for solving the big data problems,
the nonconvex stochastic ADMMs \cite{huang2016stochastic,zheng2016stochastic} have been
proposed with the VR techniques such as the SVRG \cite{johnson2013accelerating} and
the SAGA \cite{defazio2014saga}. In addition, \citet{huang2018mini}
have extended the online/stochastic ADMM \cite{ouyang2013stochastic}
to the nonconvex setting.

\vspace*{-2pt}
Although these works have studied the convergence of nonconvex stochastic ADMMs and proved
these methods have $\mathcal{O}(\frac{c}{T})$ convergence rate, where $T$ denotes number of iteration and
$c$ a constant independent on $T$, they have not provided the \textbf{optimal} incremental/stochastic first-order oracle (IFO/SFO \cite{Ghadimi2013StochasticFA}) complexity
for these methods yet. In other words, they have only proved these stochastic ADMMs have the same convergence rate to
the deterministic ADMM \cite{jiang2019structured}, but don't tell us whether these stochastic ADMMs have less IFO complexity than the deterministic ADMM,
which is a key assessment criteria of the first-order stochastic methods \cite{Reddi2016Prox}.
For example, from the existing noncovex SAGA-ADMM and SVRG-ADMM \cite{zheng2016stochastic,huang2016stochastic},
we only obtain a \textbf{rough} IFO complexity of $\mathcal{O}(n + bc\epsilon^{-1})$ for finding an $\epsilon$-approximate stationary point,
where $b$ denotes the mini-batch size. In their convergence analysis, to ensure the convergence of these methods, they need to choose a small step size $\eta$
and a large penalty parameter $\rho$.
Under this case, we maybe have $bc\geq n$, so that these stochastic ADMMs have no less IFO complexity than the deterministic ADMM.
Thus, there still exist two important problems to be addressed:
\vspace*{-4pt}
\begin{itemize}
\vspace*{-6pt}\item \emph{ Does the stochastic ADMM have less IFO complexity than the deterministic ADMM for nonconvex optimization?}
\vspace*{-6pt}\item \emph{ If the stochastic ADMM improves IFO complexity, how much can it improve?}
\end{itemize}
\vspace*{-10pt}
In the paper,  we answer the above challenging questions
with positive solutions and propose a new faster stochastic
ADMM method (\emph{i.e.,} SPIDER-ADMM) to solve the following nonconvex nonsmooth problem:
\vspace*{-3pt}
\begin{align} \label{eq:2}
\min_{x,\{y_j\}_{j=1}^m} & f(x) := \left\{
\begin{aligned}
 & \frac{1}{n}\sum_{i=1}^n f_i(x)  \ \mbox{(finite-sum)}\\
 & \mathbb{E}_{\zeta} [f(x;\zeta)]  \ \mbox{(online)}
 \end{aligned} \right.  +  \sum_{j=1}^m g_j(y_j) \nonumber \\
 \mbox{s.t.} & \quad Ax + \sum_{j=1}^mB_jy_j =c,
\end{align}
where $A\in \mathbb{R}^{l\times d}$, $B_j\in \mathbb{R}^{l\times p}$ for all $j\in [m]$, and 
$\zeta$ is a random variable following an unknown distribution.
Here $f(x): \mathbb{R}^d\rightarrow \mathbb{R}$ is a \emph{nonconvex} and smooth
function, and $g_j(y_j): \mathbb{R}^p\rightarrow \mathbb{R}$ is a convex and possibly
\emph{nonsmooth} function for all $j\in [m], \ m\geq 1$.
In machine learning, $f(x)$ can be used for losses such as activation functions of neural networks, $\sum_{j=1}^m g_j(y_j)$ can
be used for not only single structure penalty (\emph{e.g.,} sparse, low rank) but also
superposition structures penalties (\emph{e.g.}, sparse + low rank, sparse + group sparse), which are
widely applied in robust PCA \cite{candes2011robust}, subspace clustering \cite{liu2010robust}, and dirty models \cite{jalali2010dirty}.
For the problem \eqref{eq:2}, its finite-sum subproblem generally arises from the empirical loss minimization and M-estimation.
While its online subproblem comes from the expected loss minimization. To address the online subproblem, we extend the SPIDER-ADMM to
the online setting, and propose an online SPIDER-ADMM.

\subsection*{Challenges and Contributions}
Our SPIDER-ADMM methods build on the variance-reduced technique of SPIDER \cite{fang2018spider} and SpiderBoost \cite{wang2018spiderboost},
which is a variant of stochastic
recursive gradient algorithm (SARAH \cite{nguyen2017sarah,nguyen2017stochastic}) and reaches the state-of-the-art IFO complexity
as the SNVRG \cite{zhou2018stochastic}.
Although the SPIDER and SpiderBoost have shown good performances in the stochastic gradient descent (SGD) and proximal SGD methods,
applying these techniques to the nonconvex ADMM method \emph{is not a trivial task}.
There exist the following two main \textbf{challenges}:
\vspace{-6pt}
\begin{itemize}
 \vspace*{-3pt} \item Due to failure of the Fej\'{e}r monotonicity of iteration, the convergence analysis of the nonconvex ADMM is generally quite difficult \cite{wang2015convergence}.
                 With using the inexact stochastic gradient, this difficulty is greater in the nonconvex stochastic ADMM methods;
 \vspace*{-10pt} \item To obtain the optimal IFO complexity of our methods, we need to design a new effective \emph{Lyapunov} function, which can not follow
                 the existing nonconvex stochastic ADMM methods \cite{huang2016stochastic}.
 \vspace*{-10pt}
\end{itemize}
In this paper, thus, we will fill this gap between
the nonconvex ADMM and the SPIDER/SpiderBoost methods.
Our main \textbf{contributions} are summarized as follows:
\vspace{-6pt}
\begin{itemize}
\setlength\itemsep{0em}
\vspace*{-4pt}
\item[1)] We propose a faster stochastic ADMM ( \emph{i.e.,} SPIDER-ADMM ) method for nonconvex optimization based on the SPIDER/SpiderBoost.
          Moreover, we prove that the SPIDER-ADMM achieves a lower
          IFO complexity of $\mathcal{O}(n+n^{1/2}\epsilon^{-1})$
          for finding an $\epsilon$-approximate stationary point,
          which improves the deterministic ADMM by a factor $\mathcal{O}(n^{1/2})$.
 \vspace*{-4pt}
\item[2)] We extend the SPIDER-ADMM method to the online setting, and propose a faster online SPIDER-ADMM for
          nonconvex optimization. Moreover, we prove that the online SPIDER-ADMM achieves a lower
          IFO complexity of $\mathcal{O}(\epsilon^{-\frac{3}{2}})$,
          which improves the existing best results by a factor of $\mathcal{O}(\epsilon^{-\frac{1}{2}})$.
\vspace*{-4pt}
\item[3)] We provide a useful theoretical analysis framework for nonconvex stochastic ADMM methods
          with providing the optimal IFO complexity.
          Based on our new analysis framework, we also prove that the existing nonconvex SVRG-ADMM and SAGA-ADMM
          have the optimal IFO complexity of $\mathcal{O}(n+n^{2/3}\epsilon^{-1})$.
          Thus, our SPIDER-ADMM improves the existing stochastic ADMMs by a factor of $\mathcal{O}(n^{1/6})$.
\end{itemize}
\subsection*{Notations}
 \vspace*{-6pt}
Let $y_{[m]}= \{y_1,\cdots,y_m\}$ and $y_{[j:m]}= \{y_{j},\cdots,y_m\}$ for $j\in[m]=\{1,2,\cdots,m\}$.
Given a positive definite matrix $G$, $\|x\|^2_G = x^TGx$;
$\sigma_{\max}(G)$ and $\sigma_{\min}(G)$
denote the largest and smallest eigenvalues of matrix $G$, respectively;
$\kappa_G = \frac{\sigma_{\max}(G)}{\sigma_{\min}(G)} \geq 1$.
$\sigma^A_{\max}$ and $\sigma^A_{\min}$ denote the largest and smallest eigenvalues of matrix $A^T A$, respectively.
Given positive definite matrices $\{H_j\}_{j=1}^m$, let $\sigma^H_{\min} = \min_j \sigma_{\min}(H_j)$
and $\sigma^H_{\max} = \max_j \sigma_{\max}(H_j)$. $I_d$ denotes a $d\times d$ identity matrix.
\section{Preliminaries}
 \vspace*{-6pt}
In the section, we introduce some preliminaries regarding problem \eqref{eq:2}.
First, we restate the standard
$\epsilon$-approximate stationary point of the nonconvex problem \eqref{eq:2}
used in \cite{jiang2019structured,zheng2016stochastic}.

\begin{definition} \label{def:1}
Given $\epsilon>0$, the point $(x^*,y_{[m]}^*,z^*)$ is said to be an
$\epsilon$-stationary point of the problem \eqref{eq:2},
if it holds that
\begin{align}
 \mathbb{E}\big[ \mbox{dist}(0,\partial L(x^*,y^*_{[m]},z^*))^2 \big] \leq \epsilon,
\end{align}
where $L(x,y_{[m]},z)=f(x) + \sum_{j=1}^m g_j(y_j) - \langle z, Ax + \sum_{j=1}^mB_jy_j-c \rangle$,
\vspace*{-6pt}
\begin{align}
   \partial L(x,y_{[m]},z) = \left [ \begin{matrix}
     \nabla_x L(x,y_{[m]},z) \\
     \partial_{y_1} L(x,y_{[m]},z) \\
      \cdots \\
     \partial_{y_m} L(x,y_{[m]},z) \\
     -Ax-\sum_{j=1}^mB_jy_j+c
 \end{matrix}
 \right ], \nonumber
\end{align}
\vspace*{-6pt}
and $\mbox{dist}(0,\partial L)=\inf_{L'\in \partial L} \|0-L'\|.$
\end{definition}

Next, we give some standard
assumptions regarding problem \eqref{eq:2} as follows:
\begin{assumption}
Each loss function $f_i(x)$ is $L$-smooth such that
\vspace*{-6pt}
\begin{align}
\|\nabla f_i(x)-\nabla f_i(y)\| \leq L \|x - y\|, \ \forall x,y \in \mathbb{R}^d, \nonumber
\end{align}
which is equivalent to
\vspace*{-6pt}
\begin{align}
f_i(x) \leq f_i(y) + \nabla f_i(y)^T(x-y) + \frac{L}{2}\|x-y\|^2.  \nonumber
\end{align}
\end{assumption}
\vspace*{-6pt}
\begin{assumption}
Full gradient of loss function $f(x)$ is bounded, i.e., there exists a constant $\delta >0$ such that for all $x$,
it follows $\|\nabla f(x)\|^2 \leq \delta^2$.
\end{assumption}
\vspace*{-6pt}
\begin{assumption}
$f(x) $ and $g_j(y_j)$ for all $j\in [m]$ are all lower bounded, and let
$f^*=\inf_x f(x) > - \infty$ and $g_j^*=\inf_{y_j} g_j(y_j) > - \infty$.
\end{assumption}
\vspace*{-6pt}
\begin{assumption}
$A$ is a full row or column rank matrix.
\end{assumption}

Assumption 1 imposes smoothness on the individual loss functions, which is commonly used in convergence analysis of
the nonconvex algorithms \cite{Ghadimi2013StochasticFA,ghadimi2016mini}.
Assumption 2 shows full gradient of loss function have a bounded norm, which is used in
the stochastic gradient-based and ADMM-type methods \cite{boyd2011distributed,suzuki2013dual,hazan2016introduction}.
Assumptions 3 and 4 have been used in the study of nonconvex ADMMs
\cite{hong2016convergence,jiang2019structured,zheng2016stochastic}.
Assumptions 3 guarantees the feasibility of the problem \eqref{eq:2}.
Assumption 4 guarantees the matrix $A^TA$ or $AA^T$ is non-singular.
Since there exist multiple regularizers in the above problem \eqref{eq:2},
$A$ is general a full column rank matrix.
Without loss of generality, we will use the full column rank matrix $A$ below.
\section{ Fast SPIDER-ADMM Method }
In the section, we propose a new faster stochastic ADMM algorithm, \emph{i.e.,} SPIDER-ADMM, to solve
the finite-sum problem \eqref{eq:2}.
We begin with giving the augmented Lagrangian function of the problem \eqref{eq:2}:
\vspace*{-6pt}
\begin{align}
 \mathcal {L}_{\rho}(x,y_{[m]},z) = & f(x) + \sum_{j=1}^mg_j(y_j) - \langle z, Ax + \sum_{j=1}^m B_jy_j-c\rangle \nonumber \\
 &  + \frac{\rho}{2} \|Ax+ \sum_{j=1}^mB_jy_j-c\|^2,
\end{align}
where $z \in \mathbb{R}^l$ and $\rho >0$ denote the dual variable and penalty parameter, respectively.
Algorithm \ref{alg:1} gives the SPIDER-ADMM algorithmic framework.

In Algorithm \ref{alg:1}, we use the proximal method to update the variables $\{y_j\}_{j=1}^m$.
At the step 9 of Algorithm \ref{alg:1}, we update the variables $\{y_j\}_{j=1}^m$ by solving the following subproblem, for all $j\in [m]$
\begin{align}
y^{k+1}_j \!=\! \mathop{\arg\min}_{y_j\in \mathbb{R}^p} \mathcal {L}_{\rho}(x_k,y^{k+1}_{[j-1]},y_j,y^{k}_{[j+1:m]},z_k)
\!+\! \frac{1}{2}\|y_j\!-\!y^k_j\|_{H_j}^2, \nonumber
\end{align}
where
\begin{align}
\vspace*{-6pt}
& \mathcal {L}_{\rho} (x_k,y^{k+1}_{[j-1]},y_j,y^{k}_{[j+1:m]},z_k) =  f(x_k) + \sum_{i=1}^{j-1} g_i(y^{k+1}_{i}) \nonumber  \\
&+ g_j(y_{j}) \!+\! \sum_{i=j+1}^m g_j(y^{k}_{j}) \!-\! z_k^T( B_jy_j +\tilde{c}) \!+\! \frac{\rho}{2}\| B_jy_j +\tilde{c} \|^2 \nonumber
\vspace*{-6pt}
\end{align}
with $\tilde{c} = Ax_k+\sum_{i=1}^{j-1}B_iy_i^{k+1}+ \sum_{i=j+1}^{m}B_iy_i^{k}-c$ and $H_j\succ 0$.
When set $H_j = r_j I_p - \rho B_j^TB_j\succeq I_p$ with $r_j \geq \rho
\sigma_{\max}(B^T_j B_j) + 1$ for all $j\in [m]$ to linearize the term $\frac{\rho}{2}\|B_jy_j + \tilde{c}\|^2$,
then we can use the following proximal operator to update $y_j$, for all $j\in [m]$
\begin{align}
 y^{k+1}_j = \mathop{\arg\min}_{y_j\in \mathbb{R}^p} \frac{r_j}{2}\|y_j-w^k_j\|^2 +  g_j(y_j),
\end{align}
where $w^k_j=\frac{1}{r_j}\big( H_jy^k_j - \rho B^T_j\tilde{c} + B^T_jz_k\big)$.

To update $x$, we define
an approximated function over $x_k$ as follows:
\vspace*{-6pt}
\begin{align}
 &\hat{\mathcal {L}}_{\rho} (x,y^{k+1}_{[m]}, z_k,v_k) =
 f(x_k) + v_k^T(x-x_k)+ \frac{1}{2\eta}\|x-x_k\|^2_G \nonumber \\
 & + \sum_{j=1}^m g_j(y^{k+1}_{j}) - z_k^T(Ax+ \sum_{j=1}^m B_jy^{k+1}_j-c) \nonumber \\
 & + \frac{\rho}{2}\|Ax+\sum_{j=1}^mB_jy_j^{k+1}-c\|^2,
\end{align}
where $\eta>0$ is a step size; $v_k$ is a stochastic gradient over $x_k$;
$G\succ 0$ is a positive matrix.
In updating $x$, to avoid computing inverse of $\frac{G}{\eta} + \rho A^TA$,
we can set $G = r I_d - \rho \eta A^TA \succeq I_d$ with $r \geq \rho \eta
\sigma^A_{\max} + 1 $ to linearize term $\frac{\rho}{2}
\|Ax+ \sum_{j=1}^mB_jy^{k+1}_{j}-c\|^2$. Then at the step 10 of Algorithm \ref{alg:1}, we have
\begin{align}
 x_{k+1} = \frac{Gx_k}{r} - \frac{\eta v_k}{r} - \frac{\eta\rho}{r}A^T(\sum_{j=1}^m B_jy_j^{k+1} - c - \frac{z_k}{\rho}). \nonumber
\end{align}

In Algorithm \ref{alg:1}, after setting $v_0 = \nabla f(x_0)$, for each subsequent iteration $k$, we have:
\vspace*{-6pt}
\begin{align}
 v_{k} = \nabla f_{\mathcal{I}_k}(x_k) - \nabla f_{\mathcal{I}_k}(x_{k-1}) + v_{k-1},
\end{align}
where $\nabla f_{\mathcal{I}_k}(x_{k}) = \frac{1}{|\mathcal{I}_k|} \sum_{i\in\mathcal{I}_k }\nabla f_i(x_t)$.
It is easy to check $\mathbb{E}[v_k|x_0] = \nabla f(x_k)$,
\emph{i.e.,} an unbiased estimate gradient over $x_k$.
Comparing the existing SVRG-ADMM,
our SPIDER-ADMM constructs stochastic gradient $v_k$ based on
the information $x_{k-1}$ and $v_{k-1}$, while the SVRG-ADMM constructs $v_k$ based on
the information $x_{0}$ and $v_{0}$ (\emph{i.e.,} the initalization information of each outer loop).
Due to using more fresh information, thus, SPIDER-ADMM
can yield more accurate estimation of the full gradient than SVRG-ADMM.
Simultaneously, it does not require to additional computation and memory,
so it costs less memory than the existing SAGA-ADMM.
\begin{algorithm}[htb]
   \caption{ SPIDER-ADMM Algorithm }
   \label{alg:1}
\begin{algorithmic}[1]
   \STATE {\bfseries Input:} $b$, $q$, $K$, $\eta>0$ and $\rho>0 $;
   \STATE {\bfseries Initialize:} $x_0 \in \mathbb{R}^d$, $y^0_j \in \mathbb{R}^p, \ j\in [m]$ and $z_0 \in \mathbb{R}^l$;
   \FOR {$k=0,1,\cdots,K-1$}
   \IF {$\mbox{mod}(k,q) = 0$}
   \STATE{} Compute $v_k = \nabla f(x_k)$;
   \ELSE
   \STATE{} Uniformly randomly pick a mini-batch $\mathcal{I}_k$ (with replacement) from $\{1,2,\cdots,n\}$ with $|\mathcal{I}_k|=b$,
            and compute $$ v_{k} = \nabla f_{\mathcal{I}_k}(x_k) - \nabla f_{\mathcal{I}_k}(x_{k-1}) + v_{k-1};$$
   \ENDIF
   \STATE{}  $ y^{k+1}_j= \arg\min_{y_j} \big\{ \mathcal {L}_{\rho}(x_k,y^{k+1}_{[j-1]},y_j,y^{k}_{[j+1:m]},z_k) + \frac{1}{2}\|y_j-y^k_j\|_{H_j}^2 \big\}$ for all $j\in [m]$;
   \STATE{}  $ x_{k+1}= \arg\min_x \hat{\mathcal {L}}_{\rho}\big( x,y^{k+1}_{[m]},z_k,v_k \big)$;
   \STATE{}  $z_{k+1} = z_k- \rho(Ax_{k+1} + \sum_{j=1}^mB_jy^{k+1}_j - c)$;
   \ENDFOR
   \STATE {\bfseries Output \ (in theory):} Chosen uniformly random from $\{x_{k},y_{[m]}^{k},z_k\}_{k=1}^{K}$.
   \STATE {\bfseries Output \ (in practice):} $\{x_{K},y_{[m]}^{K},z_K\}$.
\end{algorithmic}
\end{algorithm}


\section{ Fast Online SPIDER-ADMM Method }
In the section, we propose an online SPIDER-ADMM to
solve the online problem \eqref{eq:2}, which is equivalent to the following stochastic constrained problem:
\begin{align} \label{eq:9}
\min_{x,\{y_j\}_{j=1}^m} & \mathbb{E}_{\zeta} [ f(x;\zeta) ] + \sum_{j=1}^m g_j(y_j),  \nonumber \\
 \mbox{s.t.} & \ Ax + \sum_{j=1}^mB_jy_j =c,
\end{align}
where $f(x)=\mathbb{E}_{\zeta} [ f(x;\zeta) ]$ denotes a population risk
over an underlying data distribution. The problem \eqref{eq:9} can be viewed as
having infinite samples, so we cannot evaluate the full gradient $\nabla f(x)$.
For solving the problem \eqref{eq:9}, so we use stochastic sampling to evaluate the full gradient.
Algorithm \ref{alg:2} shows the algorithmic framework of online SPIDER-ADMM method.
In Algorithm \ref{alg:2}, we use the mini-batch samples to estimate
the full gradient, and update the variables $\{x,y_j,z\}_{j=1}^m$, which is the same as
in Algorithm \ref{alg:1}.

\begin{algorithm}[htb]
   \caption{ Online SPIDER-ADMM Algorithm }
   \label{alg:2}
\begin{algorithmic}[1]
   \STATE {\bfseries Input:} $b_1$, $b_2$, $q$, $K$, $\eta>0$ and $\rho>0$;
   \STATE {\bfseries Initialize:} $x_0 \in \mathbb{R}^d$, $y^0_j \in \mathbb{R}^p, \ j\in [m]$ and $z_0 \in \mathbb{R}^l$;
   \FOR {$k=0,1,\cdots,K-1$}
   \IF {$\mbox{mod}(k,q) = 0$}
   \STATE{} Draw $S_1$ samples with $|S_1|=b_1$, and compute $v_k = \frac{1}{b_1}\sum_{i\in S_1}\nabla f_i(x_k)$;
   \ELSE
   \STATE{} Draw $S_2$ samples with $|S_2|=b_2=\sqrt{b_1}$, and compute
            $$v_k = \frac{1}{b_2} \sum_{i\in S_2}\big(\nabla f_i(x_k) - f_i(x_{k-1})\big) + v_{k-1};$$
   \ENDIF
   \STATE{}  $ y^{k+1}_j= \arg\min_{y_j} \big\{ \mathcal {L}_{\rho}(x_k,y^{k+1}_{[j-1]},y_j,y^{k}_{[j+1:m]},z_k) + \frac{1}{2}\|y_j-y^k_j\|_{H_j}^2 \big\} $ for all $j\in [m]$;
   \STATE{}  $ x_{k+1}= \arg\min_x \hat{\mathcal {L}}_{\rho}\big( x,y^{k+1}_{[m]},z_k,v_k \big)$;
   \STATE{}  $z_{k+1} = z_k- \rho(Ax_{k+1} + \sum_{j=1}^mB_jy^{k+1}_j - c)$;
   \ENDFOR
   \STATE {\bfseries Output \ (in theory):} Chosen uniformly random from $\{x_{k},y_{[m]}^{k},z_k\}_{k=1}^{K}$.
   \STATE {\bfseries Output \ (in practice):} $\{x_{K},y_{[m]}^{K},z_K\}$.
\end{algorithmic}
\end{algorithm}

\section{Convergence Analysis}
In the section, we study the convergence properties of both the SPIDER-ADMM and online SPIDER-ADMM.
At the same time, based on our new theoretical analysis framework, we afresh analyze the convergence properties
of existing ADMM-based nonconvex optimization algorithms,
\emph{i.e.,} SVRG-ADMM and SAGA-ADMM, and derive their optimal
IFO complexity for finding an $\epsilon$-approximate stationary point.
\subsection{Convergence Analysis of SPIDER-ADMM}
In the subsection, we study convergence properties of the SPIDER-ADMM algorithm.
The detailed proofs are provided in the Appendix \ref{Appendix:A1}.
Throughout the paper, let $n_k = \lceil k/q \rceil$ such that $(n_k-1)q \leq k \leq n_kq-1$.

\begin{lemma}
 Suppose the sequence $\{x_k,y_{[m]}^k,z_k\}_{k=1}^K$ is generated from Algorithm \ref{alg:1},
 and define a \emph{Lyapunov} function $R_k$ as follows:
 \begin{align}
 R_k \!=\! & \mathcal{L}_{\rho} (x_k,y_{[m]}^k,z_k) \!+\! (\frac{9L^2}{\sigma^A_{\min}\rho}\!+\!\frac{3\sigma^2_{\max}(G)}{\sigma^A_{\min}\eta^2\rho})\|x_k\!-\!x_{k-1}\|^2  \nonumber \\
 & \!+\! \frac{2L^2}{\sigma^A_{\min}\rho b} \sum_{i=(n_k-1)q}^{k-1}\mathbb{E}\|x_{i+1}\!-\!x_i\|^2. \nonumber
\end{align}
Let $b=q$, $\eta=\frac{2\alpha\sigma_{\min}(G)}{3L} \ (0<\alpha \leq 1)$ and $\rho = \frac{\sqrt{170}\kappa_GL}{\sigma^A_{\min}\alpha}$,
then we have
\begin{align}
 \frac{1}{K}\sum_{k=0}^{K-1} \big(\|x_{k+1} - x_{k}\|^2 + \sum_{j=1}^m \|y_j^k-y_j^{k+1}\|^2\big) \leq \frac{R_{0}-R^*}{K\gamma}, \nonumber
\end{align}
where $\gamma = \min(\chi,\sigma_{\min}^H)$ with $\chi\geq \frac{\sqrt{170}\kappa_GL}{4\alpha}$ and $R^*$ is a lower bound of the function $R_k$.
\end{lemma}

Let $\theta_k = \mathbb{E} [\|x_{k+1}-x_{k}\|^2 + \|x_{k}-x_{k-1}\|^2 + \frac{1}{q}\sum_{i=(n_k-1)q}^k \|x_{i+1}-x_i\|^2 + \sum_{j=1}^m \|y_j^k-y_j^{k+1}\|^2 ]$.
Next, based on the above lemma, we give the convergence properties of SPIDER-ADMM.

\begin{theorem} \label{th:1}
 Suppose the sequence $\{x_k,y_{[m]}^k,z_k)_{k=1}^K$ is generated from Algorithm \ref{alg:1}. Let
  \begin{align}
   &\nu_1 = m\big(\rho^2\sigma^B_{\max}\sigma^A_{\max} + \rho^2(\sigma^B_{\max})^2 + \sigma^2_{\max}(H)\big), \nonumber \\
   &\nu_2 = 3(L^2+ \frac{\sigma^2_{\max}(G)}{\eta^2}),\nu_3 = \frac{18 L^2 }{\sigma^A_{\min} \rho^2} + \frac{3\sigma^2_{\max}(G) }{\sigma^A_{\min}\eta^2\rho^2}, \nonumber
 \end{align}
 and let $b=q$, $\eta = \frac{2\alpha\sigma_{\min}(G)}{3L} \ (0<\alpha \leq 1)$, and
 $\rho = \frac{\sqrt{170}\kappa_GL}{\sigma^A_{\min}\alpha}$,
 then we have
 \begin{align}
  \frac{1}{K} \sum_{k=1}^K \mathbb{E}\big[ \mbox{dist}(0,\partial L(x_k,y_{[m]}^k,z_k))^2\big] &\leq\frac{\nu_{\max}}{K}\sum_{k=1}^{K-1} \theta_k \nonumber\\
  &\leq  \frac{3\nu_{\max}(R_{0}-R^*)}{K\gamma}, \nonumber
 \end{align}
 where $\gamma = \min(\chi,\sigma_{\min}^H)$ with $\chi\geq \frac{\sqrt{170}\kappa_GL}{4\alpha}$, $\nu_{\max}=\max\{\nu_1,\nu_2,\nu_3\}$ and
 $R^*$ is a lower bound of the function $R_k$.
 It implies that the iteration number $K$ satisfies
 \begin{align}
  K = \frac{3\nu_{\max}(R_{0}-R^*)}{\epsilon \gamma}, \nonumber
 \end{align}
 then $(x_{k^*},y_{[m]}^{k^*},z_{k^*})$ is an $\epsilon$-approximate stationary point of \eqref{eq:2}, where $k^* = \mathop{\arg\min}_{k}\theta_{k}$.
\end{theorem}

\begin{remark}
 Theorem \ref{th:1} shows that the SPIDER-ADMM has $O(1/K)$ convergence rate. Moreover, given $b=q=\sqrt{n}$, $\eta = \frac{2\alpha\sigma_{\min}(G)}{3L} \ (0<\alpha \leq 1)$ and
 $\rho = \frac{\sqrt{170}\kappa_GL}{\sigma^A_{\min}\alpha}$, the SPIDER-ADMM has the optimal
 IFO of $\mathcal{O}(n+n^{\frac{1}{2}}\epsilon^{-1})$ for finding an $\epsilon$-approximate stationary point. In particular,
 we can choose $\alpha \in (0,1]$ according to different problems to obtain appropriate step-size $\eta$ and penalty parameter $\rho$, \emph{e.g.}, set $\alpha=1$,
 we have $\eta=\frac{2\sigma_{\min}(G)}{3L}$ and $\rho=\frac{\sqrt{170}\kappa_GL}{\sigma^A_{\min}}$.
\end{remark}
\subsection{Convergence Analysis of Online SPIDER-ADMM}
In the subsection, we study convergence properties of the online SPIDER-ADMM algorithm.
The detailed proofs are provided in the Appendix \ref{Appendix:A2}.

\begin{lemma}
 Suppose the sequence $\{x_k,y_{[m]}^k,z_k\}_{k=1}^K$ is generated from Algorithm \ref{alg:2},
 and define a \emph{Lyapunov} function $\Phi_k$ as follows:
 \begin{align}
 \Phi_k = & \mathcal{L}_{\rho} (x_k,y_{[m]}^k,z_k) + (\frac{9L^2}{\sigma^A_{\min}\rho}+\frac{3\sigma^2_{\max}(G)}{\sigma^A_{\min}\eta^2\rho})\|x_k-x_{k-1}\|^2 \nonumber \\
 & + \frac{2L^2}{\sigma^A_{\min}\rho b_2} \sum_{i=(n_k-1)q}^{k-1}\mathbb{E}\|x_{i+1}-x_i\|^2. \nonumber
\end{align}
Let $b_2=q$, $\eta=\frac{2\alpha\sigma_{\min}(G)}{3L} \ (0<\alpha \leq 1)$ and $\rho= \frac{\sqrt{170}\kappa_GL}{\sigma^A_{\min}\alpha}$,
then we have
\begin{align}
& \frac{1}{K}\sum_{k=0}^{K-1} \big(\|x_{k+1} - x_{k}\|^2 + \sum_{j=1}^m \|y_j^k-y_j^{k+1}\|^2\big) \nonumber \\
 &\leq \frac{\Phi_{0}-\Phi^*}{K\gamma} + \frac{2\delta^2}{b_1L\gamma} + \frac{72\delta^2}{\sigma^A_{\min}b_1\rho\gamma}, \nonumber
\end{align}
where $\gamma = \min(\chi,\sigma_{\min}^H)$ with $\chi\geq \frac{\sqrt{170}\kappa_GL}{4\alpha}$ and $\Phi^*$ is a lower bound of the function $\Phi_k$.
\end{lemma}

Let $\theta_k = \mathbb{E} [ \|x_{k+1}-x_{k}\|^2+\|x_{k}-x_{k-1}\|^2+\frac{1}{q}\sum_{i=(n_k-1)q}^k \|x_{i+1}-x_i\|^2+ \sum_{j=1}^m \|y_j^k-y_j^{k+1}\|^2 ]$.
\begin{theorem} \label{th:2}
 Suppose the sequence $\{x_k,y_{[m]}^k,z_k)_{k=1}^K$ is generated from Algorithm \ref{alg:2}. Let
  \begin{align}
   &\nu_1 = m\big(\rho^2\sigma^B_{\max}\sigma^A_{\max} + \rho^2(\sigma^B_{\max})^2 + \sigma^2_{\max}(H)\big), \nonumber \\
   &\nu_2 = 3(L^2+ \frac{\sigma^2_{\max}(G)}{\eta^2}),\nu_3 = \frac{18 L^2 }{\sigma^A_{\min} \rho^2} + \frac{3\sigma^2_{\max}(G) }{\sigma^A_{\min}\eta^2\rho^2}, \nonumber
 \end{align}
 and let $b_2=q=\sqrt{b_1}$, $\eta = \frac{2\alpha\sigma_{\min}(G)}{3L} \ (0<\alpha \leq 1)$ and
 $\rho = \frac{\sqrt{170}\kappa_GL}{\sigma^A_{\min}\alpha}$,
 then we have
  \begin{align}
  &\frac{1}{K} \sum_{k=1}^K \mathbb{E} \big[ \mbox{dist}(0,\partial L(x_k,y_{[m]}^k,z_k))^2\big]\leq \frac{\nu_{\max}}{K}\sum_{k=1}^{K-1} \theta_k
   + \frac{w}{b_1} \nonumber \\
  & \leq  \frac{3\nu_{\max}(\Phi_{0}-\Phi^*)}{K\gamma} + \frac{6\nu_{\max}\delta^2}{b_1\gamma}(\frac{1}{L} + \frac{36}{\sigma^A_{\min}\rho})+ \frac{w}{b_1}, \nonumber
\end{align}
 where $w= 12\delta^2\max\{1, \frac{6}{\sigma^A_{\min}\rho^2}\}$, $\gamma = \min(\chi,\sigma_{\min}^H)$ with $\chi\geq \frac{\sqrt{170}\kappa_GL}{4\alpha}$,
 $\nu_{\max}=\max\{\nu_1,\nu_2,\nu_3\}$ and
 $\Phi^*$ is a lower bound of the function $\Phi_k$.
 It implies that $K$ and $b_1$ satisfy
 \begin{align}
  K \!=\! \frac{6\nu_{\max}(\Phi_{0}-\Phi^*)}{\epsilon \gamma}, \ b_1 \!=\! \frac{12\nu_{\max}\delta^2}{\epsilon\gamma}(\frac{1}{L}+\frac{36}{\sigma^A_{\min}\rho}) \!+\! \frac{2w}{\epsilon}, \nonumber
 \end{align}
 then $(x_{k^*},y_{[m]}^{k^*},z_{k^*})$ is an $\epsilon$-approximate stationary point of \eqref{eq:2}, where $k^* = \mathop{\arg\min}_{k}\theta_{k}$.
\end{theorem}
\begin{remark}
 Theorem \ref{th:2} shows that given $b_2=q=\sqrt{b_1}$, $\eta = \frac{2\alpha\sigma_{\min}(G)}{3L} \ (0<\alpha \leq 1)$,
 $\rho = \frac{\sqrt{170}\kappa_GL}{\sigma^A_{\min}\alpha}$ and $b_1=\mathcal{O}(\epsilon^{-1})$, the online SPIDER-ADMM has the optimal
 IFO of $\mathcal{O}(\epsilon^{-\frac{3}{2}})$ for finding an $\epsilon$-approximate stationary point.
\end{remark}

\subsection{Convergence Analysis of Non-convex SVRG-ADMM}
In the subsection, we extend the existing nonconvex SVRG-ADMM method \cite{huang2016stochastic,zheng2016stochastic}
to the multiple variables setting for solving the problem \eqref{eq:2}.
The SVRG-ADMM algorithm is described in Algorithm 3 given in the Appendix \ref{Appendix:A3}.
Next, we analyze convergence properties of the SVRG-ADMM algorithm,
and derive its optimal IFO complexity.

\begin{lemma}
Suppose the sequence $\{(x^{s}_t,y_{[m]}^{s,t},z^{s}_t)_{t=1}^M\}_{s=1}^S$ is generated from Algorithm \ref{alg:3},
 and define a \emph{Lyapunov} function:
 \begin{align}
 \Gamma^s_t\!=& \mathbb{E}\big[\mathcal{L}_{\rho} (x^s_t,y_{[m]}^{s,t},z^s_t) \!+\! (\frac{3\sigma^2_{\max}(G)}{\sigma^A_{\min}\eta^2\rho} \!+\! \frac{9L^2}{\sigma^A_{\min}\rho})\|x^s_{t}-x^s_{t-1}\|^2 \nonumber \\
 &+ \frac{9L^2 }{\sigma^A_{\min}\rho b}\|x^s_{t-1}-\tilde{x}^s\|^2 + c_t\|x^s_{t}-\tilde{x}^s\|^2\big], \nonumber
 \end{align}
 where the positive sequence $\{c_t\}$ satisfies, for $s =1,2,\cdots,S$
 \begin{equation*}
  c_t= \left\{
  \begin{aligned}
  & \frac{18 L^2 }{\sigma^A_{\min}\rho b} +
     \frac{L}{b} + (1+\beta)c_{t+1}, \ 1 \leq t \leq M, \\
  & 0, \ t \geq M+1.
  \end{aligned}
  \right.\end{equation*}
Let $M=n^{\frac{1}{3}}$, $b=n^{\frac{2}{3}}$, $\eta = \frac{\alpha\sigma_{\min}(G)}{5L} \ (0< \alpha \leq 1)$ and
$\rho = \frac{2\sqrt{231}\kappa_G L}{\sigma^A_{\min}\alpha}$, we have
\begin{align}
\frac{1}{T}\sum_{s=1}^S \sum_{t=0}^{M-1} & \big(\sigma_{\min}^H\sum_{j=1}^m \|y_j^{s,t}-y_j^{s,t+1}\|^2 + \chi_t \|x^s_{t+1}-x^s_t\|^2 \nonumber \\
& + \frac{L}{2b}\|x^s_t-\tilde{x}^s\|^2 \big) \leq \frac{\Gamma^1_0 - \Gamma^*}{T}  .
\end{align}
where $T=MS$, $\chi_t \geq \frac{\sqrt{231}\kappa_G L}{2\alpha} > 0$ and $\Gamma^*$ denotes a lower bound of function $\Gamma^s_t$.
\end{lemma}

Let $\theta^s_{t} = \mathbb{E} [ \|x^s_{t+1}-x^s_{t}\|^2 + \|x^s_{t}-x^s_{t-1}\|^2 + \frac{1}{b}(\|x^s_{t}-\tilde{x}^s\|^2 + \|x^s_{t-1}-\tilde{x}^s\|^2 ) + \sum_{j=1}^m\|y_j^{s,t} - y_j^{s,t+1}\|^2 ]$.
\begin{theorem} \label{th:3}
 Suppose the sequence $\{(x^{s}_t,y_{[m]}^{s,t},z^{s}_t)_{t=1}^M\}_{s=1}^S$ is generated from Algorithm \ref{alg:3} . Let
 \begin{align}
 & \nu_1 = m\big(\rho^2\sigma^B_{\max}\sigma^A_{\max} + \rho^2(\sigma^B_{\max})^2 + \sigma^2_{\max}(H)\big), \nonumber \\
 & \nu_2 = 3L^2 + \frac{3\sigma^2_{\max}(G)}{\eta^2},
 \ \nu_3 = \frac{9L^2 }{\sigma^A_{\min}\rho^2} + \frac{3\sigma^2_{\max}(G)}{\sigma^A_{\min}\eta^2\rho^2}, \nonumber
\end{align}
 and given $M=n^{\frac{1}{3}}$, $b=n^{\frac{2}{3}}$, $\eta = \frac{\alpha\sigma_{\min}(G)}{5L} \ (0 < \alpha \leq 1)$
 and $\rho = \frac{2\sqrt{231}\kappa_G L}{\sigma^A_{\min}\alpha}$,
 then we have
 \begin{align}
\frac{1}{T}\sum_{s=1}^S \sum_{t=0}^{M-1}\mathbb{E}\big[ \mbox{dist}(0,\partial L(x^s_t,y_{[m]}^{s,t},z^s_t))^2\big] \leq \frac{2\nu_{\max}(\Gamma^1_0 - \Gamma^*)}{\gamma T}, \nonumber
\end{align}
where $\gamma = \min(\sigma_{\min}^H,\frac{L}{2},\chi_t)$, $\nu_{\max}= \max(\nu_1,\nu_2,\nu_3)$ and
$\Gamma^*$ is a lower bound of function $\Gamma^s_t$.
It implies that the whole iteration number $T=MS$ satisfies
 \begin{align}
 T = \frac{2\nu_{\max}(\Gamma^1_0 - \Gamma^*)}{\epsilon \gamma}, \nonumber
 \end{align}
then $(x^{s^*}_{t^*},y_{[m]}^{s^*,t^*},z^{s^*}_{t^*})$ is an $\epsilon$-stationary point of \eqref{eq:2}, where $(t^*,s^*) = \mathop{\arg\min}_{t,s}\theta^s_{t}$.
\end{theorem}
\begin{remark}
 Theorem \ref{th:3} shows that given $M=n^{\frac{1}{3}}$, $b=n^{\frac{2}{3}}$, $\eta = \frac{\alpha\sigma_{\min}(G)}{5L} \ (0<\alpha \leq 1)$ and
 $\rho = \frac{2\sqrt{231}\kappa_GL}{\sigma^A_{\min}\alpha}$, the non-convex SVRG-ADMM has the optimal
 IFO complexity of $\mathcal{O}(n+n^{\frac{2}{3}}\epsilon^{-1})$ for finding an $\epsilon$-approximate stationary point.
\end{remark}
\subsection{Convergence Analysis of Non-convex SAGA-ADMM}
In the subsection, we extend the existing nonconvex SAGA-ADMM method \cite{huang2016stochastic}
to the multiple variables setting for solving the problem \eqref{eq:2}.
The SAGA-ADMM algorithm is described in Algorithm 4 given in the Appendix \ref{Appendix:A4}.
Next, we analyze convergence properties of non-convex SAGA-ADMM,
and derive its the optimal IFO complexity.

\begin{lemma}
 Suppose the sequence $\{x_t,y_{[m]}^t,z_t\}_{t=1}^T$ is generated from Algorithm \ref{alg:4},
 and define a \emph{Lyapunov} function
 \begin{align}
 & \Omega_t \!= \!\mathbb{E}\big[ \mathcal{L}_{\rho} (x_t,y_{[m]}^t,z_t) \!+ \! (\frac{3\sigma^2_{\max}(G)}{\sigma^A_{\min}\eta^2\rho} \!+ \! \frac{9L^2}{\sigma^A_{\min}\rho}) \|x_{t}\!-\!x_{t-1}\|^2 \nonumber \\
 & + \frac{9L^2}{\sigma^A_{\min}\rho b}\frac{1}{n} \sum_{i=1}^n\|x_{t-1}-u^{t-1}_i\|^2 + c_t\frac{1}{n} \sum_{i=1}^n\|x_{t}-u^t_i\|^2 \big], \nonumber
\end{align}
 where the positive sequence $\{c_t\}$ satisfies
 \begin{equation*}
  c_t \!=\! \left\{
  \begin{aligned}
  & \frac{18L^2 }{\sigma^A_{\min}\rho b} \!+\! \frac{L}{b} \!+\! (1-p)(1+\beta)c_{t+1}, \ 0 \leq t \leq T-1, \\
  & 0, \ t \geq T,
  \end{aligned}
  \right.\end{equation*}
where $p$ denotes probability of an index $i$ being in $\mathcal{I}_t$.
Further, let $b= n^{\frac{2}{3}}$, $\eta = \frac{\alpha\sigma_{\min}(G)}{17L} \ (0 < \alpha \leq 1)$ and $\rho = \frac{2\sqrt{2031}\kappa_G}{\sigma^A_{\min}\alpha}$
we have
\begin{align}
 \frac{1}{T} \sum_{t=1}^T & \big( \sigma_{\min}^H\sum_{j=1}^m \|y_j^t-y_j^{t+1}\|^2 + \chi_t \|x_t-x_{t+1}\|^2 \nonumber \\
 &+  \frac{L}{2b}\frac{1}{n}\sum_{i=1}^n\|x_t-u^t_i\|^2 \big) \leq \frac{\Omega_0 - \Omega^*}{T}, \nonumber
\end{align}
where $\chi_t \geq \frac{\sqrt{2031}\kappa_GL}{2\alpha} >0$ and $\Omega^*$ denotes a lower bound of function $\Omega_t$.
\end{lemma}

Let $\theta_t =  \mathbb{E} [ \|x_{t+1}-x_t\|^2 + \|x_t-x_{t-1}\|^2 + \frac{1}{b n}\sum^n_{i=1} (\|x_t-u^t_i\|^2 + \|x_{t-1}-u^{t-1}_i\|^2)+\sum_{j=1}^m \|y_j^t-y_j^{t+1}\|^2 ]$.

\begin{theorem} \label{th:4}
Suppose the sequence $\{x_t,y_{[m]}^t,z_t\}_{t=1}^T$ is generated from Algorithm \ref{alg:4}. Let
\begin{align}
  &\nu_1 =m\big(\rho^2\sigma^B_{\max}\sigma^A_{\max} + \rho^2(\sigma^B_{\max})^2 + \sigma^2_{\max}(H)\big), \nonumber \\
  & \nu_2 = 3L^2 + \frac{3\sigma^2_{\max}(G)}{\eta^2}, \ \nu_3 = \frac{9L^2 }{\sigma^A_{\min}\rho^2} + \frac{3\sigma^2_{\max}(G)}{\sigma^A_{\min}\eta^2\rho^2}, \nonumber
 \end{align}
 and given $b= n^{\frac{2}{3}}$, $\eta = \frac{\alpha\sigma_{\min}(G)}{17L} \ (0 < \alpha \leq 1)$ and $\rho = \frac{2\sqrt{2031}\kappa_G}{\sigma^A_{\min}\alpha}$,
 then we have
 \begin{align}
 \frac{1}{T}\sum_{t=1}^T \mathbb{E}\big[ \mbox{dist}(0,\partial L(x_t,y_{[m]}^t,z_t))^2\big] \leq \frac{2\nu_{\max}(\Omega_0 - \Omega^*)}{\gamma T},  \nonumber
 \end{align}
 where $\gamma = \min(\sigma_{\min}^H, \frac{L}{2}, \chi_{t})$ with $\chi_t \geq \frac{\sqrt{2031}\kappa_GL}{2\alpha} >0$, $\nu_{\max}= \max(\nu_1,\nu_2,\nu_3)$ and
 $\Omega^*$ is a lower bound of function $\Omega_t$. It implies that the iteration number $T$ satisfies
 \begin{align}
  T = \frac{2\nu_{\max}}{\epsilon \gamma}(\Omega_0- \Omega^*),  \nonumber
 \end{align}
 then $(x_{t^*},y_{[m]}^{t^*},z_{t^*})$ is an $\epsilon$-approximate stationary point of \eqref{eq:2}, where $t^* = \mathop{\arg\min}_{ 1\leq t\leq T}\theta_{t}$.
\end{theorem}
\begin{remark}
 Theorem \ref{th:4} shows that given $b=n^{\frac{2}{3}}$, $\eta = \frac{\alpha\sigma_{\min}(G)}{17L} \ (0<\alpha \leq 1)$ and
 $\rho = \frac{2\sqrt{2031}\kappa_GL}{\sigma^A_{\min}\alpha}$, the non-convex SAGA-ADMM has the optimal
 IFO of $\mathcal{O}(n+n^{\frac{2}{3}}\epsilon^{-1})$ for finding an $\epsilon$-approximate stationary point.
\end{remark}

\begin{remark}
\textbf{Our contributions} on convergence analysis of both the non-convex SVRG-ADMM and SAGA-ADMM are given as follows:
\begin{itemize}
\item We extend both the existing non-convex SVRG-ADMM and SAGA-ADMM to the multi-block setting
      for solving the problem \eqref{eq:2};
\item We not only give its optimal IFO complexity of $\mathcal{O}(n+n^{\frac{2}{3}}\epsilon^{-1})$,
      but also provide the \emph{specific and simple} choice on the step-size $\eta$ and penalty parameter $\rho$.
\end{itemize}
\end{remark}
\begin{table}
  \centering
  \caption{Real datasets } \label{tab:2}
  \begin{tabular}{c|c|c|c}
  \hline
  datasets & $\#samples$ & $\#features$ & $\#classes$ \\ \hline
  \emph{a9a}   & 32,561  & 123 & 2\\
  \emph{w8a} & 64,700  & 300 & 2\\
  \emph{ijcnn1} & 126,702 & 22 & 2 \\
  \emph{covtype.binary} & 581,012 & 54 & 2\\ \hline
  \emph{letter} & 15,000 & 16 & 26 \\
  \emph{sensorless} & 58,509 &48 & 11 \\
  \emph{mnist} & 60,000 & 780 & 10 \\
  \emph{covtype} & 581,012  & 54 & 7 \\
  \hline
  \end{tabular}
\end{table}

\section{ Experiments }
In this section, we will compare the proposed algorithm (SPIDER-ADMM) with the existing non-convex algorithms \big( nc-ADMM \cite{jiang2019structured},
nc-SVRG-ADMM \cite{huang2016stochastic,zheng2016stochastic}, nc-SAGA-ADMM \cite{huang2016stochastic} and nc-SADMM \cite{huang2018mini} \big) on two applications:
1) Graph-guided binary classification; 2) Multi-task learning.
In the experiment, we use some publicly available datasets\footnote{
These data are from the LIBSVM website (www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/).},
which are summarized in Table \ref{tab:2}. All algorithms are implemented in MATLAB,
and all experiments are performed on a PC with an Intel i7-4790 CPU and 16GB memory.
\begin{figure}[htbp]
\centering
\subfigure[\emph{a9a}]{\includegraphics[width=0.23\textwidth]{figs/a9a_obj}}
\subfigure[\emph{w8a}]{\includegraphics[width=0.23\textwidth]{figs/w8a_obj}}
\subfigure[\emph{ijcnn1}]{\includegraphics[width=0.23\textwidth]{figs/ijcnn1_obj}}
\subfigure[\emph{covtype.binary}]{\includegraphics[width=0.23\textwidth]{figs/covtypeb_obj}}
\caption{Objective value \emph{versus} CPU time of the \emph{nonconvex} graph-guided binary classification model on some real datasets.} \label{fig:1}
\vspace{-1em}
\end{figure}
\subsection{ Graph-Guided Binary Classification }
In the subsection, we focus on the binary classification task.
Specifically, given a set of training samples $(a_i,b_i)_{i=1}^n$,
where $a_i\in \mathbb{R}^d$, $b_i \in \{-1,1\}$,
then we solve the following nonconvex empirical loss minimization problem:
\begin{align} \label{eq:80}
 \min_{x\in \mathbb{R}^d} \ \frac{1}{n}\sum_{i=1}^n f_i(x) + \lambda\|Ax\|_1,
\end{align}
where $f_i(x)=\frac{1}{1+\exp(b_i a_i^Tx)}$ is the nonconvex sigmoid loss function.
We use the nonsmooth regularizer \emph{i.e.,}
graph-guided fused lasso \cite{kim2009multivariate}, and $A$ decodes the sparsity pattern of graph,
which is obtained by sparse precision matrix estimation \cite{friedman2008sparse}.
To solve the problem \eqref{eq:80},
we give an auxiliary variable $y$ with the constraint $y = Ax$.
In the experiment, we fix the parameter $\lambda=10^{-5}$, and use the same initial solution $x_0$
from the standard normal distribution for all algorithms.

Figure \ref{fig:1} shows that the objective values of our SPIDER-ADMM method faster decrease than those of other methods,
as CPU time consumed increases.
Thus, these results demonstrate that our method has a relatively faster convergence rate than other methods.

\begin{figure}[htbp]
\centering
\subfigure[\emph{letter}]{\includegraphics[width=0.23\textwidth]{figs/letter_obj}}
\subfigure[\emph{sensorless}]{\includegraphics[width=0.23\textwidth]{figs/sensorless_obj}}
\subfigure[\emph{mnist}]{\includegraphics[width=0.23\textwidth]{figs/mnist_obj}}
\subfigure[\emph{covtype}]{\includegraphics[width=0.23\textwidth]{figs/covtype_obj}}
\caption{Objective value \emph{versus} CPU time of the \emph{nonconvex} multi-task learning on some real datasets.} \label{fig:2}
\vspace{-1em}
\end{figure}
\subsection{ Multi-Task Learning }
In this subsection, we focus on the multi-task learning task with sparse and low-rank structures.
Specifically, given a set of training samples $(a_i,b_i)_{i=1}^n$,
where $a_i\in \mathbb{R}^d$ and $b_i \in \{1,2,\cdots,c\}$, then let $D\in \mathbb{R}^{n\times c}$ with $D_{ij}=1$ if $j=b_i$,
and $D_{ij}=0$ otherwise.
This multi-task learning is equivalent to solving the following nonconvex problem:
\begin{equation}
 \min_{X \in \mathbb{R}^{c\times d}} \ \frac{1}{n}\sum_{i=1}^n f_i(X) + \lambda_1\sum_{ij}\kappa(|X_{ij}|) + \lambda_2\|X\|_* ,
\end{equation}
where $f_i(X) = \log(\sum_{j=1}^c\exp(X_{j,.}a_i)) - \sum_{j=1}^c D_{ij}X_{j,.}a_i$ is a multinomial logistic loss
function, $\kappa(|X_{ij}|)=\beta\log(1+\frac{|X_{ij}|}{\alpha})$ is the
nonconvex log-sum penalty function \cite{candes2008enhancing}.
Next, we change the above problem into the following form:
\begin{align} \label{eq:82}
 \min & \quad \frac{1}{n}\sum_{i=1}^n \bar{f}_i(X) + \lambda_1\kappa_0\|Y_1\|_1 + \lambda_2\|Y_2\|_* \\
 \mbox{s.t.} & \quad AX + B_1Y_1 + B_2Y_2 = 0, \nonumber
\end{align}
where $\bar{f}_i(X) = f_i(X) + \lambda_1\big( \sum_{ij}\kappa(|X_{ij}|)-\kappa_0\|X\|_1 \big)$, and $\kappa_0=\kappa'(0)$.
Here $A = [I_c;I_c] \in \mathbb{R}^{2c\times c}$, $B_1=[-I_c;0] \in \mathbb{R}^{2c\times c}$ and $B_2=[0;-I]$.
By the Proposition 2.3 in \citet{yao2016efficient},  $\bar{f}_i(X)$ is nonconvex and smooth.
In the experiment, we fix the parameters $\lambda_1 = 10^{-5}$ and $\lambda_2 = 10^{-4}$,
and use the same initial solution $x_0$ from the standard normal distribution for all algorithms.

Figure \ref{fig:2} shows that objective values of our SPIDER-ADMM faster decrease than those of the other methods,
as CPU time consumed increases.
Similarly, these results also demonstrate that
our method has a relatively faster convergence rate than other methods.


\section{Conclusion}
In the paper, we proposed a faster stochastic ADMM method (\emph{i.e.,} SPIDER-ADMM)
for nonconvex optimization. Moreover, we proved that the SPIDER-ADMM achieves a lower IFO complexity of $\mathcal{O}(n+n^{1/2}\epsilon^{-1})$.
Further, we extended the SPIDER-ADMM to the online setting, and
proposed a faster online ADMM method (\emph{i.e.,} online SPIDER-ADMM).
As one of major contribution of this paper, we provided a new theoretical analysis framework for the nonconvex stochastic ADMM methods
with providing an optimal IFO complexity.
Based on our new theoretical analysis framework, we studied the unsolved optimal IFO complexity of
the existing non-convex SVRG-ADMM and SAGA-ADMM methods, and also proved that they reach an IFO complexity of $\mathcal{O}(n+n^{2/3}\epsilon^{-1})$.
In the future work, we can apply the stage-wise stochastic momentum technique \cite{chen2018universal} to accelerate our algorithms.

\section*{Acknowledgments}
We thank the anonymous reviewers for their helpful comments.
F.H. and H.H. were partially supported by U.S. NSF IIS 1836945,  IIS 1836938,  DBI 1836866,  IIS 1845666,  IIS 1852606,  IIS 1838627, IIS 1837956.
S.C. was partially supported by the NSFC under Grant No. 61806093 and No. 61682281,
and the Key Program of NSFC under Grant No. 61732006.



\nocite{langley00}

