---
abstract: 'In recent years, stochastic gradient descent (SGD) and its variants have been the dominant optimization methods for training deep neural networks. However, SGD suffers from limitations such as the lack of theoretical guarantees, vanishing gradients, excessive sensitivity to input, and difficulties solving highly non-smooth constraints and functions. To overcome these drawbacks, alternating minimization-based methods for deep neural network optimization have attracted fast-increasing attention recently. As an emerging and open domain, however, several new challenges need to be addressed, including: 1) there is no guarantee of global convergence under mild, practical conditions, and 2) cubic time complexity in the size of feature dimensions. We therefore propose a novel Deep Learning Alternating Minimization (DLAM) algorithm to deal with these two challenges. Our innovative inequality-constrained formulation infinitely approximates the original problem with non-convex equality constraints, enabling our proof of global convergence of the DLAM algorithm under mild, practical conditions. The time complexity is successfully reduced from $O(d^3)$ to $O(d^2)$ via a dedicated algorithm design for subproblems that is enhanced by iterative quadratic approximations and backtracking. Experiments on benchmark datasets demonstrate the effectiveness of our proposed DLAM algorithm.'
author:
- Junxiang Wang
- Fuxun Yu
- Xiang Chen
- Liang Zhao
bibliography:
- 'example\_paper.bib'
title: Global Convergence of the Alternating Minimization Algorithm for Deep Neural Network Problems
---

Introduction
============

Stochastic gradient descent (SGD) and its variants such as adaptive momentum estimation (Adam) [@kingma2014adam] have become popular optimization methods for training deep neural networks. These methods split a dataset into multiple batches and then optimize them sequentially by gradient descent in each epoch. SGD has two main advantages: not only is it simple to implement, it can also be applied in online settings where new coming training data are used to train models. However, while many researchers have provided solid theoretical guarantees on the convergence of SGD [@kingma2014adam; @j.2018on; @sutskever2013importance], the assumptions of their proofs cannot be applied to problems involving deep neural networks, which are highly nonsmooth and nonconvex. Aside from the lack of theoretical guarantees, several additional drawbacks restrict the applications of SGD. It suffers from the gradient vanishing problem, meaning that the error signal diminishes as the gradient is backpropagated, which prevents the neural networks from utlizing further training [@taylor2016training], and the gradient of the activation function is highly sensitive to the input (i.e. poor conditioning), so a small change in the input can lead to a dramatic change in the gradient.\
To tackle these intrinsic drawbacks of gradient descent optimization methods, alternating minimization methods have started to attract attention as a potential way to solve deep learning problems. Here, the loss function of a deep neural network is reformulated as a nested function associated with multiple linear and nonlinear transformations across multi-layers. This nested structure is then decomposed into a series of linear and nonlinear equality constraints by introducing auxiliary variables. The linear and nonlinear equality constraints generate multiple subproblems, which can be minimized alternately. Some recent alternating minimization methods have focused on applying the Alternating Direction Method of Multipliers (ADMM) [@taylor2016training] and Block Coordinate Descent (BCD) [@zeng2018global], with empirical evaluations demonstrating good scalability in terms of the number of layers and high accuracy on the test sets, especially for neural networks that are very deep, thanks to parallelism [@taylor2016training]. These methods also avoid gradient vanishing problems and allow for non-differentiable activation functions such as binarized neural networks [@courbariaux2015binaryconnect], as well as allowing for complex non-smooth regularization and the constraints that are increasingly important for deep neural architectures that are required to satisfy practical requirements such as interpretability, energy-efficiency, and cost awareness [@carreira2014distributed].\
However, as an emerging domain, alternating minimization for deep model optimization suffers from a number of unsolved challenges including: **1. The lack of global convergence guarantee with mild and practical conditions.** Most existing alternating minimization methods for deep learning cannot provide convergence guarantees [@taylor2016training]; the few that do impose guarantee request for conditions that are unrealistic or hard to satisfy for most practical applications [@lau2018proximal]. This is because in order to enable the alternating minimization methods to disentangle the nested functions into subproblems, additional equality constraints must be added to coordinate these subproblems. However, these equality constraints are inherently nonconvex due to the nonlinearity of the activation functions (e.g., sigmoids), making it prohibitively difficult to obtain a global minimum and convergence [@taylor2016training; @lau2018proximal]. **2. Cubic time complexity in the size of the feature dimensions.** Existing methods are generally very time-consuming. This is because in order to quantify the auxiliary variables, existing alternating minimization methods typically require matrix inversion computation, which has $O(d^3)$ time complexity, where $d$ is the dimension of features, and is thus prohibitively costly [@taylor2016training; @boyd2011distributed].\
In order to simultaneously address these technical problems, we propose a new formulation of the deep neural network problem, along with a novel Deep Learning Alternating Minimization (DLAM) algorithm. The proposed framework is highly generic and sufficiently flexible to be utilized in common fully-connected deep neural network models, as well as being easily extendable to other models such as convolutional neural networks [@krizhevsky2012imagenet] and recurrent neural networks [@mikolov2010recurrent]. Specifically, we, for the first time, transform the original deep neural network optimization problem into an inequality-constrained problem that can be infinitely approximate to the original one. Applying this innovation to a inequality-constraint based transformation ensures the convexity of all subproblems, and hence easily ensures global minima. The operation of matrix inversion is avoided by the quadratic approximation technique and a backtracking algorithm, which also speeds up the convergence of the DLAM algorithm. Moreover, while existing methods require typically strict and complex conditions, such as Kurdyka-Łojasiewicz (KL) properties [@lau2018proximal] to prove convergence, our proposed methods requires simple and mild conditions to guarantee convergence and covers most of the commonly-used loss functions and activation functions. Last but not least, the choice of hyperparameters has no effect on the convergence of our DLAM algorithm theoretically. Our contributions in this paper include:

-   We propose a novel formulation for deep neural network optimization. The deeply nested activation functions are disentangled into separate functions innovatively coordinated by inequality constraints that are inherently convex.

-   We present a novel and efficient DLAM algorithm. A quadratic approximation technique and a backtracking algorithm are utilized to avoid matrix inversion, thus reducing the computational cost to $O(d^2)$ considerably. Every subproblem has a closed-form solution, further boosting the efficiency.

-   We investigate several attractive convergence properties of the DLAM algorithm under mild conditions. The model assumptions are very mild, ensuring that most deep learning problems will satisfy our assumptions. The new DLAM algorithm is guaranteed to converge to a critical point.

-   We conduct experiments on benchmark datasets to validate our proposed DLAM algorithm. Experiments on two benchmark datasets show that the new algorithm performs well compared with SGD or its variants and ADMM.

The rest of paper is organized as follows. In Section \[sec:related work\], we summarize recent research related to this topic. In Section \[sec:algorithm\], we present the problem formulation and the new DLAM algorithm. In Section \[sec:convergence\], we introduce the main convergence results for the DLAM algorithm. Section \[sec:experiment\] reports the results of the extensive experiments conducted to validate the convergence and effectiveness of the new DLAM. Section \[sec:conclusion\] concludes by summarizing the research.

Related Work {#sec:related work}
============

All of the existing works on optimization methods in deep neural network problems falls into two major classes: stochastic gradient descent methods and alternating minimization methods. This research related to both is discussed in this section.\
**Stochastic gradient descent methods:** The renaissance of SGD can be traced back to 1951, when Robbins and Monro published the first paper[@robbins1951textordfemininea]. The famous back-propagation algorithm was introduced by Rumelhart et al. [@rumelhart1986learning]. Many variants of SGD methods have since been presented, including the use of Polyak momentum, which accelerates the convergence of iterative methods [@polyak1964some], and research by Sutskever et al., who highlighted the importance of Nesterov momentum and initialization [@sutskever2013importance]. Many well-known SGD methods that incorporate with adaptive learning rates have been proposed by the deep learning community, including AdaGrad [@duchi2011adaptive], RMSProp [@tielemandivide], Adam [@kingma2014adam] and AMSGrad [@j.2018on].\
**Alternating minimization methods for deep learning:** Previous work on the application of alternating minimization algorithms to deep learning problems can be categorized into two main types. The first research strand proposes the use of alternating minimization algorithms for specific applications. For example, Taylor et al. presented an Alternating Direction Method of Multipliers (ADMM) algorithm to transform a fully-connected neural network problem into an equality-constrained problem, where many subproblems split by ADMM can be solved in parallel [@taylor2016training], while Zhang et al. handled very deep supervised hashing (VDSH) problems by utlizing an ADMM algorithm to overcome issues related to vanishing gradients and poor computational efficiency [@zhang2016efficient]. Zhang and Bastiaan trained a deep neural network by utlizing ADMM with a graph [@zhang2017training] and Askari et al. introduced a new framework for multilayer feedforward neural networks and solved the new framework using block coordinate descent (BCD) methods [@askari2018lifted]. Others have proposed novel alternating minimization methods and proved their convergence results. For instance, Carreira and Wang suggested a method involving the use of auxiliary coordinates (MAC) to replace a nested neural network with a constrained problem without nesting [@carreira2014distributed]. Lin and Yao, and Lau et al. both proposed BCD algorithms, proving its convergence via the Kurdyka-Łojasiewicz (KL) property [@lau2018proximal], while Choromanska et al. proposed a BCD algorithm for training deep feedforward neural networks based on the concept of co-activation memory [@choromanska2018beyond] and a BCD algorithm with R-linear convergence was proposed by Zhang and Brand to train Tikhonov regularized deep neural networks [@zhang2017convergent]. However, most of these researchers focused on specific applications of neural networks rather than their general formulations. Even though several did discuss general neural network problems and provide theoretical guarantees, the assumptions involved are hard to satisfy in practice.

The DLAM algorithm {#sec:algorithm}
==================

In this section, we present our novel DLAM algorithm. Section \[sec:problem\] provides the new algorithm’s formulation and Section \[sec:quadratic approxmation\] shows how the DLAM algorithm and quadratic approximation technique function together to solve all the subproblems.

      Notations                                Descriptions
  ----------------- -------------------------------------------------------------------
         $L$                                 Number of layers.
        $W_l$                     The weight vector in the $l$-th layer.
        $b_l$                       The intercept in the $l$-th layer.
        $z_l$        The temporary variable of the linear mapping in the $l$-th layer.
     $h_l(z_l)$           The nonlinear activation function in the $l$-th layer.
        $a_l$                         The output of the $l$-th layer.
         $x$                      The input matrix of the neural network.
         $y$                           The predefined label vector.
     $R(z_l,y)$                   The risk function in the $l$-th layer.
   $\Omega_l(W_l)$             The regularization term in the $l$-th layer.
   $\varepsilon_l$      The tolerance of the nonlinear mapping in the $l$-th layer.

  : Notation Used in This Paper

\[tab:notation\]

Inequality Approximation for Deep Learning {#sec:problem}
------------------------------------------

 The important notation used in this paper is listed in Table \[tab:notation\]. A typical fully-connected deep neural network consists of $L$ layers, each of which are defined by a linear mapping and a nonlinear activation function. A linear mapping is composed of a weight vector $W_l\in \mathbb{R}^{n_l\times n_{l-1}}$ , where $n_l$ is the number of neurons on the $l$-th layer and an intercept $b_l\in \mathbb{R}^{n_l}$; a nonlinear mapping is defined by an activation function $h_l(\bullet)$. Given an input $a_{l-1}\in \mathbb{R}^{n_{l-1}}$ from the $(l-1)$-th layer, the $l$-th layer outputs $a_l=h_l(W_la_{l-1}+b_l)$. By introducing an auxiliary variable $z_l$ as the temporary result of the linear mapping, the deep neural network problem is formulated mathematically as follows:

\[prob: problem 1\] $$\begin{aligned}
     &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \min\nolimits_{a_l,W_l,b_l,z_l} R(z_L;y)+\sum\nolimits_{l=1}^L\Omega_l(W_l)\\
     &s.t.\ z_l\!=W_la_{l-1}\!+\!b_l(l\!=\!1\!,\!\cdots\!,\!L)\!,\ a_l\!=\!h_l\!(z_l)\! \ \!(l\!=\!1,\cdots,L\!-\!1)\end{aligned}$$

where $a_0=x\in\mathbb{R}^d$ is the input of the deep neural network, $d$ is the number of feature dimensions, and $y$ is a predefined label vector. $R(z_L;y)$ is the risk function for the $L$-th layer, which is convex and proper, and $\Omega_l(W_l)$ is a regularization term on the $l$-th layer, which is also convex and proper. The equality constraint $a_l=h_l(z_l)$ is the most challenging to handle here, because common activation functions such as tanh and smooth sigmoid are nonlinear. This makes them nonconvex constraints and hence it is difficult to obtain a global minimum when updating $z_l$ [@taylor2016training]. To deal with this challenge, we innovatively transform the original nonconvex constraints into convex inequality constraints, which can be infinitely approximate to Problem \[prob: problem 1\]. To do this, we introduce a tolerance $\varepsilon_l>0$ and reformulate Problem \[prob: problem 1\] to reach the following form: $$\begin{aligned}
     \min\nolimits_{W_l,b_l,z_l,a_l} &R(z_L;y)\!+\!\sum\nolimits_{l=1}^L\Omega_l(W_l)\\&\!+\!\sum\nolimits_{l\!=\!1}^{L\!-\!1}\mathbb{I}\!(h_l\!(z_l)\!-\!\varepsilon_l\!\leq\! a_l\!\leq\! h_l\!(z_l)\!+\!\varepsilon_l) \\& s.t. z_l=W_la_{l-1}+b_l(l\!=\!1\!,\!\cdots\!,\!L)\end{aligned}$$ $\mathbb{I}(h_l(z_l)-\varepsilon_l\leq a_l\leq h_l(z_l)+\varepsilon_l)$ is an indicator function such that the value is 0 if $h_l(z_l)-\varepsilon_l\leq a_l\leq h_l(z_l)+\varepsilon_l$ and $\infty$ otherwise. For the linear constraint $z_l\!=W_la_{l-1}\!+\!b_l$, this can be transformed into a penalty term in the objective function to minimize the difference between $z_l$ and $W_la_{l-1}\!+\!b_l$. The formulation is shown as follows:

\[prob: problem 2\] $$\begin{aligned}
     &\min\nolimits_{W_l,b_l,z_l,a_l} F(\textbf{W},\textbf{b},\textbf{z},\textbf{a})=R(z_L;y)\!+\!\sum\nolimits_{l=1}^L\Omega_l(W_l)\\&\!+\!\sum\nolimits_{l\!=\!1}^L\!\phi\!(a_{l-1}\!,\!W_l\!,\!b_l\!,\!z_l)\!+\!\sum\nolimits_{l\!=\!1}^{L\!-\!1}\! \mathbb{I}(h_l\!(z_l)\!-\!\varepsilon_l\!\leq\! a_l\!\leq\! h_l\!(z_l)\!+\!\varepsilon_l)\end{aligned}$$

The penalty term is defined as $\phi(a_{l-1},W_l,b_l,z_l)=(\rho/2)\Vert z_l-W_la_{l-1}-b_l\Vert^2_2$, where $\rho>0$ a penalty parameter. $\textbf{W}=\{W_l\}_{l=1}^{L}$, $\textbf{b}=\{b_l\}_{l=1}^{L}$, $\textbf{z}=\{z_l\}_{l=1}^{L}$,$\textbf{a}=\{a_l\}_{l=1}^{L-1}$. The reason for introducing $\varepsilon_l$ is that it allows us to project the nonlinear constraint to a convex $\varepsilon_l$-ball, thus transforming the nonconvex Problem \[prob: problem 1\] into the multi-convex Problem \[prob: problem 2\], which is much easier to solve. Here, a multi-convex problem means this problem is convex with regard to one variable while fixing others. For example, Problem \[prob: problem 2\] is convex with regard to $\textbf{z}$ when $\textbf{W}$, $\textbf{b}$, and $\textbf{a}$ are fixed. As $\rho\rightarrow \infty$ and $\varepsilon_l\rightarrow 0$, Problem \[prob: problem 2\] approaches Problem \[prob: problem 1\].

Quadratically-Approximated Alternative Optimization {#sec:quadratic approxmation}
---------------------------------------------------

In this section, we present the DLAM algorithm developed to solve Problem \[prob: problem 2\], shown in Algorithm \[algo:DLAM\]. Lines 4, 5, 7, and 10 update $W_l$, $b_l$,$z_l$ and $a_l$, respectively, and the four relevant subproblems are discussed in detail below:\

\[algo:DLAM\] $y$, $a_0=x$. $a_l,W_l,b_l, z_l(l=1,\cdots,L)$. Initialize $\rho$, $k=0$. Update $W_l^{k+1}$ using Algorithm \[algo:theta update\]. Update $b^{k+1}_l$ in Equation . Update $z_l^{k+1}$ in Equation . Update $z_l^{k+1}$ in Equation . Update $a^{k+1}_{l}$ using Algorithm \[algo:tau update\]. $k\leftarrow k+1$. convergence. Output $a_l,W_l,b_l,z_l$.

\
**1. Update $W_l$**\
The variables $W_l(l=1,\cdots,L)$ are updated as follows: $$\begin{aligned}
    W^{k+1}_l\leftarrow\arg\min\nolimits_{W_l} \phi(a^{k+1}_{l-1},W_l,b^k_l,z^k_l)+\Omega_l(W_l)
    \label{eq: update W original}\end{aligned}$$ Because $W_l$ and $a_{l-1}$ are coupled in $\phi(\bullet)$, solving $W_l$ requires an inversion operation of $a^{k+1}_{l-1}$, which is computationally expensive. In order to handle this challenge, we define $P^{k+1}_l(W_l;\theta^{k+1}_l)$ as a quadratic approximation of $\phi$ at $W^k_{l}$, which is mathematically reformulated as follows [@beck2009fast]: $$\begin{aligned}
P^{k\!+\!1}_l(W_l\!;\!\theta^{k\!+\!1}_l)&\!=\!\phi(a^{k\!+\!1}_{l\!-\!1}\!,\!W^k_{l}\!,\!z^k_{l}\!,\!b^k_{l})\!+\!<\!\nabla_{W^k_l} \phi\!,\!W_l\!-\!W^k_l\!>\\&+\Vert\theta_l^{k+1}\circ (W_l-W_l^k)^{\circ 2}\Vert_{1}/2\end{aligned}$$ where $\theta_l^{k+1}>0$ is a parameter vector, $\circ$ denotes Hadamard product (the elementwise product), and $a^{\circ b}$ denotes $a$ to the Hadamard power of $b$ and $\Vert\bullet\Vert_{1}$ is the $\ell_{1}$ norm. $<\bullet,\bullet>$ is a Frobenius inner product.$\nabla_{W^k_l}\phi=\rho(W_{l}^ka^{k+1}_{l-1}+b^k_{l}-z^k_{l})(a^{k+1}_{l-1})^T(l=1,\cdots,L)$ is the gradient of $\phi$ with regard to $W_l$ at $W^k_l$. Obviously, $P^{k+1}_l(W^k_l;\theta^{k+1}_l)=\phi(a^{k+1}_{l-1},W^k_{l},b_{l}^k,z^k_{l})$. Rather than minimizing the original subproblem in Equation , we instead minimize the following: $$\begin{aligned}
 &W^{k+1}_l \leftarrow \arg\min\nolimits_{W_l} P^{k+1}_l(W_l;\theta_l^{k+1})+\Omega_l (W_l)     \label{eq:update W}\end{aligned}$$ For $\Omega_l(W_l)$, common regularization terms like $\ell_1$ or $\ell_2$ regularizations lead to closed-form solutions. As for the choice of $\theta_l^{k+1}$, the backtracking algorithm is shown in Algorithm \[algo:theta update\]. Specifically, for a given $\theta^{k+1}_l$, we minimize Equation to obtain $W^{k+1}_l$ until the condition in Line 3 is satisfied. The backtracking algorithm always terminates because as $\theta^{k+1}_l\rightarrow \infty$, $W^{k+1}_l\rightarrow W^{k}_l$, and $W^{k}_l$ satisfies the condition in Line 3. The time complexity of Algorithm \[algo:theta update\] is $O(d^2)$, where $d$ is the dimension of the neurons or features.

\[algo:theta update\]$a_{l-1}^{k+1}$, $W^k_{l}$, $b^k_l$, $z^k_l$, $\rho$, some constant $\gamma>1$. $\theta^{k+1}_l$,$W^{k+1}_l$. Initialize $\alpha$. update $\zeta$ in Equation where $\theta^{k+1}_l=\alpha$. $\alpha\leftarrow \alpha\gamma$.\
update $\zeta$ in Equation where $\theta^{k+1}_l=\alpha$.\
Output $\theta^{k+1}_l \leftarrow \alpha $.\
Output $W^{k+1}_{l}\leftarrow \zeta$.

\
**2. Update $b_l$**\
The variables $b_l(l=1,\cdots,L)$ are updated as follows: $$\begin{aligned}
    b_l^{k+1}\leftarrow \arg\min\nolimits_{b_l} 
     \phi(a^{k+1}_{l-1},W^{k+1}_l,b_l,z^k_l).\end{aligned}$$ The above subproblem has a closed-form solution $b^{k+1}_l=z^{k}_l-W^{k+1}_la^{k+1}_{l-1}$. However, the value of $b^{k+1}_l$ is subject to fluctuation as the signs of either $W_l^{k+1}$ or $a^{k+1}_{l-1}$ may change, which slows down the convergence of $b^{k+1}_l$. We therefore define $U^{k+1}_l(b_l;L_b)$ as a quadratic approximation of $\phi$ at $b^k_l$, which is formulated mathematically as follows [@beck2009fast]: $$\begin{aligned}
    U^{k+1}_l(b_l;L_b)&=\phi(a^{k+1}_{l-1},W^{k+1}_l,b^k_l,z^k_l)+(\nabla_{b^k_l} \phi)^T(b_l-b^k_l)\\&+(L_b/2)\Vert b_l-b^k_l\Vert^2_2.\end{aligned}$$ where $L_b\geq \rho$ is a parameter and $\nabla_{b^k_l} \phi=\rho(b^k_l+W^{k+1}_la^{k+1}_{l-1}-z^k_l)$. Here, $L_b\geq \rho$ is required for the convergence analysis [@beck2009fast]. Without loss of generality, we set $L_b=\rho$. We can now solve the following subproblem: $$\begin{aligned}
    b^{k+1}_l\leftarrow \arg\min\nolimits_{b_l} U^{k+1}_l(b_l;\rho) \label{eq:update b}\end{aligned}$$ The solution to Equation is $$\begin{aligned}
    b^{k+1}_l\leftarrow b^{k}_l-\nabla_{b^k_l} \phi/\rho
    \end{aligned}$$ This indicates that although $b^{k+1}_l$ is closely related to $b^{k}_l$, it is more resistant to the impact of a sign change for either $W^{k+1}_l$ or $a^{k+1}_{l-1}$.\
**3. Update $z_l$**\
The variables $z_l(l=1,\cdots,L)$ are updated as follows: $$\begin{aligned}
      z^{k\!+\!1}_l\! &\leftarrow\arg\min\nolimits_{z_l}\! \phi(a^{k\!+\!1}_{l\!-\!1},W^{k\!+\!1}_l,b^{k\!+\!1}_l,z_l)\\&\!+\!\mathbb{I}(h_l(z_l)\!-\!\varepsilon_l\leq a^{k}_l\leq h_l(z_l)\!+\!\varepsilon_l) (l< \!L)\\
    z^{k+1}_L &\leftarrow\arg\min\nolimits_{z_L} \phi(a^{k+1}_{L-1},W^{k+1}_L,b^{k+1}_L,z_L)+R(z_L;y)\end{aligned}$$ As when updating $b_l$, we define $V^{k+1}_l(z_l;L_z)$ as a quadratic approximation of $\phi$ at $z^k_l$, which is formulated mathematically as follows: $$\begin{aligned}
    V^{k\!+\!1}_l(z_l\!;\!L_z)&=\!\phi(a_{l\!-\!1}^{k\!+\!1},W^{k\!+\!1}_l,b^{k\!+\!1}_l,z^k_l)\!+\!(\nabla_{z^k_l}\phi)^T(z_l\!-\!z^k_l)\\&\!+\!(L_z/2)\Vert z_l\!-\!z^k_l\Vert^2_2\end{aligned}$$ where $L_z\geq \rho$ is a parameter and $\nabla_{z^k_l}\phi=\rho(z_l^k-W^{k+1}_l a^{k+1}_{l-1}-b^{k+1}_l)$. Without loss of generality, we set $L_z=\rho$. Obviously, $V^{k+1}_l(z^k_l;\rho)=\phi(a^{k+1}_{l-1},W^{k+1}_l,b^{k+1}_l,z^k_l)$. Therefore, we solve the following problems: $$\begin{aligned}
 \nonumber z^{k\!+\!1}_l\! &\leftarrow\arg\min\nolimits_{z_l}\! V^{k+1}_l(z_l;\rho)\\& +\!\mathbb{I}(h_l(z_l)\!-\!\varepsilon_l\leq a^{k}_l\leq h_l(z_l)\!+\!\varepsilon_l) (l< \!L) \label{eq:update z}\\
    z^{k+1}_L &\leftarrow\arg\min\nolimits_{z_L} V^{k+1}_L(z_L;\rho)+R(z_L;y) \label{eq:update zl}   \end{aligned}$$ As for $z_l(l=1,\cdots,l-1)$, the solution is $$\begin{aligned}
z_l^{k+1}\leftarrow \min(\max(B^{k+1}_1,z^{k}_l-\nabla\phi_{z^k_l}/\rho), B^{k+1}_2).   \end{aligned}$$ where $B^{k+1}_1$ and $B^{k+1}_2$ represent the lower bound and the upper bound of the set $\{ z_l|h_l(z_l)\!-\!\varepsilon_l\leq a^{k}_l\leq h_l(z_l)\!+\!\varepsilon_l\}$. Equation is easy to solve using the Fast Iterative Soft Thresholding Algorithm (FISTA) [@beck2009fast].\
**4. Update $a_l$**\
The variables $a_l(l=1,\cdots,L-1)$ are updated as follows: $$\begin{aligned}
    a^{k+1}_{l}&\leftarrow \arg\min\nolimits_{a_{l}} \phi(a_{l},W^k_{l+1},b^k_{l+1},z^k_{l+1})\\&+\mathbb{I}(h_l(z^{k+1}_l)-\varepsilon_l\leq a_l\leq h_l(z^{k+1}_l)+\varepsilon_l)\end{aligned}$$ As when solving $W^{k+1}_l$, the quadratic approximation of $\phi$ at $a^k_l$ is defined as $$\begin{aligned}
Q^{k\!+\!1}_l(a_{l};\!\tau^{k\!+\!1}_l)&=\!\phi(a^k_{l}\!,\!W^k_{l\!+\!1}\!,\!b^k_{l\!+\!1}\!,\!z^k_{l\!+\!1})\!+\!(\nabla_{a^k_{l}} \phi)^T(a_{l}\!-\!a^k_{l})\\&+\Vert\tau_l^{k+1}\circ (a_{l}-a_{l}^k)^{\circ 2}\Vert_{1}/2\end{aligned}$$ and this allows us to solve the following problem instead: $$\begin{aligned}
    \nonumber a^{k+1}_l&\leftarrow\arg\min\nolimits_{a_l} Q^{k+1}_l(a_l;\tau^{k+1}_l)\\&+\mathbb{I}(h_l(z^{k+1}_l)-\varepsilon_l\leq a_l\leq h_l(z^{k+1}_l)+\varepsilon_l) \label{eq:update a}\end{aligned}$$ where $\tau_l^{k+1}>0$ is a parameter vector. $\nabla_{a^k_{l}}\phi=\rho(W_{l+1}^{k})^T(W_{l+1}^{k}a^{k}_{l}+b^k_{l+1}-z^k_{l+1})(l=1,\cdots,L-1)$ is the gradient of $\phi$ with regard to $a_l$ at $a^k_l$. Obviously, $Q^{k+1}_l(a^k_{l};\tau^{k+1}_l)=\phi(a^k_{l},W^k_{l+1},b_{l+1}^k,z^k_{l+1})$. Because $Q^{k+1}_l(a_{l};\tau^{k+1}_l)$ is a quadratic function with respect to $a_{l}$, the solution can be obtained by $$\begin{aligned}
a_{l}^{k+1}\leftarrow a^k_{l}-\nabla_{a^k_{l}}\phi/\tau_{l}^{k+1}\end{aligned}$$ given a suitable $\tau_l^{k+1}$. Now the main focus is how to choose $\tau_l^{k+1}$. Similar to Algorithm \[algo:theta update\], the backtracking algorithm for finding a suitable $\tau_l^{k+1}$ is shown in Algorithm \[algo:tau update\]. The time complexity of Algorithm \[algo:tau update\] is $O(d^2)$, where $d$ is the dimension of neurons or features.

\[algo:tau update\]$a_{l}^k$, $W^k_{l+1}$, $z^{k+1}_l$, $z^k_{l+1}$, $b^k_{l+1}$, $\rho$, some constant $\eta>1$. $\tau^{k+1}_l$,$a^{k+1}_{l}$. Pick up $t$ such that $\beta=a^k_l-\nabla_{a^k_l}\phi/t$ and $h_{l}(z^{k+1}_{l})-\varepsilon_l\leq \beta\leq h_{l}(z^{k+1}_{l})+\varepsilon_l$. $t\leftarrow t\eta$.\
$\beta\leftarrow a^k_{l}-\nabla\phi_{a^k_{l}}/t$.\
Output $\tau^{k+1}_l \leftarrow t $.\
Output $a^{k+1}_{l}\leftarrow \beta$.

Convergence Analysis {#sec:convergence}
====================

In this section, we present the main convergence analyses for the DLAM algorithm. Specifically, Section \[sec:assumptions\] introduces the assumptions necessary to guarantee convergence. The main convergence properties of the new DLAM algorithm are presented in Section \[sec:convergence property\].

Assumptions {#sec:assumptions}
-----------

First recall the concept of coercivity as follows [@wang2015global]:

Suppose $f(x_1,\cdots,x_m)$ is a function with respect to $(x_1,\cdots,x_m)\in G$ where $G$ is the domain, then $f(x_1,\cdots,x_m)$ is coercive if $(x_1,\cdots,x_m)\in G$ and $\Vert (x_1,\cdots,x_m)\Vert\rightarrow \infty$ leads to $f\rightarrow \infty$.

Then we propose a new concept called *multi-coercitvity* based on the above definition of coercivity, which is defined as:

Suppose $g(u_1,\cdots,u_m)$ is a function with respect to $(u_1,\cdots,u_m)\in G$ where $G$ is the domain, then $g(u_1,\cdots,u_m)$ is coercive with regard to $u_1$ if $(u_1,\cdots,u_m)\in G$ and $\Vert u_1\Vert\rightarrow \infty$ while fixing $u_i(i=2,3,\cdots,m)$ leads to $g\rightarrow \infty$. If $g$ is coercive with regard to all variables $u_i(i=1,2,\cdots,m)$, then $g$ is multi-coercive.

Multi-coercivity is milder than coercivity because although a coercive function must be multi-coercive, the reverse is not the case. For example, $f_1(x,y)=x+y$ is coercive with regard to $x$ and $y$, respectively, and therefore is multi-coercive. However, $f_1$ is not coercive because when $\Vert (x,y)\Vert\rightarrow\infty$ and $(x,y)$ follows the line $x+y=0$, $f_1=0$. Given this definition of multi-coercivity, we can make the following assumption:

$F(\textbf{W},\textbf{b},\textbf{z},\textbf{a})$ is multi-coercive over the set $S=\{(\textbf{W},\textbf{b},\textbf{z},\textbf{a}):h_l(z_l)-\varepsilon_l\leq a_l\leq h_l(z_l)+\varepsilon_l(l=1\cdots,L-1)\}$. \[ass:assumption 1\]

With the introduction of the multi-coercivity assumption, the condition of our problem becomes much milder than is the case with the solutions proposed by previous researchers and also enables our framework to cover most of the common loss functions utilized in neural networks; cross entropy and square loss are both multi-coercive, for example. The next assumption guarantees that all subproblems have global minima. Before stating the second assumption, recall the definition of quasilinearity [@boyd2004convex]:

A function $f(x)$ is quasiconvex if for any sublevel set $S_\alpha(f)=\{x|f(x)\leq \alpha\}$ is a convex set. Likewise, A function $f(x)$ is quasiconcave if for any superlevel set $S_\alpha(f)=\{x|f(x)\geq \alpha\}$ is a convex set. A function $f(x)$ is quasilinear if it is both quasiconvex and quasiconcave.

Given this definition of quasilinearity, we can make the following assumption:

\[ass:assumption 3\] Activation functions $h_l(z_l)(l=1,\cdots,n)$ are quasilinear functions.

Assumption \[ass:assumption 3\] ensures that the nonlinear constraint $a_l=h_l(z_l)$ in Problem \[prob: problem 1\] is projected in a convex set. Fortunately, most of the widely used nonlinear activation functions, including tanh [@zamanlooy2014efficient], smooth sigmoid [@glorot2010understanding], and the rectified linear unit(Relu) [@maas2013rectifier] that are used in deep neural networks are quasilinear. They therefore fit neatly into our framework and incorporate several important theoretical properties.

Key Convergence Properties {#sec:convergence property}
--------------------------

We introduce several important main convergence properties possessed by the DLAM algorithm in this section. If Assumptions 1 and 2 hold, then Properties 1-3 stated below are satisfied. These three properties are proven to be possessed by the DLAM algorithm, and are key for demonstrating the theoretical merits of DLAM; the proofs of them are provided in the supplementary materials. Finally, the global convergence and convergence rate of the DLAM are proved based on Properties 1-3. The three convergence properties are as follows:

For any $\rho> 0$ and $\varepsilon_l>0$, starting from any $(\textbf{W}^0,\textbf{b}^0,\textbf{z}^0,\textbf{a}^0)$ such that $h_l(z^0_l)-\varepsilon_l\leq a^0_l\leq h_l(z^0_l)+\varepsilon_l(l=1,\cdots,L)$, $\{(\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{a}^k)\}$ is bounded, and $F(\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{a}^k)$ is lower bounded. \[pro:property 1\]

Property \[pro:property 1\] guarantees that all variables and objective functions during iterations are bounded. The proof of Property \[pro:property 1\] requires Lemma \[lemma:lemma 3\] and Assumption \[ass:assumption 1\], and the proof is elaborated in Theorem \[theorem: property 1\] in the supplementary materials.

\[pro:property 2\] For any $\rho>0$ and $\varepsilon_l>0$, we have $$\begin{aligned}
\nonumber &F(\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{a}^k)-F(\textbf{W}^{k+1},\textbf{b}^{k+1},\textbf{z}^{k+1},\textbf{a}^{k+1})\\\nonumber&\geq\sum\nolimits_{l=1}^{L} \Vert\theta_l^{k\!+\!1}\circ (W^{k\!+\!1}_l-W_l^k)^{\circ 2}\Vert_{1}/2 \\\nonumber &+\!(\rho/2)\!\sum\nolimits_{l\!=\!1}^L\Vert b^{k\!+\!1}_l\!-\!b^k_l\Vert^2_2\!+\!(\rho/2)\!\sum\nolimits_{l\!=\!1}^L\Vert z^{k\!+\!1}_l\!-\!z^k_l\Vert^2_2\\&+\sum\nolimits_{l=1}^{L-1} \Vert\tau_l^{k\!+\!1}\circ (a^{k\!+\!1}_l-a_l^k)^{\circ 2}\Vert_{1}/2 \label{eq: property2} \end{aligned}$$

Property \[pro:property 2\] depicts the monotonic decrease of the objective value during iterations. The proof of Property \[pro:property 2\] is detailed in Theorem \[theorem: property 2\] in the supplementary materials.

\[pro:property 3\] There exists a constant $C^{k+1}>0$ and $g\in \partial F(\textbf{W}^{k+1},\textbf{b}^{k+1},\textbf{z}^{k+1},\textbf{a}^{k+1})$ such that $$\begin{aligned}
    \nonumber\Vert g\Vert\! &\leq\! C^{k\!+\!1}\!(\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{a}^{k\!+\!1}\!-\!\textbf{a}^{k}\Vert)\label{eq: property3}\end{aligned}$$

Property \[pro:property 3\] ensures that the subgradient of the objective function is bounded by variables. The proof of Property \[pro:property 3\] requires Property \[pro:property 1\] and the proof process is elaborated in Theorem \[theorem: property 3\] in the supplementary materials. We will now move on to present the global convergence of the DLAM algorithm using the following three theorems. The first theorem states that Properties \[pro:property 1\]-\[pro:property 3\] are guaranteed.

\[thero: theorem 1\] For any $\rho>0$ and $\varepsilon_l>0$, if Assumptions \[ass:assumption 1\] and \[ass:assumption 3\] are satisfied, then Properties \[pro:property 1\]-\[pro:property 3\] hold.

This theorem can be concluded by the proof of Theorems \[theorem: property 1\], \[theorem: property 2\] and \[theorem: property 3\] in the supplementary materials.

The next theorem presents the global convergence of the DLAM algorithm.

\[Global Convergence\] \[thero: theorem 2\] For the variables $(\textbf{W},\textbf{b},\textbf{z},\textbf{a})$ in Problem \[prob: problem 2\], starting from any $(\textbf{W}^0,\textbf{b}^0,\textbf{z}^0,\textbf{a}^0)$ such that $h_l(z^0_l)-\varepsilon_l\leq a^0_l\leq h_l(z^0_l)+\varepsilon_l(l=1,\cdots,L)$ , it has at least a limit point $(\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{a}^*)$, and any limit point $(\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{a}^*)$ is a critical point. That is, $0\in \partial F(\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{a}^*)$.

Because $(\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{a}^k)$ is bounded, there exists a subsequence $(\textbf{W}^s,\textbf{b}^s,\textbf{z}^s,\textbf{a}^s)$ such that $(\textbf{W}^s,\textbf{b}^s,\textbf{z}^s,\textbf{a}^s)\rightarrow (\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{a}^*)$ where $(\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{a}^*)$ is a limit point. By Properties \[pro:property 1\] and \[pro:property 2\], $F(\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{a}^k)$ is non-increasing and lower bounded and hence converged. By Property \[pro:property 2\], we prove that $\Vert \textbf{W}^{k+1}-\textbf{W}^k\Vert\rightarrow 0$, $\Vert \textbf{b}^{k+1}-\textbf{b}^k\Vert\rightarrow 0$, $\Vert \textbf{z}^{k+1}-\textbf{z}^k\Vert\rightarrow 0$ and $\Vert \textbf{a}^{k+1}-\textbf{a}^k\Vert\rightarrow 0$ as $k\rightarrow \infty$ . We infer there exists $g^k\in \partial F(\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{a}^k)$ such that $\Vert g^k\Vert \rightarrow 0$ as $k\rightarrow \infty$ based on Property \[pro:property 3\]. Specifically, $\Vert g^s\Vert \rightarrow 0$ as $s\rightarrow \infty$. According to the definition of general subgradient (Defintion 8.3 in [@rockafellar2009variational]), we have $0\in \partial F(\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{a}^*)$. In other words, the limit point $(\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{a}^*)$ is a critical point of $F$.

Theorem \[thero: theorem 2\] shows that our proposed DLAM algorithm converges globally no matter what $\rho$ and $\varepsilon_l$ are chosen. This ensures that our DLAM algorithm is parameter-restriction free, so the choice of hyperparameters has no effect on its convergence.\
The next theorem shows that the convergence rate of DLAM is $o(1/k)$, which is shown as follows:

For a sequence $(\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{a}^k)$, define $c_k=\min\nolimits_{0\leq i\leq k}(\sum\nolimits_{l=1}^{L} \Vert\theta_l^{i\!+\!1}\circ (W^{i\!+\!1}_l-W_l^i)^{\circ 2}\Vert_{1}/2+\!(\rho/2)\!\sum\nolimits_{l\!=\!1}^L\Vert b^{i\!+\!1}_l\!-\!b^i_l\Vert^2_2\!+\!(\rho/2)\!\sum\nolimits_{l\!=\!1}^L\Vert z^{i\!+\!1}_l\!-\!z^i_l\Vert^2_2+\sum\nolimits_{l=1}^{L-1} \Vert\tau_l^{i\!+\!1}\circ (a^{i\!+\!1}_l-a_l^i)^{\circ 2}\Vert_{1}/2)$, then the convergence rate of $c_k$ is $o(1/k)$. \[thero: theorem 3\]

The proof of this theorem is in Appendix \[sec:convergence rate proof\] in the supplementary materials.

Experiments {#sec:experiment}
===========

In this section, the dlADMM algorithm is evaluated by several benchmark datasets. Effectiveness, efficiency and convergence properties of dlADMM are compared with state-of-the-art methods. All experiments were conducted on 64-bit Ubuntu16.04 LTS with Intel(R) Xeon processor and GTX1080Ti GPU.

Experiment Setup
----------------

### Dataset

In this experiment, two benchmark datasets were used for comparison: MNIST [@lecun1998gradient] and Fashion MNIST [@xiao2017fashion]. The MNIST dataset has ten classes of handwritten-digit images, which was firstly introduced by Lecun et al. in 1998 [@lecun1998gradient]. It contains 55,000 training samples and 10,000 test samples with 196 features each, which is provided by the Keras library [@chollet2017deep]. Unlike the MNIST dataset, the Fashion MNIST dataset has ten classes of assortment images on the website of Zalando, which is Europe’s largest online fashion platform [@xiao2017fashion]. The Fashion-MNIST dataset consists of 60,000 training samples and 10,000 test samples with 784 features each.

### Experiment Settings

We set up two different architectures of multi-layer neural networks in the experiment. Two network structures contained two hidden layers with $100$ and $500$ hidden units each, respectively. The rectified linear unit (Relu) was used for the activation function for both network structures. The loss function was set as the deterministic cross-entropy loss. $\rho$ was set to $10^{-4}$. $\varepsilon$ was initialized as $10$ and updated adaptively as follows: if $R(z^k_L;y)>10\varepsilon^k$, $\varepsilon^{k+1}=\max(2\varepsilon^k,1)$; if $\varepsilon^k>10R(z^k_L;y)$, $\varepsilon^{k+1}=\min(\varepsilon^k/2,0.01)$, which balances between the loss function $R(z_L;y)$ and $\varepsilon$. The number of iteration was set to $150$. In the experiment, one iteration means one epoch.

### Comparison Methods

Since this paper focuses on fully-connected deep neural networks, SGD and its variants and ADMM 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 the accuracy of the training dataset. The baselines are described as follows:\
1. Stochastic Gradient Descent (SGD) [@bottou2010large]. The SGD and its variants are the most popular deep learning optimizers, whose convergence has been studied extensively in the literature.\
2. Adaptive gradient algorithm (Adagrad) [@duchi2011adaptive]. Adagrad is an improved version of SGD: rather than fixing the learning rate during iteration, it adapts the learning rate to the hyperparameter.\
3. Adaptive learning rate method (Adadelta) [@zeiler2012adadelta]. As an improved version of the Adagrad, the Adadelta is proposed to overcome the sensitivity to hyperparameter selection.\
4. Alternating Direction Method of Multipliers (ADMM) [@taylor2016training]. ADMM is a powerful convex optimization method because it can split a objective function into a series of subproblems, which are coordinated to get global solutions. It is scalable to large-scale datasets and supports parallel computations.\

Experimental Results
--------------------

In this section, experimental results of DLAM algorithm are analyzed against comparison methods.

### Convergence

First, we show that our proposed DLAM algorithm converges for both the MNIST dataset and the Fashion MNIST dataset.

![Convergence curves of DLAM algorithm on MNIST and Fashion MNIST datasets for two neural network structures: DLAM algorithm converged.[]{data-label="fig:convergence"}](1a.pdf){width="\columnwidth"}

(a). Objective value on

the $100\times100$ neural network

![Convergence curves of DLAM algorithm on MNIST and Fashion MNIST datasets for two neural network structures: DLAM algorithm converged.[]{data-label="fig:convergence"}](1b.pdf){width="\columnwidth"}

(b). Objective value on

the $500\times500$ neural network

The convergence of DLAM algorithm is shown in Figure \[fig:convergence\]. The X axis and Y axis denote number of iterations and the logorithm of objective value, respectively. Overall, the objective value decreased monotonically during iteration whatever network structures and datasets we choose. Specifically, the objective value droped tremendously at the early stage, and then converged smoothly towards the critical point of the problem. We also found that the objective value for the Fashion MNIST dataset decreased more quickly than that for the MNIST dataset.

### Performance

![Training accuracy of all methods for the MNIST and Fashion MNIST datasets on two neural network structures: DLAM algorithm performed competitively.[]{data-label="fig:training accuracy"}](2a.pdf){width="\textwidth"}

(a). The MNIST dataset

on the $100\times100$ neural network.

![Training accuracy of all methods for the MNIST and Fashion MNIST datasets on two neural network structures: DLAM algorithm performed competitively.[]{data-label="fig:training accuracy"}](2b.pdf){width="\textwidth"}

(b). The Fashion MNIST dataset on

the $100\times100$ neural network.

![Training accuracy of all methods for the MNIST and Fashion MNIST datasets on two neural network structures: DLAM algorithm performed competitively.[]{data-label="fig:training accuracy"}](2c.pdf){width="\textwidth"}

(c). The MNIST dataset

on the $500\times500$ neural network.

![Training accuracy of all methods for the MNIST and Fashion MNIST datasets on two neural network structures: DLAM algorithm performed competitively.[]{data-label="fig:training accuracy"}](2d.pdf){width="\textwidth"}

(d). The Fashion MNIST dataset on

the $500\times500$ neural network.

![Test accuracy of all methods for the MNIST and Fashion MNIST datasets on two neural network structures: DLAM algorithm performed competitively.[]{data-label="fig:test accuracy"}](3a.pdf){width="\textwidth"}

(a). The MNIST dataset

on the $100\times100$ neural network.

![Test accuracy of all methods for the MNIST and Fashion MNIST datasets on two neural network structures: DLAM algorithm performed competitively.[]{data-label="fig:test accuracy"}](3b.pdf){width="\textwidth"}

(b). The Fashion MNIST dataset on

the $100\times100$ neural network.

![Test accuracy of all methods for the MNIST and Fashion MNIST datasets on two neural network structures: DLAM algorithm performed competitively.[]{data-label="fig:test accuracy"}](3c.pdf){width="\textwidth"}

(c). The MNIST dataset

on the $500\times500$ neural network.

![Test accuracy of all methods for the MNIST and Fashion MNIST datasets on two neural network structures: DLAM algorithm performed competitively.[]{data-label="fig:test accuracy"}](3d.pdf){width="\textwidth"}

(d). The Fashion MNIST dataset on

the $500\times500$ neural network.

Figure \[fig:training accuracy\] and Figure \[fig:test accuracy\] show the curves of the training accuracy and test accuracy of our proposed DLAM algorithm and baselines, respectively. Overall, both the training accuracy and the test accuracy of our proposed DLAM outperformed all baselines for the MNIST dataset, while those of our proposed DLAM algortihm performed competitively for the Fashion MNIST dataset. Specifically, the curves of our DLAM algorithm soared to $0.7$ at the early stage, and then raised steadily towards to $0.8$ or more. The curves of the SGD-related methods, SGD, Adadelta, and Adagrad, moved more slowly than our proposed DLAM algorithm. The curves of the ADMM also rocketed to around $0.8$, but decreased slightly later on.

### Efficiency

In this subsection, the relationship between running time per iteration of our proposed DLAM algorithm and two potential factors, namely, the value of $\rho$, the size of training sample was explored. The running time was calculated by the average of 150 iterations. The computational result for the MNIST dataset and Fashion MNIST dataset on the $100\times100$ neural network is shown in Table \[tab:running time\]. The number of training samples of the MNIST dataset ranged from 11,000 to 55,000, with an increase of 11,000 each time, whereas The number of training samples of the Fashion MNIST dataset ranged from 12,000 to 60,000, with an increase of 12,000 each time. The value of $\rho$ ranged from 0.0001 to 1, with multiplying by 10 each time. Generally, the running time increased as the training sample and the value of $\rho$ became larger.

[c|c|c|c|c|c]{}\
&11000&22000 &33000 &44000 &55000\
0.0001& 0.1692& 0.3216& 0.5010& 0.7164& 0.9413\
0.001& 0.2061& 0.4328& 0.6951& 0.9792& 1.2442\
0.01& 0.3334& 0.6516& 1.0277& 1.3956& 1.7783\
0.1& 0.4795& 0.9428& 1.4524& 1.959& 2.4410\
1& 0.7684& 1.4810& 2.2626& 3.0299& 3.7504\
\
&12,000&24,000&36,000&48,000&60,000\
0.0001& 0.2500& 0.5081& 0.8492& 1.1911& 1.5092\
0.001& 0.2980& 0.5980& 0.9595& 1.3265& 1.6744\
0.01& 0.4199& 0.8028& 1.2787& 1.7535& 2.2025\
0.1& 0.5758& 1.0928& 1.7230& 2.3261& 2.9234\
1& 0.8795& 1.6464& 2.5580& 3.4492& 4.2902\

Conclusion
==========

Even though stochastic gradient descent (SGD) is a popular method to train deep neural networks, alternating minimization methods have attracted increasing attention from a great deal of researchers recently as they have several advantages including solid theoretical guarantees and avoiding gradient vanishing problems. In this paper, we propose a novel formulation of the original deep neural network problem and a novel Deep Learning Alternating Minimization (DLAM) algorithm. Specifically, the nonlinear constraint is projected into a convex set so that all subproblems are solvable. At the same time, the quadratic approximation technique and the backtracking algorithm are applied to boost up scalability. Furthermore, several mild assumptions are established to prove the global convergence of our DLAM algorithm. Experiments on real-world datasets demonstrate the convergence, effectiveness, and efficiency of our DLAM algorithm. \[sec:conclusion\]

**Appendix**

Preliminary Lemmas for Proving Three Properties {#sec:proofs}
===============================================

In this section, we give preliminary lemmas which are useful for the proofs of three properties. The proofs of Lemmas \[lemma:lemma 1\] and \[lemma:lemma 2\] both require Assumption \[ass:assumption 3\]. The proof of Lemma \[lemma:lemma 3\] requires Lemmas \[lemma:lemma 1\] and \[lemma:lemma 2\]. To simplify the notation, $\textbf{W}^{k+1}_{\leq l}=\{\{W^{k+1}_i\}_{i=1}^{l},\{W^{k}_i\}_{i=l+1}^{L}\}$, $\textbf{b}^{k+1}_{\leq l}=\{\{b^{k+1}_i\}_{i=1}^{l},\{b^{k}_i\}_{i=l+1}^{L}\}$, $\textbf{z}^{k+1}_{\leq l}=\{\{z^{k+1}_i\}_{i=1}^{l},\{z^{k}_i\}_{i=l+1}^{L}\}$ and $\textbf{a}^{k+1}_{\leq l}=\{\{a^{k+1}_i\}_{i=1}^{l},\{a^{k}_i\}_{i=l+1}^{L-1}\}$.

\[lemma:lemma 1\] Equation holds if and only if there exists $s\in \partial\Omega_l(W^{k+1}_l)$, the subgradient of $\Omega_l(W^{k+1}_l)$ such that $$\begin{aligned}
    \nabla_{W^k_l} \phi +\theta^{k+1}_l \circ(W^{k+1}_l-W^k_l)+s=0\end{aligned}$$ Likewise, Equation holds if and only if there exists $r\in \partial\mathbb{I}(h_l(z^{k+1}_l)-\varepsilon_l\leq a^{k}_l\leq h_l(z^{k+1}_l)+\varepsilon_l)$ such that $$\begin{aligned}
    \nabla_{z^k_l} \phi +\rho(z^{k+1}_l-z^k_l)+r=0\end{aligned}$$ Equation holds if and only if there exists $u\in\partial R(z^{k+1}_L;y)$ such that $$\begin{aligned}
  \nabla_{z^k_L} \phi +\rho(z^{k+1}_L-z^k_L)+u=0  \end{aligned}$$ Equation holds if and only if there exists $v\in \partial\mathbb{I}(h_l(z^{k+1}_l)-\varepsilon_l\leq a^{k+1}_l\leq h_l(z^{k+1}_l)+\varepsilon_l)$ such that $$\begin{aligned}
 \nabla_{a^k_l} \phi +\rho(a^{k+1}_l-a^k_l)+v=0   \end{aligned}$$

These can be obtained by directly applying the optimality conditions of Equation , Equation , Equation and Equation , respectively.

For Equation , Equation and Equation , if $L_b\geq \rho$ and $L_z\geq \rho$,then the following inequalities hold: $$\begin{aligned}
    &U^{k+1}_{l}(b^{k+1}_l;L_b)\geq \phi(a^{k+1}_{l-1},W^{k+1}_l,b^{k+1}_l,z^{k}_l) \label{eq:lipschitz b}\\
    &V^{k+1}_{l}(z^{k+1}_l;L_z)\geq \phi(a^{k+1}_{l-1},W^{k+1}_l,b^{k+1}_l,z^{k+1}_l)
    \label{eq:lipschitz z}\end{aligned}$$ \[lemma:lemma 2\]

Because $\phi(a_{l-1},W_l,b_l,z_l)$ is differentiable continuous with respect to $b_l$ and $z_l$ with Lipschitz coefficient $\rho$ (the definition of Lipschitz differentiablity can be found in [@beck2009fast]), we directly apply Lemma 2.1 in [@beck2009fast] to $\phi$ to obtain Equation and Equation , respectively.

\[lemma:lemma 3\] It holds that for $\forall k\in \mathbb{N}$ and $l=1,2,\cdots,L$, $$\begin{aligned}
\nonumber& F(\textbf{W}^{k\!+\!1}_{\leq l\!-\!1}\!,\!\textbf{b}^{k\!+\!1}_{\leq l\!-\!1}\!,\!\textbf{z}^{k\!+\!1}_{\leq l\!-\!1}\!,\!\textbf{a}^{k\!+\!1}_{\leq l\!-\!1})\!-\!F(\textbf{W}^{k\!+\!1}_{\leq l}\!,\!\textbf{b}^{k\!+\!1}_{\leq l\!-\!1}\!,\!\textbf{z}^{k\!+\!1}_{\leq l\!-\!1},\textbf{a}^{k\!+\!1}_{\leq l\!-\!1})\!\\&\geq\! \Vert\theta_l^{k\!+\!1}\!\circ\! (W^{k\!+\!1}_l-W_l^k)^{\circ 2}\Vert_{1}/2.\label{eq:W optimality}\\
\nonumber &F(\textbf{W}^{k\!+\!1}_{\leq l},\textbf{b}^{k\!+\!1}_{\leq l\!-\!1},\textbf{z}^{k\!+\!1}_{\leq l\!-\!1},\textbf{a}^{k\!+\!1}_{\leq l\!-\!1})\!-\!F(\textbf{W}^{k\!+\!1}_{\leq l},\textbf{b}^{k\!+\!1}_{\leq l},\textbf{z}^{k\!+\!1}_{\leq l\!-\!1},\textbf{a}^{k\!+\!1}_{\leq l\!-\!1})\\&\geq (\rho/2)\Vert b_l^{k+1}-b_l^{k}\Vert^2_2.\label{eq:b optimality}\\
\nonumber&F(\textbf{W}^{k\!+\!1}_{\leq l},\textbf{b}^{k\!+\!1}_{\leq l},\textbf{z}^{k\!+\!1}_{\leq l\!-\!1},\textbf{a}^{k\!+\!1}_{\leq l\!-\!1})\!-\!F(\textbf{W}^{k\!+\!1}_{\leq l},\textbf{b}^{k\!+\!1}_{\leq l},\textbf{z}^{k\!+\!1}_{\leq l},\textbf{a}^{k\!+\!1}_{\leq l\!-\!1})\!\\&\geq\! (\rho/2)\Vert z_l^{k\!+\!1}\!-\!z_l^k\Vert^2_2. \label{eq:z optimality}\\
\nonumber &F(\textbf{W}^{k\!+\!1}_{\leq l},\textbf{b}^{k\!+\!1}_{\leq l},\textbf{z}^{k\!+\!1}_{\leq l},\textbf{a}^{k\!+\!1}_{\leq l-1})\!-\!F(\textbf{W}^{k\!+\!1}_{\leq l},\textbf{b}^{k\!+\!1}_{\leq l},\textbf{z}^{k\!+\!1}_{\leq l},\textbf{a}^{k\!+\!1}_{\leq l})\\&\geq \Vert\tau_l^{k\!+\!1}\circ (a^{k\!+\!1}_l-a_l^k)^{\circ 2}\Vert_{1}/2.
\label{eq:a optimality}\end{aligned}$$

Essentially, all inequalities can be obtained by applying optimality conditions of updating $W^{k+1}_l$, $b^{k+1}_l$, $z^{k+1}_l$ and $a^{k+1}_l$, respectively. We only prove Equation and Equation since Equation and Equation follow the same routine of Equation and Equation , respectively.\
Firstly, we focus on proving Equation . The stopping criterion of Algorithm \[algo:theta update\] shows that $$\begin{aligned}
    \phi(a^{k+1}_{l-1},W^{k+1}_l,b^k_l,z^k_l)\leq P^{k+1}_l(W^{k+1}_l;\theta^{k+1}_l).\label{ineq:ineq1}\end{aligned}$$ Because $\Omega_{W_l}(W_l)$ is convex, according to the definition of subgradient, we have $$\begin{aligned}
  \Omega_l(W^k_{l})\geq  \Omega_l(W^{k+1}_l)+s^T(W_l^{k}-W_l^{k+1}) \label{ineq:ineq2}\end{aligned}$$ where $s$ is defined in the premise of Lemma \[lemma:lemma 1\]. Therefore, we have $$\begin{aligned}
    &F(\textbf{W}^{k\!+\!1}_{\leq l\!-\!1},\textbf{b}^{k\!+\!1}_{\leq l\!-\!1},\textbf{z}^{k\!+\!1}_{\leq l\!-\!1},\textbf{a}^{k\!+\!1}_{\leq l\!-\!1})\!-\!F(\textbf{W}^{k\!+\!1}_{\leq l},\textbf{b}^{k\!+\!1}_{\leq l\!-\!1},\textbf{z}^{k\!+\!1}_{\leq l\!-\!1},\textbf{a}^{k\!+\!1}_{\leq l\!-\!1})\\&=\phi(a^{k+1}_{l-1},W^{k}_l,b^k_l,z^k_l)+\Omega_l(W^k_l)-\phi(a^{k+1}_{l-1},W^{k+1}_l,b^k_l,z^k_l)\\&-\Omega_l(W^{k+1}_{l})\ \text{(Definition of $F$ in Problem \ref{prob: problem 2})}\\&\geq \Omega_l(W^k_l)-\Omega_l(W^{k+1}_{l})-(\nabla_{W^k_l}\phi)^T(W^{k+1}_l-W^k_l)\\&-\Vert\theta_l^{k+1}\circ (W^{k+1}_l-W_l^k)^{\circ 2}\Vert_{1}/2(\text{ Equation \eqref{ineq:ineq1}})\\&\geq s^T(W^{k}_l-W^{k+1}_l)-(\nabla_{W^k_l} \phi)^T(W^{k+1}_l-W^k_l)-\\&\Vert\theta_l^{k+1}\circ (W^{k+1}_l-W_l^k)^{\circ 2}\Vert_{1}/2(\text{ Equation \eqref{ineq:ineq2}})\\ &=(s\!+\!\nabla\! \phi^T_{W^k_l})(W^k_l\!-\!W^{k\!+\!1}_l)\!-\!\Vert\theta_l^{k\!+\!1}\circ (W^{k\!+\!1}_l\!-\!W_l^k)^{\circ 2}\Vert_{1}/2\\&=\Vert\theta_l^{k\!+\!1}\circ (W^{k\!+\!1}_l\!-\!W_l^k)^{\circ 2}\Vert_{1}/2\ (\text{Lemma \ref{lemma:lemma 1}}).\end{aligned}$$ Secondly, we focus on proving Equation . For $l<L$, because $\mathbb{I}(h_l(z_l)-\varepsilon_l\leq a^k_l\leq h_l(z_l)+\varepsilon_l)$ is convex with regard to $z_l$, according to the definition of subgradient, we have $$\begin{aligned}
\nonumber &\mathbb{I}(h_l(z^{k}_l\!)\!-\!\varepsilon_l\leq a^k_l\leq h_l\!(z^{k}_l)\!+\!\varepsilon_l)\\&\geq  \!\mathbb{I}(h_l\!(z^{k\!+\!1}_l)\!-\!\varepsilon_l\!\leq\! a^k_l\!\leq\! h_l\!(z^{k\!+\!1}_l)\!+\!\varepsilon_l)\!+\!r^T\!(z_l^{k}\!-\!z_l^{k\!+\!1})   \label{ineq:ineq3}\end{aligned}$$ where $r$ is defined in Lemma \[lemma:lemma 1\]. $$\begin{aligned}
&F(\textbf{W}^{k+1}_{\leq l},\textbf{b}^{k+1}_{\leq l},\textbf{z}^{k+1}_{\leq l-1},\textbf{a}^{k+1}_{\leq l-1})-F(\textbf{W}^{k+1}_{\leq l},\textbf{b}^{k+1}_{\leq l},\textbf{z}^{k+1}_{\leq l},\textbf{a}^{k+1}_{\leq l-1})\\&=\phi(a^{k\!+\!1}_{l\!-\!1},W^{k\!+\!1}_l,b^{k\!+\!1}_l,z^k_l)\!+\!\mathbb{I}(h_l(z^k_l)\!-\!\varepsilon_l\leq a^k_l\!\leq\! h_l(z^k_l)\!+\!\varepsilon_l)\\&\!-\! \phi\!(a^{k\!+\!1}_{l\!-\!1},W^{k\!+\!1}_l,b^{k\!+\!1}_l,z^{k\!+\!1}_l)\!-\!\mathbb{I}(h_l\!(z^{k\!+\!1}_l)\!-\!\varepsilon_l\!\leq\! a^k_l\!\leq\! h_l(z^{k\!+\!1}_l)\!+\!\varepsilon_l)\\&\text{(Definition of $F$ in Problem \ref{prob: problem 2})}\\&\geq \!-(\nabla \!_{z^k_l}\phi)^T\!(z^{k\!+\!1}_l\!-\!z^k_l)\!-\!(\rho/2)\Vert z^{k\!+\!1}_l\!-\!z^k_l\Vert^2_2\\&+\!\mathbb{I}(h_l\!(z^k_l)\!-\!\varepsilon_l\!\leq\! a^k_l\!\leq\! h_l\!(z^k_l)\!+\!\varepsilon_l)\!-\!\mathbb{I}\!(h_l\!(z^{k+1}_l)\!-\!\varepsilon_l\!\leq\! a^k_l\!\leq\! h_l\!(z^{k\!+\!1}_l)\!+\!\varepsilon_l)\\&\text{(Equation \eqref{eq:lipschitz z})}\\&\geq -(\nabla_{z^k_l} \phi)^T(z^{k+1}_l-z^k_l)-(\rho/2)\Vert z^{k+1}_l-z^k_l\Vert^2_2+r^T(z^k_l-z^{k+1}_l) \\&\text{(Equation \eqref{ineq:ineq3})}\\&=-(\nabla_{z^k_l} \phi)^T(z^{k+1}_l-z^k_l)-(\rho/2)\Vert z^{k+1}_l-z^k_l\Vert^2_2\\&+(\nabla_{z^k_l}\phi+\rho(z_l^{k+1}-z_l^k))^T(z^{k+1}_l-z_l^k) \ \text{(Lemma \ref{lemma:lemma 1})}\\&=(\rho/2)\Vert z_l^{k+1}-z_l^k\Vert^2_2.\end{aligned}$$ For $z_L$, the same routine applies.

Proof of Theorem \[thero: theorem 1\]
=====================================

Proving Theorem \[thero: theorem 1\] is equal to proving jointly Theorems \[theorem: property 1\], \[theorem: property 2\] and \[theorem: property 3\], which are elaborated in the following.

\[theorem: property 1\] Given that Assumptions \[ass:assumption 1\] and \[ass:assumption 3\] hold, the DLAM satisfies Property \[pro:property 1\].

The lower boundness of $F(\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{a}^k)$ can be directly obtained by Lemma \[lemma:lemma 3\]. Moreover, we set $l=L$ and $k:=k+1$ in Equation to get $$\begin{aligned}
F(\textbf{W}^k,\!\textbf{b}^k,\!\textbf{z}^k,\!\textbf{a}^k)\!\leq\! F(\textbf{W}^{k} ,\!\textbf{b}^{k},\!\textbf{z}^{k}_{\leq L-1},\!\textbf{a}^{k})\!-\!(\rho/2)\Vert z^k_L\!-\!z^{k-1}_L\Vert^2_2.\end{aligned}$$ Because $F(\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{a}^k)$ is multi-coercive and $F(\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{a}^k)$ is upper bounded by $F(\textbf{W}^{k} ,\textbf{b}^{k},\textbf{z}^{k}_{\leq L-1},\textbf{a}^{k})-(\rho/2)\Vert z^k_L-z^{k-1}_L\Vert^2_2$, $(\textbf{W}^k,\textbf{b}^k,\textbf{z}^k,\textbf{a}^k)$ is also bounded.

\[theorem: property 2\] Given that Assumptions \[ass:assumption 1\] and \[ass:assumption 3\] hold, the DLAM satisfies Property \[pro:property 2\].

This can be obtained by adding Equation , Equation and Equation from $l=1$ to $l=L$ and Equation from $l=1$ to $l=L-1$.

\[theorem: property 3\] Given that Assumptions \[ass:assumption 1\] and \[ass:assumption 3\] hold, the DLAM satisfies Property \[pro:property 3\].

We know that $
    \partial F(\textbf{W}^{k+1}_1,\textbf{b}^{k+1},\textbf{z}^{k+1},\textbf{a}^{k+1})=(\partial_{\textbf{W}^{k+1}} F,\nabla_{\textbf{b}^{k+1}} F,\partial_{\textbf{z}^{k+1}} F,\partial_{\textbf{a}^{k+1}} F)
$ where $\partial_{\textbf{W}^{k+1}} F=\{\partial_{W^{k+1}_l} F\}_{l=1}^L$, $\nabla_{\textbf{b}^{k+1}} F=\{\nabla_{b^{k+1}_l} F\}_{l=1}^L$, $\partial_{\textbf{z}^{k+1}} F=\{\partial_{z^{k+1}_l} F\}_{l=1}^L$ and $\partial_{\textbf{a}^{k+1}} F=\{\partial_{a^{k+1}_l} F\}_{l=1}^{L-1}$. To prove Property \[pro:property 3\], we need to give a upper bound of $\partial_{W^{k+1}_l} F$, $\partial_{b^{k+1}_l} F$, $\partial_{z^{k+1}_l} F$ and $\partial_{a^{k+1}_l} F$ by a linear combination of $\Vert W_l^{k+1}-W_l^k\Vert$, $\Vert b_l^{k+1}-b_l^k\Vert$, $\Vert z_l^{k+1}-z_l^k\Vert$ and $\Vert a_l^{k+1}-a_l^k\Vert$.\
For $W^{k+1}_l$, $$\begin{aligned}
    \partial_{W^{k+1}_l} &F\\&=\partial\Omega_l(W^{k+1}_l)\!+\!\nabla_{W_l^{k\!+\!1}} \phi(a^{k\!+\!1}_{l\!-\!1},W^{k\!+\!1}_l,b^{k\!+\!1}_l,z^{k\!+\!1}_l)\\&\text{(Definition of $F$ in Equation \eqref{prob: problem 2})}\\&=\nabla_{W_l^{k+1}} \phi(a^{k+1}_{l-1},W^{k+1}_l,b^{k+1}_l,z^{k+1}_l)\\&-\nabla_{W_l^{k}} \phi(a^{k+1}_{l-1},W^{k}_l,b^{k}_l,z^{k}_l)-\theta^{k+1}_l\circ(W^{k+1}_l-W^{k}_l)\\&+\partial\Omega_l(W^{k+1}_l)+\nabla_{W_l^{k}} \phi(a^{k+1}_{l-1},W^{k}_l,b^{k}_l,z^{k}_l)\\&+\theta^{k+1}_l\circ(W^{k+1}_l-W^{k}_l)\\
    &=\rho(W^{k+1}_l-W^{k}_l)a^{k+1}_{l-1}(a^{k+1}_{l-1})^T+\rho(b^{k+1}_l-b^k_l)(a^{k+1}_{l-1})^T\\&-\rho(z^{k+1}_l-z^k_l)(a^{k+1}_{l-1})^T-\theta^{k+1}_l\circ(W^{k+1}_l-W^k_l)\\&+\partial\Omega_l(W^{k+1}_l)+\nabla_{W_l^{k}} \phi(a^{k+1}_{l-1},W^{k}_l,b^{k}_l,z^{k}_l)\\&+\theta^{k+1}_l\circ(W^{k+1}_l-W^{k}_l)\end{aligned}$$ Because $$\begin{aligned}
&\Vert\rho(W^{k+1}_l-W^{k}_l)a^{k+1}_{l-1}(a^{k+1}_{l-1})^T+\rho(b^{k+1}_l-b^k_l)(a^{k+1}_{l-1})^T\\&-\rho(z^{k+1}_l-z^k_l)(a^{k+1}_{l-1})^T-\theta^{k+1}_l\circ(W^{k+1}_l-W^k_l)\Vert\\
&\leq\rho\Vert(W^{k+1}_l-W^{k}_l)a^{k+1}_{l-1}(a^{k+1}_{l-1})^T\Vert+\rho\Vert(b^{k+1}_l-b^k_l)(a^{k+1}_{l-1})^T\Vert\\&\!+\!\rho\!\Vert(z^{k\!+\!1}_l\!-\!z^k_l)\!(\!a^{k\!+\!1}_{l\!-\!1})^T\!\Vert\!+\!\Vert\theta^{k\!+\!1}_l\!\circ\!(W^{k\!+\!1}_l\!-\!W^k_l)\Vert \text{(triangle inequality)}\\&\leq \rho\Vert W^{k+1}_l-W^{k}_l\Vert\Vert a^{k+1}_{l-1}\Vert\Vert a^{k+1}_{l-1}\Vert+\rho\Vert b^{k+1}_l-b^k_l\Vert\Vert a^{k+1}_{l-1}\Vert\\&+\rho\Vert z^{k+1}_l-z^k_l\Vert\Vert a^{k+1}_{l-1}\Vert+\Vert\theta^{k+1}_l\Vert\Vert W^{k+1}_l-W^k_l\Vert\\& \text{(Cauchy-Schwarz inequality)}\end{aligned}$$ and the optimality condition of Equation yields $$\begin{aligned}
    0\in \partial\Omega_l(W^{k\!+\!1}_l)\!+\!\nabla_{W_l^{k}}\! \phi(a^{k\!+\!1}_{l\!-\!1},W^{k}_l,b^{k}_l,z^{k}_l)\!+\!\theta^{k\!+\!1}_l\!\circ\!(W^{k\!+\!1}_l\!-\!W^{k}_l)\end{aligned}$$ Because $a^{k+1}_{l-1}$ is bounded by Property \[pro:property 1\], $\Vert\partial_{W_l^{k+1}} F\Vert$ can be upper bounded by a linear combination of $\Vert W^{k+1}_l-W^k_l\Vert$, $\Vert b^{k+1}_l-b^k_l\Vert$ and $\Vert z^{k+1}_l-z^k_l\Vert$.\
For $b^{k+1}_l$, $$\begin{aligned}
    \nabla F_{b^{k+1}_l}&=\nabla_{b^{k+1}_l} \phi(a^{k+1}_{l-1},W^{k+1}_{l-1},b^{k+1}_l,z^{k+1}_l)\\&=\nabla_{b^{k\!+\!1}_l}\! \phi(a^{k\!+\!1}_{l\!-\!1},W^{k\!+\!1}_{l\!-\!1},b^{k\!+\!1}_l,z^{k\!+\!1}_l)\\&-\nabla_{b^{k}_l} \phi(a^{k+1}_{l-1},W^{k+1}_{l-1},b^{k}_l,z^{k}_l)-\rho(b^{k+1}_l-b^k_l)\\&\text{($\nabla_{b^{k}_l} \phi(a^{k+1}_{l-1},W^{k+1}_{l-1},b^{k}_l,z^{k}_l)+\rho(b^{k+1}_l-b^k_l))=0$ by} \\&\text{the optimality condition of Equation \eqref{eq:update b})}\\&= \rho(z^{k+1}_l-z^{k}_l).
    \end{aligned}$$ Therefore, $\Vert\nabla_{b^{k+1}_l} F\Vert$ is linearly independent on $\Vert z^{k+1}_l-z^k_l\Vert$.\
For $z^{k+1}_l(l< L)$, $$\begin{aligned}
    &\partial_{z^{k+1}_l} F\\&=\nabla_{z^{k+1}_l} \phi(a^{k+1}_{l-1},W^{k+1}_{l-1},b^{k+1}_l,z^{k+1}_l)\\&+\partial \mathbb{I}(h_l(z^{k+1}_l)-\varepsilon_l\leq a^{k+1}_l\leq h_l(z^{k+1}_l)+\varepsilon_l)\\&= \nabla\!_{z^{k\!+\!1}_l} \phi(a^{k\!+\!1}_{l\!-\!1},W^{k\!+\!1}_{l\!-\!1},b^{k\!+\!1}_l,z^{k\!+\!1}_l)\!-\!\nabla\!_{z^{k}_l} \phi(a^{k\!+\!1}_{l\!-\!1},W^{k\!+\!1}_{l\!-\!1},b^{k\!+\!1}_l,z^{k}_l)\!\\&-\!\rho(z^{k+1}_l-z^k_l)+\partial \mathbb{I}(h_l(z^{k+1}_l)-\varepsilon_l\leq a^{k+1}_l\leq h_l(z^{k+1}_l)+\varepsilon_l)\\&\!-\!\partial \mathbb{I}\!(h_l\!(z^{k\!+\!1}_l)\!-\!\varepsilon_l\!\leq\! a^{k}_l\!\leq\! h_l\!(z^{k\!+\!1}_l)\!+\!\varepsilon_l)\!+\! \nabla\!_{z^{k}_l} \phi\!(a^{k\!+\!1}_{l\!-\!1},W^{k\!+\!1}_{l\!-\!1},b^{k\!+\!1}_l,z^{k}_l)\\&+\partial \mathbb{I}(h_l(z^{k+1}_l)-\varepsilon_l\leq a^{k}_l\leq h_l(z^{k+1}_l)+\varepsilon_l)+\rho(z^{k+1}_l-z^k_l)  \\&=\nabla\!_{z^{k\!+\!1}_l} \phi(a^{k\!+\!1}_{l\!-\!1},W^{k\!+\!1}_{l\!-\!1},b^{k\!+\!1}_l,z^{k\!+\!1}_l)\!-\!\nabla_{z^{k}_l} \phi(a^{k\!+\!1}_{l\!-\!1},W^{k+1}_{l-1},b^{k+1}_l,z^{k}_l)\\&-\rho(z^{k+1}_l-z^k_l) \\&\text{($0\in \partial \mathbb{I}(\bullet)$ and $0\!\in\!\nabla\!_{z^{k}_l} \phi\!(a^{k\!+\!1}_{l\!-\!1},W^{k\!+\!1}_{l\!-\!1},b^{k\!+\!1}_l,z^{k}_l)\!+\rho(z_{l}^{k\!+\!1}\!-\!z_{l}^k)$}\\&\text{$+\!\partial \mathbb{I}(h_l(z^{k\!+\!1}_l)\!-\!\varepsilon_l\leq a^{k}_l\leq h_l(z^{k\!+\!1}_l)\!+\!\varepsilon_l)$ by Equation \eqref{eq:update z})}\\&=0
    \end{aligned}$$ Likewise, we can prove $\partial_{z^{k+1}_L}  F=0$ by the optimality condition of Equation .\
For $a^{k+1}_{l}$, $$\begin{aligned}
        &\partial_{a^{k+1}_l} F\\&=\nabla_{a^{k+1}_l} \phi(a^{k+1}_l,W^{k+1}_{l+1},z^{k+1}_{l+1},b^{k+1}_{l+1})\\&+\partial\mathbb{I}(h_l(z^{k+1}_l)-\varepsilon_l\leq a^{k+1}_l\leq h_l(z^{k+1}_l)+\varepsilon_l)\\&=\nabla_{a^{k+1}_l} \phi(a^{k+1}_l,W^{k+1}_{l+1},z^{k+1}_{l+1},b^{k+1}_{l+1})\\&+\partial\mathbb{I}(h_l(z^{k+1}_l)-\varepsilon_l\leq a^{k+1}_l\leq h_l(z^{k+1}_l)+\varepsilon_l)\\&-\nabla_{a^{k}_l} \phi(a^{k}_l,W^{k}_{l+1},z^{k}_{l+1},b^{k}_{l+1})-\tau^{k+1}_l\circ(a^{k+1}_l-a^k_l)\\&+\nabla_{a^{k}_l} \phi(a^{k}_l,W^{k}_{l+1},z^{k}_{l+1},b^{k}_{l+1})+\tau^{k+1}_l\circ(a^{k+1}_l-a^k_l)\\&=(W^{k+1}_{l+1})^T(W^{k+1}_{l+1}a^{k+1}_l-z^{k+1}_{l+1}-b^{k+1}_{l+1})\\&-(W^{k}_{l+1})^T(W^{k}_{l+1}a^{k}_l-z^{k}_{l+1}-b^{k}_{l+1})-\tau^{k+1}_l\circ(a^{k+1}_l-a^k_l)\\&+\nabla_{a^{k}_l} \phi(a^{k}_l,W^{k}_{l+1},z^{k}_{l+1},b^{k}_{l+1})+\tau^{k+1}_l\circ(a^{k+1}_l-a^k_l)\\&+\partial\mathbb{I}(h_l(z^{k+1}_l)-\varepsilon_l\leq a^{k+1}_l\leq h_l(z^{k+1}_l)+\varepsilon_l)
    \end{aligned}$$ Because $$\begin{aligned}
        &\Vert (W^{k\!+\!1}_{l\!+\!1})^T(W^{k\!+\!1}_{l\!+\!1}a^{k\!+\!1}_l\!-\!z^{k\!+\!1}_{l\!+\!1}\!-\!b^{k\!+\!1}_{l\!+\!1})\!-\!(W^{k}_{l\!+\!1})^T(W^{k}_{l\!+\!1}a^{k}_l\!-\!z^{k}_{l\!+\!1}\!-\!b^{k}_{l\!+\!1})\!\\&-\!\tau^{k\!+\!1}_l\circ(a^{k\!+\!1}_l\!-\!a^k_l)\Vert\\&=\Vert ((W^{k+1}_{l+1})^T W^{k+1}_{l+1}a^{k+1}_{l}-(W^{k}_{l})^T W^{k}_{l+1}a^{k}_{l})\\&+( (W^{k+1}_{l+1})^Tz^{k+1}_{l+1}-(W^{k}_{l+1})^Tz^{k}_{l+1})\\&+\!( (W^{k\!+\!1}_{l\!+\!1})^Tb^{k\!+\!1}_{l\!+\!1}\!-\!(W^{k}_{l\!+\!1})^Tb^{k}_{l\!+\!1})\!+\!( \tau^{k\!+\!1}_l\!\circ(a^{k\!+\!1}_l\!-\!a^k_l))\Vert\\&\leq\Vert (W^{k+1}_{l+1})^T W^{k+1}_{l+1}a^{k+1}_{l}-(W^{k}_{l})^T W^{k}_{l+1}a^{k}_{l}\Vert\\&+\Vert (W^{k+1}_{l+1})^Tz^{k+1}_{l+1}-(W^{k}_{l+1})^Tz^{k}_{l+1}\Vert\\&+\Vert (W^{k+1}_{l+1})^Tb^{k+1}_{l+1}-(W^{k}_{l+1})^Tb^{k}_{l+1}\Vert+\Vert \tau^{k+1}_l\circ(a^{k+1}_l-a^k_l)\Vert\\& \text{(triangle inequality)}\\&=\Vert(W^{k\!+\!1}_{l\!+\!1}\!-\!W^k_{l\!+\!1})^TW^{k\!+\!1}_{l\!+\!1}a^{k\!+\!1}_l\!+\!(W^k_{l\!+\!1})^T(W^{k+1}_{l+1}-W^k_{l+1})a^{k+1}_l\\&+(W^k_{l+1})^TW^k_{l+1}(a^{k+1}_l-a^{k}_l)\Vert\\&+\Vert (W^{k+1}_{l+1})^T(z_{l+1}^{k+1}-z_{l+1}^k)+(W_{l+1}^{k+1}-W_{l+1}^k)^Tz_{l+1}^k\Vert\\&+\Vert (W^{k+1}_{l+1})^T(b_{l+1}^{k+1}-b_{l+1}^k)+(W_{l+1}^{k+1}-W_{l+1}^k)^Tb_{l+1}^k\Vert\\&+\Vert \tau^{k+1}_l\circ(a^{k+1}_l-a^k_l)\Vert\\&\leq \Vert(W^{k\!+\!1}_{l\!+\!1}\!-\!W^k_{l\!+\!1})^TW^{k\!+\!1}_{l\!+\!1}a^{k\!+\!1}_l\Vert\!+\!\Vert(W^k_{l\!+\!1})^T(W^{k\!+\!1}_{l\!+\!1}-W^k_{l\!+\!1})a^{k\!+\!1}_l\Vert\\&+\Vert(W^k_{l\!+\!1})^TW^k_{l\!+\!1}(a^{k\!+\!1}_l\!-\!a^{k}_l)\Vert\!+\!\Vert (W^{k+1}_{l+1})^T(z_{l+1}^{k+1}-z_{l+1}^k)\Vert\\&+\Vert(W_{l+1}^{k+1}-W_{l+1}^k)^Tz_{l+1}^k\Vert+\Vert (W^{k+1}_{l+1})^T(b_{l+1}^{k+1}-b_{l+1}^k)\Vert\\&+\Vert(W_{l+1}^{k+1}-W_{l+1}^k)^Tb_{l+1}^k\Vert+\Vert \tau^{k+1}_l\circ(a^{k+1}_l-a^k_l)\Vert\\& (\text{triangle inequality})\\&\leq \Vert W^{k\!+\!1}_{l+1}\!-\!W^k_{l\!+\!1}\Vert\Vert W^{k\!+\!1}_{l\!+\!1}\Vert a^{k\!+\!1}_l\Vert\!+\!\Vert W^k_{l\!+\!1}\Vert\Vert W^{k\!+\!1}_{l\!+\!1}-W^k_{l\!+\!1}\Vert \Vert a^{k\!+\!1}_l\Vert\\&+\Vert W^k_{l+1}\Vert \Vert W^k_{l+1}\Vert \Vert a^{k+1}_l-a^{k}_l\Vert+\Vert W^{k+1}_{l+1}\Vert\Vert z_{l+1}^{k+1}-z_{l+1}^k\Vert\\&+\Vert W_{l+1}^{k+1}-W_{l+1}^k\Vert \Vert z_{l+1}^k\Vert+\Vert W^{k+1}_{l+1}\Vert\Vert b_{l+1}^{k+1}-b_{l+1}^k\Vert\\&+\Vert W_{l+1}^{k+1}-W_{l+1}^k\Vert \Vert b_{l+1}^k\Vert+\Vert \tau^{k+1}_l\Vert\Vert a^{k+1}_l-a^k_l\Vert \\& (\text{Cauchy-Schwarz inequality})
        \end{aligned}$$ and the optimality condition of Equation yields $$\begin{aligned}
         \\&0\in \nabla_{a^{k}_l} \phi(a^{k}_l,W^{k}_{l+1},z^{k}_{l+1},b^{k}_{l+1})+\tau^{k+1}_l\circ(a^{k+1}_l-a^k_l)\\&+\partial\mathbb{I}(h_l(z^{k+1}_l)-\varepsilon_l\leq a^{k+1}_l\leq h_l(z^{k+1}_l)+\varepsilon_l)
     \end{aligned}$$ Because $a^{k+1}_{l}$,$W^k_{l+1}$,$W^{k+1}_{l+1}$, $z^k_{l+1}$ and $b^k_{l+1}$ are bounded, $\Vert\partial_{a_l^{k+1}} F\Vert$ can be upper bounded by a linear combination of $\Vert W^{k+1}_l-W^k_l\Vert$,$\Vert b^{k+1}_l-b^k_l\Vert$,$\Vert z^{k+1}_l-z^k_l\Vert$ and $\Vert a^{k+1}_l-a^k_l\Vert$.\

Proof of Theorem \[thero: theorem 3\] {#sec:convergence rate proof}
=====================================

To prove this theorem, we will first show that $c_k$ satisfies two conditions: (1). $c_k\geq c_{k+1}$. (2). $\sum\nolimits_{k=0}^\infty c_k$ is bounded. We then conclude the convergence rate of $o(1/k)$ based on these two conditions. Specifically, first, we have $$\begin{aligned}
    c_k&\!=\!\min\nolimits_{0\leq i\leq k}\!(\sum\nolimits_{l=1}^{L} \Vert\theta_l^{i\!+\!1}\circ (W^{i\!+\!1}_l-W_l^i)^{\circ 2}\Vert_{1}/2 \\ &+\!(\rho/2)\!\sum\nolimits_{l\!=\!1}^L\Vert b^{i\!+\!1}_l\!-\!b^i_l\Vert^2_2\!+\!(\rho/2)\!\sum\nolimits_{l\!=\!1}^L\Vert z^{i\!+\!1}_l\!-\!z^i_l\Vert^2_2\\&+\sum\nolimits_{l=1}^{L-1} \Vert\tau_l^{i\!+\!1}\circ (a^{i\!+\!1}_l-a_l^i)^{\circ 2}\Vert_{1}/2)\\&\geq\!\min\nolimits_{0\leq i\leq k+1}\!(\sum\nolimits_{l=1}^{L} \Vert\theta_l^{i\!+\!1}\circ (W^{i\!+\!1}_l-W_l^i)^{\circ 2}\Vert_{1}/2 \\ &+\!(\rho/2)\!\sum\nolimits_{l\!=\!1}^L\Vert b^{i\!+\!1}_l\!-\!b^i_l\Vert^2_2\!+\!(\rho/2)\!\sum\nolimits_{l\!=\!1}^L\Vert z^{i\!+\!1}_l\!-\!z^i_l\Vert^2_2\\&+\sum\nolimits_{l=1}^{L-1} \Vert\tau_l^{i\!+\!1}\circ (a^{i\!+\!1}_l-a_l^i)^{\circ 2}\Vert_{1}/2)\\&= c_{k+1}\end{aligned}$$ Therefore $c_k$ satisfies the first condition. Second, $$\begin{aligned}
    &\sum\nolimits_{k\!=\!0}^\infty c_k\\\!&=\sum\nolimits_{k=0}^\infty\min\nolimits_{0\leq i\leq k}\!(\sum\nolimits_{l=1}^{L} \Vert\theta_l^{i\!+\!1}\circ (W^{i\!+\!1}_l-W_l^i)^{\circ 2}\Vert_{1}/2 \\&+\!(\rho/2)\!\sum\nolimits_{l\!=\!1}^L\Vert b^{i\!+\!1}_l\!-\!b^i_l\Vert^2_2\!+\!(\rho/2)\!\sum\nolimits_{l\!=\!1}^L\Vert z^{i\!+\!1}_l\!-\!z^i_l\Vert^2_2\\&+\sum\nolimits_{l=1}^{L-1} \Vert\tau_l^{i\!+\!1}\circ (a^{i\!+\!1}_l-a_l^i)^{\circ 2}\Vert_{1}/2)\\&\leq \sum\nolimits_{k=0}^\infty( \sum\nolimits_{l=1}^{L} \Vert\theta_l^{k\!+\!1}\circ (W^{k\!+\!1}_l-W_l^k)^{\circ 2}\Vert_{1}/2 \\ &+\!(\rho/2)\!\sum\nolimits_{l\!=\!1}^L\Vert b^{k\!+\!1}_l\!-\!b^k_l\Vert^2_2\!+\!(\rho/2)\!\sum\nolimits_{l\!=\!1}^L\Vert z^{k\!+\!1}_l\!-\!z^k_l\Vert^2_2\\&+\sum\nolimits_{l=1}^{L-1} \Vert\tau_l^{k\!+\!1}\circ (a^{k\!+\!1}_l-a_l^k)^{\circ 2}\Vert_{1}/2) \\&\leq F(\textbf{W}^0,\textbf{b}^0,\textbf{z}^0,\textbf{a}^0)-F(\textbf{W}^*,\textbf{b}^*,\textbf{z}^*,\textbf{a}^*)\text{(Property \ref{pro:property 2})}\end{aligned}$$ So $\sum\nolimits_{k=0}^\infty c_k$ is bounded and $c_{k}$ satisfies the second condition. Finally, it has been proved that the sufficient conditions of convergence rate $o(1/k)$ are: (1) $c_k\geq c_{k+1}$, and (2) $\sum\nolimits_{k=0}^\infty c_k$ is bounded, and (3) $c_k\geq0$ (Lemma 1.2 in [@deng2017parallel]). Since we have proved the first two conditions and the third one $c_k \geq 0$ is obvious, the convergence rate of $o(1/k)$ is proven.
