\section{Introduction}
Let $(M,g)$ be a compact simple $d$-dimensional Riemannian manifold with boundary $\partial M$. Recall that a compact Riemannian manifold with boundary is called simple if the boundary $\partial M$ is strictly convex, and the exponential map $\exp_p:T_pM\supset\exp_p^{-1}(M)\rightarrow M$ is a diffeomorphism for any $p\in M$.
We use
the notation $TM$ for the tangent bundle of $M$,
\color{black}
and  $S M$ for the unit sphere bundle, defined as
$SM=\{(x,v)\in TM;\,|v|_g=1\}$. Denote $\partial_+(SM)=\{(x,v)\in
SM;\, x\in\partial M,\,\langle v,\nu\rangle_g < 0\}$, where $\nu$ is the unit outer normal vector field to the boundary.

The geodesic ray transform $I_mf $
of a smooth symmetric $m$-tensor $f$ is given by the following formula
\begin{equation}
\label{raytransform}
I_mf(x,v)=\int_{0}^{\tau(x,v)}  f_{i_1\cdots i_m}(\gamma(t))\dot\gamma(t)^{i_1}\cdots\dot\gamma(t)^{i_m}\mathrm{d}t,
\end{equation}
where $(x,v)\in\partial_+(SM)$ and $\gamma(t)=\gamma_{x,v}(t)$ is the unit-speed geodesic
given by the initial conditions $(\gamma(0),\dot{\gamma}(0))=(x,v)$. Here $\tau(x,v)$ is the time when the geodesic $\gamma_{x,v}$ exits the manifold $M$.

The geodesic ray transform for the case $m=2$ arises from the linearization of the \textit{boundary rigidity} and \textit{lens rigidity} problem. We refer to \cite{merry2011inverse} for detailed explanations. If one only considers the rigidity problems in the same conformal class, the problem reduces to the case $m=0$. For the case $m=1$, an important application is the ultrasound Doppler tomography \cite{juhlin1992principles}.
The case $m=4$ is related with the elastic \textit{qP}-wave travel-time tomography in weakly anisotropic medium \cite{vcerveny1982linearized,Shara}. See \cite{paternain2014tensor,ilmavirta2019integral} for an overview of tensor tomography problems.\\

It is clear that if $h$ is a symmetric $(m-1)$-tensor with $h\vert_{\partial M}=0$, then $I_m(\mathrm{d}^sh)=0$. Here $\mathrm{d}^s=\sigma\nabla$, where $\nabla$ is the Levi-Civita connection with respect to $g$ and $\sigma$ is the symmetrization. The operator $I_m$ is said to be $s$-injective if $I_mf=0$ implies that $f=\mathrm{d}^sh$ for some $h$ satisfying $h\vert_{\partial M}=0$.

The $s$-injectivity of $I_m$ has been studied extensively under the simplicity assumptions on $(M,g)$. For $d\geq 2$, the injectivity of $I_0$ was established in \cite{Muk2,Muk1}, and $s$-injectivity of
  $I_1$ in \cite{AR}. In dimension $d=2$, it was proved in \cite{paternain2013tensor} that $I_m$ is $s$-injective for arbitrary $m$, while for the case $m=2$ it was proved earlier in \cite{Shara1} following the proof for the boundary rigidity problem in \cite{PU2}. In higher dimensions $d\geq 3$, the problem for $m\geq 2$ remains to be open. It was however shown in \cite{SU3} that the $s$-injectivity of $I_2$ is true for a generic simple metric $g$. Under stronger geometric conditions, such as an explicit upper bound of the sectional curvature, the $s$-injectivity is also proved \cite{PS,Shara,dairbekov2006integral,paternain2015invariant}. We will summarize the work in \cite{dairbekov2006integral,paternain2015invariant}, which establish the $s$-injectivity under some condition sufficiently implied by an explicit curvature bound \cite[Proposition 1.4]{dairbekov2006integral}. Let us introduce such condition in the following.

Following \cite{pestov2003well,paternain2015invariant}, we say that $(M,g)$ is free of $\beta$-conjugate points if, for any geodesic $\gamma$ and $J$ a nonzero vector field along $\gamma$ satisfying the equation
\begin{equation}\label{jacobi_eq}
D_t^2J+\beta R(J,\dot{\gamma})\dot{\gamma}=0
\end{equation}
with $J(0)=0$, then $J$ only vanishes at $t=0$. Here $D_t$ is the covariant derivative and $R$ is the Riemann curvature tensor. It is worth mentioning that if $(M,g)$ is free of $\beta_0$-conjugate points, it is also free of $\beta$-conjugate points for any $\beta\in [0,\beta_0)$. Note that a simple manifold is free of $1$-conjugate points, since the equation \eqref{jacobi_eq} with $\beta=1$ is the usual Jacobi equation. If $(M,g)$ has non-positive curvature, then it is free of $\beta$-conjugate points for any $\beta\geq 0$. See \cite{paternain2015invariant} for more discussions.

Let $\beta^{(1)}_{d,m}=\frac{m(d+m)}{d+2m-1}$ and $\beta^{(2)}_{d,m}=\frac{m(d+m-1)}{d+2m-2}$. Note that $\beta^{(1)}_{d,m}\geq \beta^{(2)}_{d,m}\geq 1$, where the equalities only hold when $m=1$. The following theorem is proved in \cite{dairbekov2006integral}:
\begin{theorem}{\cite[Theorem 1.1]{dairbekov2006integral}}\label{thm1}
Let $d\geq 2$, $m\geq 1$. Suppose $(M,g)$ is a compact simple Riemannian manifold free of $\beta^{(1)}_{d,m}$-conjugate points. Then $I_m$ is $s$-injective.
\end{theorem}
Then in \cite{paternain2015invariant}, the result is improved to:
\begin{theorem}{\cite[Theorem 1.6]{paternain2015invariant}}\label{thm2}
Let $d\geq 2$, $m\geq 1$. Suppose $(M,g)$ is a compact simple Riemannian manifold free of $\beta^{(2)}_{d,m}$-conjugate points. Then $I_m$ is $s$-injective.
\end{theorem}

The purpose of these notes is to provide a concise and unified presentation of the proofs in the two articles \cite{dairbekov2006integral,paternain2015invariant}. \\

Another powerful method to study the geodesic ray transform in dimension $d\geq 3$ is introduced by Uhlmann and Vasy \cite{UV}, where they proved that a modified version of the normal operator $I_0^*I_0$, localized near a convex boundary point, is elliptic as a ``scattering pseudodifferential operator". This yields the injectivity of $I_0$ under the condition that $(M,g)$ admits a strictly convex function. This condition is extensively studied (see \cite{paternain2019geodesic} and the references therein). The results were then generalized to case $m=1,2$ \cite{stefanov2018inverting} and $m=4$ \cite{de2019inverting}.
\section{Preliminaries}
In this section, we introduce the notations and fundamental integral identities needed for the proof of the main theorems.

It is known that any smooth symmetric $m$-tensor $f$ admits the \textit{Helmholtz decompostion} as
\[
f=f^s+\mathrm{d}^sp
\]
with $p\vert_{\partial M}=0$, and $\delta f^s=0$, where $\delta$ is the divergence operator ($-\delta$ is the formal adjoint of $\mathrm{d}^s$). The tensor fields $f^s$ and $\mathrm{d}^sp$ are called the \textit{solenoidal} and \textit{potential} parts of the tensor $f$. To show the $s$-injectivity of $I_m$, one only needs to show that $I_m$ is injective on solenoidal tensors. So from now on, we assume $\delta f=0$.

The tensor $f$ can be associated with a smooth function on $SM$, still denoted by $f$, in the following way
\[
f(x,v)=f_{i_1i_2\cdots i_m}(x)v^{i_1}v^{i_2}\cdots v^{i_m}.
\]
For any $(x,v)\in SM$, let $\varphi_t(x,v)=(\gamma_{x,v}(t),\dot{\gamma}_{x,v}(t))$ denote the geodesic flow of the Riemannian metric $g$ acting on $SM$. Define
\[
u(x,v)=\int_0^{\tau(x,v)}f(\varphi_t(x,v))\mathrm{d}t.
\]
Denote $X:C^\infty(SM)\rightarrow C^\infty(SM)$ to be the geodesic vector field associated with $\varphi_t$, i.e.,
\[
Xu(x,v)=\partial_tu(\varphi_t(x,v))\vert_{t=0}.
\]
Then if $I_mf=0$, we have
\begin{equation}\label{transport_eq}
Xu=-f,\quad\text{in }SM, \quad u\vert_{\partial(SM)}=0.
\end{equation}

The proof of Theorems \ref{thm1} and \ref{thm2} is based on the energy estimates for the transport equation \eqref{transport_eq} on the unit sphere bundle $SM$. We will use some concepts and notations introduced in \cite{paternain2015invariant} (see also \cite{ilmavirta2019integral} for more details). Consider the natural projection $\pi:\, SM\rightarrow M$. Let $\mathcal{V}$ denote the vertical subbundle of $TSM$ given by $\mathcal{V}=\mathrm{ker}\,\mathrm{d}\pi$. With respect to the Sasaki metric on $SM$, we have the orthogonal splitting of $TSM$ as
\[
TSM=\mathbb{R}X\oplus \mathcal{H}\oplus\mathcal{V}.
\]
For a smooth function $u\in C^\infty(SM)$, using the above splitting, we can write the gradient $\nabla_{SM} u$ as
\[
\nabla_{SM} u=((Xu)X,\,\hnabla u,\,\vnabla u).
\]
Denote by $\mathcal{Z}$ the set of smooth functions $Z:\, SM\rightarrow TM$ such that $Z(x,v)\in T_xM$ and $\langle Z(x,v),v\rangle=0$ for all $(x,v)\in SM$. Then $\vnabla u,\hnabla u\in\mathcal{Z}$. We have $X:C^\infty(SM)\rightarrow C^\infty(SM)$ and $\vnabla,\hnabla:C^\infty(SM)\rightarrow \mathcal{Z}$. We can also let $X$ act on $\mathcal{Z}$ as follows
\[
XZ(x,v):=\frac{DZ(\varphi_t(x,v))}{\mathrm{d}t}\Big\vert_{t=0}.
\]
Note that $X:\mathcal{Z}\rightarrow\mathcal{Z}$ as discussed in \cite{paternain2015invariant}. Let $N$ denote the subbundle of $\pi^*TM$, such that $N_{(x,v)}=\{v\}^\perp$. Then $\mathcal{Z}$ coincides with the smooth sections of the bundle $N$. Using the natural $L^2$-inner product on $N$, and the $L^2$-inner product on $C^\infty(SM)$, one can introduce the formal adjoints of $\vnabla$ and $\hnabla$ as $-\vdiv:\mathcal{Z}\rightarrow C^\infty(SM)$ and $-\hdiv:\mathcal{Z}\rightarrow C^\infty(SM)$. Also $X^*=-X$ whenever acting on $C^\infty(SM)$ or $\mathcal{Z}$. Let $R(x,v):\{v\}^\perp=\{v\}^\perp$ be the operator given by $R(x,v)w=R_x(w,v)v$, where $R_x$ is the Riemannian curvature tensor.\\

 The following commutator formulas are needed to derive a set of integral identities.
\begin{lemma}{\cite[Lemma 2.1]{paternain2015invariant}}\label{commutator}
The following commutator formulas hold:
\[
[X,\vnabla]=-\hnabla,\quad\quad [X,\hnabla]=R\vnabla,\quad\quad \hdiv\vnabla-\vdiv\hnabla=(d-1)X
\]
on $C^\infty(SM)$, and
\[
[X,\vdiv]=-\hdiv,\quad\quad [X,\hdiv]=-\vdiv R
\]
on $\mathcal{Z}$.
\end{lemma}

The following Pestov identity plays a crucial role in the proof of the main result.
\begin{proposition}{\cite[Proposition 2.2 and Remark 2.3]{paternain2015invariant}}
For any $u\in C^\infty(SM)$ with $u\vert_{\partial(SM)}=0$,
\begin{equation}\label{pestov}
\|\vnabla X u\|^2=\|X\vnabla u\|^2-(R\vnabla u,\vnabla u)+(d-1)\|Xu\|^2.
\end{equation}
\end{proposition}
For the proof of the above identity, as well as the integral identities below, one needs to repeatedly use integration by parts: for $u\in C^\infty(SM)$ and $Z\in\mathcal{Z}$ we have
\[
(\vnabla u,Z)=-(u,\vdiv Z);
\]
if $u,w\in C^\infty(SM)$ and $Z\in\mathcal{Z}$ with $u\vert_{\partial(SM)}=0$, we have
\[
(Xu,w)=-(u,Xw),\quad\quad
(\hnabla u,Z)=-(u,\hdiv Z);
\]
if $Z,W\in \mathcal{Z}$ with $Z\vert_{\partial (SM)}=0$, we have
\[
(XZ,W)=-(Z,XW).
\]
For more details, we refer to \cite{paternain2015invariant}.

Consider the vertical Laplacian
\[
\Delta:=-\vdiv\vnabla:\,C^\infty(SM)\rightarrow C^\infty(SM).
\]
Since $\Delta u(x,v)$ coincides with the Laplacian of the function $v\mapsto u(x,v)$ on the manifold $(S_xM,g_x)$, the Hilbert space $L^2(SM)$ admits an orthogonal decomposition
\[
L^2(SM)=\bigoplus_{k\geq 0}H_k(SM),
\]
where $\Delta u_k=k(k+d-2)u_k$ for any $u_k\in \Omega_k:=H_k(SM)\cap C^\infty(SM)$. For any $u\in L^2(SM)$, we write this decomposition as
\[
u=\sum_{k=0}^\infty u_k,\quad u_k\in H_k(SM).
\]

Note that $\Omega_k$ can be identified with the smooth \textit{trace-free} symmetric $k$-tensors. Therefore for $f$ associated with a smooth symmetric $m$-tensor, the tensor decomposition simplifies to
\[
f=\sum_{k=0}^{[m/2]}f_{m-2k}\in\bigoplus_{k=0}^{[m/2]}\Omega_{m-2k}.
\]
For more details see \cite{dairbekov2011conformal}.
We have
\begin{equation}\label{est_f}
\begin{split}
\|\vnabla f\|^2=(\vnabla f,\vnabla f)=(\Delta f,f)&=\sum_{k=0}^{[m/2]}(m-2k)(m-2k+d-2)\|f_{m-2k}\|^2\\
&\leq m(m+d-2)\sum_{k=0}^{[m/2]}\|f_{m-2k}\|^2\\
&= m(m+d-2)\|f\|^2.
\end{split}
\end{equation}
The above inequality resembles \cite[Lemma 5.2]{dairbekov2006integral} (cf. \cite[Lemma 4.5.3]{Shara} also).
Since $f$ is divergence-free, we have the following identity:
\begin{lemma}\label{divergencefree}
If $f$ is a smooth symmetric $m$-tensor satisfying $\delta f=0$, then
\begin{equation}\label{div_identity}
\hdiv\vnabla f+mXf=0.
\end{equation}
\end{lemma}
The proof of this lemma follows from tedious calculations, which will be given in the appendix.

\begin{proposition}
For $u$ satisfying $Xu=-f$ with $\delta f=0$ and $u\vert_{\partial(SM)}=0$, we have
\begin{equation}\label{eq2}
\|X\vnabla u\|^2-\|\vnabla Xu\|^2-\|\hnabla u\|^2-2m\|Xu\|^2=0.
\end{equation}
\end{proposition}
\begin{proof}
We start with the identity
\[
X\vnabla u=\vnabla Xu-\hnabla u,
\]
which follows Lemma \ref{commutator}.
Then, using \eqref{div_identity}, we have
\[
\begin{split}
\|X\vnabla u\|^2=&\|\vnabla Xu\|^2+\|\hnabla u\|^2-2(\vnabla Xu,\, \hnabla u)\\
=&\|\vnabla Xu\|^2+\|\hnabla u\|^2+2(\hdiv\vnabla Xu,u)\\
=&\|\vnabla Xu\|^2+\|\hnabla u\|^2-2m(XXu,u)\\
=&\|\vnabla Xu\|^2+\|\hnabla u\|^2+2m\|Xu\|^2.
\end{split}
\]
\end{proof}

We will also use
\begin{proposition}
For $u$ satisfying $Xu=-f$ with $\delta f=0$ and $u\vert_{\partial(SM)}=0$, we have
\begin{equation}\label{identity3}
\|\hnabla u\|^2+(d-1+2m)\|Xu\|^2-(R\vnabla u,\,\vnabla u)=0.
\end{equation}
\end{proposition}
\begin{proof}
The identity \eqref{identity3} follows directly from \eqref{pestov} and \eqref{eq2}. But we still give a ``from scratch" proof here. We start from \eqref{div_identity} and calculate using the commutator formulas in Lemma \ref{commutator}:
\begin{equation*}
\begin{split}
0=&(u,\, \hdiv\vnabla f+mXf)\\
=&-(u,\, \hdiv\vnabla Xu+mXXu)\\
=&(\hnabla u,\, \vnabla Xu)+m(Xu,\,Xu)\\
=&(\hnabla u, \, X\vnabla u +\hnabla u)+m(Xu,\,Xu)\\
=&\|\hnabla u\|^2+m\|Xu\|^2+(\hnabla u,\, X\vnabla u)\\
=&\|\hnabla u\|^2+m\|Xu\|^2-(X\hnabla u,\, \vnabla u)\\
=&\|\hnabla u\|^2+m\|Xu\|^2-(\hnabla Xu,\, \vnabla u)-(R\vnabla u,\,\vnabla u)\\
=&\|\hnabla u\|^2+m\|Xu\|^2+(\vdiv\hnabla Xu,\, u)-(R\vnabla u,\,\vnabla u)\\
=&\|\hnabla u\|^2+m\|Xu\|^2+(\hdiv\vnabla Xu,\, u)+(d-1)(Xu,\, Xu)-(R\vnabla u,\,\vnabla u)\\
=&\|\hnabla u\|^2+m\|Xu\|^2-(mX Xu,\, u)+(d-1)(Xu,\, Xu)-(R\vnabla u,\,\vnabla u)\\
=&\|\hnabla u\|^2+(d-1+2m)\|Xu\|^2-(R\vnabla u,\,\vnabla u).
\end{split}
\end{equation*}
\end{proof}
\begin{remark}
Notice that the formulas \eqref{est_f}, \eqref{div_identity}, \eqref{eq2}, \eqref{identity3} look different than their counterparts \textnormal{Lemma} $5.2$, \textnormal{Lemma} $5.1\,(2)$, $(5.14)$, $(5.6)$ in \cite{dairbekov2006integral}. This is because that we use the operators $\vnabla,\hnabla,\vdiv,\hdiv$ as defined in \cite{paternain2015invariant}, which are different than those used in \cite{dairbekov2006integral}.
\end{remark}

We also need to following lemma (cf. \cite[Remark 11.3]{paternain2015invariant} and \cite[Lemma 5.3]{dairbekov2006integral}).
\begin{lemma}\label{lemma1}
Assume $(M,g)$ is free of $\beta$-conjugate points, then
\[
\|XZ\|^2-\beta(RZ,Z)\geq 0
\]
for every $Z\in\mathcal{Z}$ with $Z\vert_{\partial(SM)}=0$. The equality holds only when $Z=0$.
\end{lemma}

\section{Proof of Theorem \ref{thm1} and Theorem \ref{thm2}}
Now we are ready to present the proof of the main theorems. We start with the transport equation \eqref{transport_eq} with $\delta f=0$, and show that $f=0$.
\begin{proof}[Proof of Theorem \ref{thm1}]
Take $\eqref{eq2}+\eqref{identity3}\times\beta^{(1)}_{d,m}$, we obtain
\[
\begin{split}
0=&\|X\vnabla u\|^2-\beta^{(1)}_{d,m} (R\vnabla u,\,\vnabla u)+(\beta^{(1)}_{d,m}-1)\|\hnabla u\|^2-\|\vnabla Xu\|^2-(2m-\beta^{(1)}_{d,m}(2m+d-1))\|Xu\|^2\\
=&\|X\vnabla u\|^2-\beta^{(1)}_{d,m}(R\vnabla u,\,\vnabla u)+(\beta^{(1)}_{d,m}-1)\|\hnabla u\|^2-\|\vnabla Xu\|^2+m(m+d-2)\|Xu\|^2.
\end{split}
\]
Notice that $\beta^{(1)}_{d,m}-1\geq 0$ and recall that (cf. \eqref{est_f})
\[
m(m+d-2)\|Xu\|^2-\|\vnabla Xu\|^2\geq 0.
\]
We have
\begin{equation}\label{ineqR}
\|X\vnabla u\|^2-\beta^{(1)}_{d,m}(R\vnabla u,\,\vnabla u)\leq 0,
\end{equation}
and therefore $\vnabla u=0$ by Lemma \ref{lemma1} if $(M,g)$ is free of $\beta^{(1)}_{d,m}$-conjugate points. Then using \eqref{identity3}, we get $f=Xu=0$.
\end{proof}

Next, we improve the result to Theorem \ref{thm2} using one more trick.

\begin{proof}[Proof of Theorem \ref{thm2}]
We will use the following inequality:
\begin{equation}\label{ineq3}
\begin{split}
0\leq \|\gamma\vnabla Xu+\hnabla u\|^2=&\gamma^2\|\vnabla Xu\|^2+\|\hnabla u\|^2+2\gamma( \vnabla Xu,\hnabla u)\\
=&\gamma^2\|\vnabla Xu\|^2+\|\hnabla u\|^2-2\gamma (\hdiv\vnabla Xu,u)\\
=&\gamma^2\|\vnabla Xu\|^2+\|\hnabla u\|^2-2\gamma m\|X u\|^2.
\end{split}
\end{equation}
for any real number $\gamma$.
Taking $\eqref{eq2}+\eqref{identity3}\times\beta-\eqref{ineq3}\times (\beta-1)$ with $\beta\geq 1$, we obtain
\[
0\geq\|X\vnabla u\|^2-\beta (R\vnabla u,\,\vnabla u)-((\beta-1)\gamma^2+1)\|\vnabla Xu\|^2+(-2m+2m\beta+d\beta-\beta+2\gamma (\beta-1)m)\|Xu\|^2
\]
for any real $\gamma$ and $\beta\geq 1$.
Notice that if
\[
m(m+d-2)((\beta-1)\gamma^2+1)\leq (-2m+2m\beta+d\beta-\beta+2\gamma (\beta-1)m),
\]
or equivalently,
\[
\beta\geq 1+\frac{m(m+d-2)-d+1}{2m\gamma-\gamma^2m(m+d-2)+2m+d-1}:=\beta_{d,m}(\gamma),
\]
we have
\[
-((\beta-1)\gamma^2+1)\|\vnabla Xu\|^2+(-2m+2m\beta+d\beta-\beta+2\gamma (\beta-1)m)\|Xu\|^2\geq 0.
\]
Here $2m\gamma-\gamma^2m(m+d-2)+2m+d-1$ (which is quadratic in $\gamma$) needs to be positive.
Notice that $\beta_{d,m}(\gamma)\in[1,\infty)$ attains the minimal value $\frac{m(m+d-1)}{2m+d-2}=\beta^{(2)}_{d,m}$ at $\gamma=\frac{1}{m+d-2}:=\gamma_{d,m}$. Thus
\[
\|X\vnabla u\|^2-\beta^{(2)}_{d,m}(R\vnabla u,\,\vnabla u)\leq 0.
\]
Notice that the above inequality improves \eqref{ineqR}. Arguing as in the proof of Theorem \ref{thm1}, one can finish the proof.
\end{proof}


\section{Introduction}

Directed neural networks with hundreds or even thousands of layers have become indispensable tools for supervised learning tasks involving large datasets [\citenum{he2016deep,he2016identity}]. Training such a multi-layer network typically requires forward pass of the input data, followed by backpropagation [\citenum{rumelhart1985learning}] of the prediction error to execute parameter updates. However, even with the use of powerful GPUs, training modern networks remains a formidable endeavor. As such, various parallelization techniques including, but not limited to data-parallelism [\citenum{iandola2016firecaffe,shallue2018measuring}], model-parallelism [\citenum{dean2012large}], and a combination of both [\citenum{harlap2018pipedream,paine2013gpu}] have been proposed in the literature to reduce training runtimes. Unfortunately, none of the above methods have tackled the scalability barrier created by the inherent serial propagation of feature and error signal through the network architecture [\citenum{gunther2020layer,parpas2019predict,choromanska2018beyond,gotmare2018decoupling}].

One way of breaking the serial nature of forward-backward propagation is to apply synthetic gradients to build decoupled neural interfaces [\citenum{jaderberg2017decoupled}], where the true gradient is approximated by a learned model so that parameter of each building module can be locally updated without performing the full serial backpropagation. However, such a method may fail in training deep convolutional networks since the construction of synthetic loss function has little relation to the target objective function [\citenum{miyato2017synthetic}]. In another recent approach, by regarding the supervised learning task as an optimal control problem [\citenum{weinan2017proposal,liu2019deep}], parareal method for solving differential equations [\citenum{gunther2020layer,parpas2019predict}] is employed to replace the conventional forward and backward propagations with iterative multigrid schemes. However, in the first-state-then-adjoint approach, the feature map is recorded at each module and then used in a subsequent process to solve the adjoint equation, which can incur significant storage overhead [\citenum{gander2008nonlinear,gotschel2019efficient}].

On the other hand, the quadratic penalty method [\citenum{gotmare2018decoupling,choromanska2018beyond}] has been applied for decoupling backpropagation, although it suffers from well-known limitations such as ill-conditioning [\citenum{nocedal2006numerical}] and poor for comparing image representations [\citenum{gotmare2018decoupling}]. Similar approaches are proposed in [\citenum{taylor2016training,zeng2019convergence}] which loosen the exact layer-wise connections by using the alternating direction method of multipliers, however, such methods are not scalable to deep convolutional networks [\citenum{lee2018local}] nor do adaptive to the stack of layers.

In this work, we propose a novel layer-parallel training method that breaks the intrinsically serial propagation of feature and gradient within the network itself. More specifically, the algorithm is constructed as follows (see also \autoref{fig-big-picture})
\begin{itemize}
\item[$\bullet$] discrete to continuum: the deep layer limit of ResNet training procedure coincides with an ODE-constrained terminal control problem [\citenum{thorpe2018deep,weinan2017proposal}], in which the underlying network with fixed parameters is interpreted as a numerical discretization of differential equations [\citenum{lu2017beyond,chen2018neural}];
\item[$\bullet$] relaxation with constraint violations: traditional adjoint methods [\citenum{troltzsch2010optimal}] for the solution of terminal control problem can be quite time-consuming due to the necessity of solving forward-backward differential equations, therefore, the penalty and augmented Lagrangian [\citenum{nocedal2006numerical}], also known as the parareal methods  [\citenum{maday2002parareal,gotschel2019efficient}] are employed to increase concurrency;
\item[$\bullet$] continuum to discrete: a stable finite difference scheme [\citenum{gholami2019anode}] is then applied to discretize the relaxed terminal control problem, which results in a decoupling of the full serial backpropagation [\citenum{hecht1992theory}] after absorbing the adjoint variables into layer-wise synthetic loss functions.
\end{itemize}
The main contribution and the potential impact of our work can be summarized as follows:
\begin{itemize}
\item[$\bullet$] though the proposed layer-wise synthetic loss function requires no label information, its dependence on auxiliary variables allows the acceleration of training without loss of accuracy;
\item[$\bullet$] by performing control updates simultaneously with the solution of adjoint equations, the storage overhead can be reduced when compared to the existing parareal approaches [\citenum{gunther2020layer}];
\item[$\bullet$] in contrast to the penalty method [\citenum{gotmare2018decoupling}], the augmented Lagrangian method is applied to avoid the issue of ill-conditioning, \textit{i.e.}, it is not necessary to send the penalty coefficient to infinity;
\item[$\bullet$] dynamical system view makes it easier to inspire ideas from mathematical sciences that do not have analogs in deep learning, \textit{e.g.}, stability analysis [\citenum{haber2017stable}] and adaptive stacking of layers for the decomposition of neural networks [\citenum{maday2020adaptive}].
\end{itemize}
The rest of this paper is organized according to the structure depicted in \autoref{fig-big-picture}, followed by the experimental results to validate the effectiveness and efficiency of the proposed algorithms.

\begin{figure*}[t]
\begin{adjustbox}{max totalsize={\textwidth}{\textheight},center}
\begin{tikzpicture}[shorten >=1pt,auto,node distance=3.5cm,thick,main node/.style={rectangle,draw,font=\normalsize},decoration={brace,mirror,amplitude=7},every text node part/.style={align=center}]
\node[main node] (1) {\parbox{4.96cm}{\centering $\operatorname*{arg\,min}\!\left\{ \textnormal{supervised loss} \,\bigg|\, \parbox{1.0cm}{\centering deep \\ ResNet} \right\}$}};
\node[main node] (2) [right =3.2cm of 1] {\parbox{5.2cm}{\centering $\operatorname*{arg\,min}\!\left\{ \textnormal{terminal loss} \,\big|\, \textnormal{neural ODE} \right\}$}};
\node[main node] (3) [below =0.8cm of 2] {\parbox{6.02cm}{\centering $\operatorname*{arg\,min}\!\left\{ \parbox{2.4cm}{\centering terminal loss + \\ intermediate loss} \,\bigg|\, \parbox{1.8cm}{\centering piecewise \\ neural ODEs} \right\}$}};
\node[main node] (4) [left =1.78cm of 3] {\parbox{6.96cm}{\centering $\operatorname*{arg\,min}\!\left\{ \parbox{3.37cm}{\centering supervised loss + \\ layer-wise synthetic loss} \,\bigg|\, \parbox{1.78cm}{\centering shallow \\ sub-ResNets} \right\}$}};
\path[every node/.style={font=\sffamily\small},->,>=stealth']
	(1) edge node[above,midway] {\textnormal{discrete to}} node[below,midway] {\textnormal{continuum}} (2)
	(2) edge node[left,midway] {\textnormal{relaxation with}} node[right,midway] {\textnormal{constraint violations}} (3)
    (3) edge node[above,midway] {\textnormal{continuum}} node[below,midway] {\textnormal{to discrete}} (4);
\end{tikzpicture}
\end{adjustbox}
\vspace{-0.5cm}
\caption{A diagram describing the construction of parallel training algorithm for residual learning.}
\label{fig-big-picture}
\vspace{-0.5cm}
\end{figure*}

\section{Supervised Learning with Residual Networks} \label{section-ResNet}

Based on the variational analysis using $\Gamma$-convergence [\citenum{thorpe2018deep}], the deep layer limit of ResNet training procedure (notation description is put into supplementary material due to space limitations), \textit{i.e.},
\begin{equation}
	\operatorname*{arg\,min}_{\{W_\ell\}_{\ell=0}^{L-1}} \left\{ \mathbb{E}_{y\in \Omega} \! \Big[ \lVert T(X_L) - h(y) \rVert \,\big|\, X_0=S(y) \Big] \, \Big| \,  X_{\ell+1} = X_\ell + F(X_\ell,W_\ell) \right\}
	\label{ResNet-Training-Task}
\end{equation}
with labeled data $\{y,h(y)\}_{y\in\Omega}$ and $0\leq \ell\leq L-1$,  takes on the form
\begin{equation}
	\operatorname*{arg\,min}_{\omega_t} \left\{ \mathbb{E}_{y\in \Omega} \! \Big[ \lVert T(x_1) - h(y) \rVert \,\big|\, x_0=S(y) \Big] \, \Big| \, dx_t = f(x_t,w_t)dt \right\}
	\label{ODE-Training-Task}
\end{equation}
where the input and output layers, namely, $S$ and $T$, are assumed to be fixed [\citenum{haber2018learning}] for the ease of illustration. By introducing the Lagrange functional with multiplier $p_t$ [\citenum{gunzburger2003perspectives}], the iterative method with step size $\eta>0$ for finding the extremal is given by (see supplementary material for more details)
\begin{empheq}[left = \empheqlbrace]{align}
	\displaystyle & \ dx_t = f(x_t,w_t)dt, & & x_0=S(y), & & (\textnormal{state equation}) \label{ODE-State-Equation}\\
	\displaystyle & \ dp_t = -p_t \frac{\partial f(x_t,w_t)}{\partial x}dt, & & p_1=\frac{\partial \varphi(x_1)}{\partial x}, & & (\textnormal{adjoint equation}) \label{ODE-Adjoint-Equation}\\
    \displaystyle & \ w_t \leftarrow w_t - \eta \left( p_t \frac{\partial f(x_t,w_t)}{\partial w} \right), & &  0\leq t\leq 1, & &  (\textnormal{control updates}) \label{ODE-Weights-Update}
\end{empheq}
which is consistent with the standard backpropagation method by taking the limit as $L\to\infty$.

Clearly, the iterative methods for solving not only the discrete problem \eqref{ResNet-Training-Task} but also the continuous problem \eqref{ODE-Training-Task} result in several forms of $locking$ [\citenum{jaderberg2017decoupled}], namely,
\begin{itemize}
\item[$\bullet$] forward locking: no module can process its incoming data before the previous node in the directed forward network have executed;
\item[$\bullet$] backward locking: no module can capture the loss change w.r.t. its state before the previous node in the backward network have executed;
\item[$\bullet$] update locking: no module parameters can be updated before all the dependent nodes have executed in both the forwards and backwards modes.
\end{itemize}
Although recent hardware developments have gradually increased the capability of data-parallelism [\citenum{iandola2016firecaffe,shallue2018measuring}], model-parallelism [\citenum{dean2012large}], and a combination of both [\citenum{paine2013gpu}] for neural network training, none of them could overcome the scalability barrier caused by the serial forward-backward propagation through the network itself [\citenum{gunther2020layer}]. As such, breaking the aforesaid locking issues is another promising approach to speed up the training time, which can be highly desirable for practical implementations.

\section{Penalty and Augmented Lagrangian Methods} \label{section-parareal-terminal-control}
To avoid the necessity of solving forward-backward differential equations (\ref{ODE-State-Equation}-\ref{ODE-Adjoint-Equation}) for control updates \eqref{ODE-Weights-Update} at each iteration, constraint violations through the penalty and augmented Lagrangian methods [\citenum{nocedal2006numerical}] are introduced in this section in order to increase the concurrency, \textit{i.e.}, the second path in \autoref{fig-big-picture}, which are also known as the parareal or multiple shooting methods [\citenum{maday2002parareal,carraro2015multiple}].

\subsection{Forward Pass using Auxiliary Variables}

\begin{figure}[t]
\vspace{-0.4cm}
\begin{adjustbox}{max totalsize={\textwidth}{\textheight},center}
\begin{tikzpicture}[xscale=1]
\draw [black] [thick] (0.8,-0.2) -- (10.4,-0.2);
\foreach \x in {0,...,24}
      \draw [black] [thick] (0.8+0.4*\x,-.3) -- (0.8+0.4*\x,-.1);
\node[above] at (0.8,-.1) {$t_0$};
\node[above] at (1.6,-.1) {$\cdots$};
\node[above] at (2.4,-.1) {$t_n$};
\node[above] at (3.2,-.1) {$\cdots$};
\node[above] at (4,-.1) {$\cdots$};
\node[above] at (4.8,-.1) {$\cdots$};
\node[above] at (5.6,-.1) {$t_{kn}$};
\node[above] at (6.4,-.1) {$\cdots$};
\node[above] at (7.2,-.1) {$t_{kn+n}$};
\node[above] at (8,-.1) {$\cdots$};
\node[above] at (8.8,-.1) {$\cdots$};
\node[above] at (9.6,-.1) {$\cdots$};
\node[above] at (10.4,-.1) {$t_L$};

\draw [-,thick, blue] (0.8,-2) to [out=60,in=160] (4,-1.2)
to [out=340,in=160] (7.2,-1) to [out=340,in=220] (10.4,-0.6);

\draw [black] [thick] (0.8,-2.5) -- (10.4,-2.5);
\foreach \x in {0,...,6}
      \draw [fill=black] [black] [thick] (0.8+1.6*\x,-2.5) circle [radius=.05];
\node[below] at (0.8,.-2.55) {$s_0$};
\node[below] at (2.4,.-2.55) {$s_1$};
\node[below] at (4,.-2.55) {$\cdots$};
\node[below] at (5.6,.-2.55) {$s_k$};
\node[below] at (7.2,.-2.55) {$s_{k+1}$};
\node[below] at (8.8,-2.55) {$\cdots$};
\node[below] at (10.4,.-2.55) {$s_K$};

\draw [fill=blue] [blue] [thick] (0.8,-2) circle [radius=.05];
\draw [fill=red] [blue] [thick] (2.4,-0.98) circle [radius=.05];
\draw [red] [thick] (2.4,.-0.5) circle [radius=.05];
\draw [-,thick,red] (0.8,-1.95) to [out=90,in=200] (2.35,.-0.5);

\draw [fill=blue] [blue] [thick] (4,.-1.2) circle [radius=.05];
\draw [fill=red] [red] [thick] (2.4,-1.6) circle [radius=.05];
\draw [red] [thick] (4,.-2.1) circle [radius=.05];
\draw [-,thick,red] (2.4,-1.6) to [out=60,in=200] (4,.-2.05);
\draw [black] [dotted,thick] (2.4,-0.56) -- (2.4,-0.92) ;
\draw [black] [dotted,thick] (2.4,-1.03) -- (2.4,-1.55) ;
\draw [black] [dotted,thick] (4,-0.55) -- (4,-1.15) ;
\draw [black] [dotted,thick] (4,-2.05) -- (4,-1.23) ;

\draw [fill=red] [red] [thick] (4,-0.6) circle [radius=.05];
\draw [red] [thick] (5.6,.-1.95) circle [radius=.05];
\draw [-,thick,red] (4,-0.6) to [out=340,in=200] (5.6,.-1.9);

\draw [fill=blue] [blue] [thick] (5.6,.-1.1) circle [radius=.05];
\draw [fill=blue] [blue] [thick] (7.2,.-1) circle [radius=.05];
\draw [fill=red] [red] [thick] (5.6,.-1.65) circle [radius=.05];
\draw [-,thick,red] (5.6,.-1.65) to [out=340,in=180] (7.2,.-0.55);
\draw [red] [thick] (7.2,.-0.6) circle [radius=.05];
\draw [black] [dotted,thick] (5.6,.-1.6) -- (5.6,.-1.15) ;
\draw [black] [dotted,thick] (5.6,.-1.9) -- (5.6,.-1.7) ;

\draw [fill=blue] [blue] [thick] (8.8,.-1.27) circle [radius=.05];
\draw [fill=red] [red] [thick] (7.2,.-1.65) circle [radius=.05];
\draw [red] [thick] (8.8,.-1) circle [radius=.05];
\draw [-,thick,red] (7.2,.-1.65) to [out=340,in=180] (8.75,.-1);
\draw [black] [dotted,thick] (7.2,.-1.6) -- (7.2,.-1.05);
\draw [black] [dotted,thick] (7.2,.-0.65) -- (7.2,.-0.95) ;

\draw [fill=blue] [blue] [thick] (10.4,.-0.6) circle [radius=.05];
\draw [fill=red] [red] [thick] (8.8,.-1.8) circle [radius=.05];
\draw [red] [thick] (10.4,.-1) circle [radius=.05];
\draw [-,thick,red] (8.8,.-1.8) to [out=40,in=260] (10.4,.-1.05);
\draw [black] [dotted,thick] (8.8,.-1.32) -- (8.8,.-1.8);
\draw [black] [dotted,thick] (8.8,.-1.06) -- (8.8,.-1.22);
\draw [black] [dotted,thick] (10.4,.-0.65) -- (10.4,.-0.95);
\end{tikzpicture}
\end{adjustbox}
\vspace{-0.6cm}
\caption{The blue/red lines describe, respectively, the forward propagation of data with true/virtual initial values and module parameters.}
\label{fig-multifidelity-forward-propagation}
\vspace{-0.4cm}
\end{figure}

To employ $K\in\mathbb{N}_+$ processors for the solution of forward pass equation \eqref{ODE-State-Equation}, we introduce a partition of the time interval $[0,1]$ into several disjoint intervals as shown in \autoref{fig-multifidelity-forward-propagation}, \textit{i.e.},
\begin{equation*}
	0 = s_0  < \ldots < s_k < s_{k+1} < \ldots < s_{K}=1.
\end{equation*}
Now we are ready to define states $\{x_t^k\}_{k=0}^{K-1}$ such that the underlying dynamic evolves according to the initial value problem that originates from auxiliary variable $\lambda_k$ with control variable $w_t^k$, \textit{i.e.},
\begin{equation}
	x^k_{s_k^+}=\lambda_k, \qquad  d x^k_t = f(x^k_t,w^k_t)dt\ \ \ \textnormal{on}\ (s_k,s_{k+1}].
	\label{Parallel-ODE-State-Equation}
\end{equation}
Here, $x^k_{s_k^+}$ and $x^k_{s_k^-}$ refer to the right and left limits of the possibly discontinuous function $x^k_t$ at the interface $t=s_k$. It is obvious that the solution to state equation \eqref{ODE-State-Equation} satisfies the equality $x_t = x^k_t$ for any $t\in[s_k,s_{k+1}]$ and $0\leq k\leq K-1$ if and only if
\begin{equation*}
	w_t^k = w_t |_{(s_k,s_{k+1})} \qquad \textnormal{and} \qquad \lambda_k = x_{s_k^+},
\end{equation*}
or, equivalently, the auxiliary variable $\lambda_k = x^{k-1}_{s_{k}^-}$ with $x^{-1}_{s_{0}^-}=x_0$ to replace the second equality. In other words, with the notations above, the terminal control problem \eqref{ODE-Training-Task} can be rewritten as
\begin{equation}
	\operatorname*{arg\,min}_{\{w^k_t\}_{k=0}^{K-1}} \left\{ \varphi(x^{K-1}_{s_K^-})  \, \Big| \, x^{k-1}_{s_{k}^-}=\lambda_k \ \, \textnormal{and}\ \, x^k_{s_k^+}=\lambda_k, \   d x^k_t = f(x^k_t,w^k_t)dt\ \ \ \textnormal{on}\ (s_k,s_{k+1}] \right\}
	\label{ODE-Training-Task-ReWritten}
\end{equation}

\subsection{Penalty Method for Control Updates}
Note that the exact connection in \eqref{ODE-Training-Task-ReWritten} is loosened when incorporating external auxiliary variables, a penalty function $\psi(\lambda,x) = \lVert \lambda-x \rVert$ (the metric will be specified later) measuring the discrepancy between $\lambda_k$ and $x^{k-1}_{s_{k}^-}$ is often introduced in order to satisfy the equality-constraint \eqref{ODE-State-Equation} approximately, which leads to an optimal control problem with both terminal and intermediate loss functions
\begin{equation}
	\operatorname*{arg\,min}_{\{w_t^k,\lambda_k\}_{k=0}^{K-1}} \bigg\{ \varphi(x^{K-1}_{s_K^-}) + \sum_{k=0}^{K-1} \frac{\beta}{\#} \psi(\lambda_k,x^{k-1}_{s_{k}^-}) \, \Big| \, x^k_{s_k^+}=\lambda_k, \   d x^k_t = f(x^k_t,w^k_t)dt\ \ \ \textnormal{on}\ (s_k,s_{k+1}] \bigg\}
	\label{Parallel-ODE-Training-Task}
\end{equation}
where $\beta>0$ is a penalty coefficient, $\#$ the total number of elements in the input variables of penalty function, and $\psi(\lambda_0,x^{-1}_{s_0^-}) = \psi(x_0,x_0) = 0$ by definition. Accordingly, the Lagrange functional
\begin{equation*}
	\mathcal{L}_\textnormal{P}( x_t^k, p_t^k, w^k_t, \lambda_k ) = \varphi(x^{K-1}_{s_K^-})  + \sum_{k=0}^{K-1} \frac{\beta}{\#} \psi(\lambda_k,x^{k-1}_{s_{k}^-}) + \int_{s_k}^{s_{k+1}} p^k_t \big( f(x_t^k,w^k_t) - \dot{x}_t^k \big)dt
\end{equation*}
implies that adjoint variable $p_t^k$ satisfies the backward differential equations [\citenum{liberzon2011calculus}]
\begingroup
\renewcommand*{\arraystretch}{3}
\vspace{-0.2cm}
\begin{equation}
\begin{array}{l}
	\displaystyle dp_t^{k} = - p_t^{k}\frac{\partial f(x_t^{k},w^k_t)}{\partial x} dt \ \ \ \textnormal{on}\ [s_k,s_{k+1}),\\
	\displaystyle p^k_{s_{k+1}^-} = \left(1-\delta_{k,K-1}\right) \frac{\beta}{\#} \frac{\partial \psi(\lambda_{k+1},x^k_{s_{k+1}^-})}{\partial x} +\delta_{k,K-1} \, \frac{\partial \varphi(x^{k}_{s_{k+1}^-})}{\partial x},
\end{array}
	\label{Parallel-ODE-Adjoint-Equation}
\end{equation}
\endgroup
for any $0\leq k\leq K-1$ (detailed derivations are provided in the supplementary material). Here, the notation $\delta_{k,K-1}$, or $\delta$ for short, represents the Kronecker Delta function throughout this work.

On the other hand, the update rule for control variables follows from the steepest descent method
\begin{equation}
	w_t^k \leftarrow w_t^k - \eta \left( p_t^{k} \frac{\partial f(x_t^k,w^k_t)}{\partial w}\right) \ \ \ \textnormal{on} \ [s_k,s_{k+1}],
	\label{Parallel-ODE-Weights-Update}
\end{equation}
for $0\leq k\leq K-1$, while various gradient-based algorithms [\citenum{gotschel2019efficient,nocedal2006numerical}] can be applied for the correction of auxiliary variables at each iteration, \textit{e.g.},
\begin{equation}
	\lambda_0 \equiv x_0\ \qquad \textnormal{and}\ \qquad \lambda_k \leftarrow \lambda_k - \eta \left( \frac{\beta}{\#} \frac{\partial \psi(\lambda_k,x^{k-1}_{s_{k}^-})}{\partial \lambda} + p^k_{s_k^+} \right)\ \ \textnormal{for}\ 1\leq k\leq K-1.
	\label{Parallel-ODE-Slack-Variables-Update}
\end{equation}
As a result, the penalty method for solving the minimization problem \eqref{ODE-Training-Task} at each iteration consists of
\begin{equation*}
	\bullet \ \textnormal{local operations \eqref{Parallel-ODE-State-Equation}, \eqref{Parallel-ODE-Adjoint-Equation} and \eqref{Parallel-ODE-Weights-Update}} \ \ \ \ \ \ \bullet \ \textnormal{communicate to correct auxiliary variables}\ \eqref{Parallel-ODE-Slack-Variables-Update}
\end{equation*}
which breaks the forward, backward and update locking issues mentioned before.

Note that by choosing $\psi(\lambda,x) = \lVert \lambda-x \rVert_{\ell_2}^2$, it can be inferred from \eqref{Parallel-ODE-Slack-Variables-Update} that constraint violations of the approximate minimizer to problem \eqref{Parallel-ODE-Training-Task} satisfy
\begin{equation}
	  \lambda_k - x^{k-1}_{s_{k}^-}  \approx - \frac{\#}{2\beta}p^k_{s_k^+}\ \ \textnormal{for}\ 1\leq k\leq K-1
	  \label{Constraint-Violation-Penalty}
\end{equation}
which implies that the constraint violations are penalized more severely as the penalty coefficient becomes larger, thereby forcing the minimizer of \eqref{Parallel-ODE-Training-Task} close to the feasible region of our original problem \eqref{ODE-Training-Task}. Although such a quadratic penalty method has been extensively used [\citenum{gotmare2018decoupling,choromanska2018beyond}] due to its simplicity and intuitive appeal, it suffers from well-known limitations such as ill-conditioning for large values of penalty coefficient [\citenum{nocedal2006numerical}] and poor quality for feature extraction [\citenum{gotmare2018decoupling}]. In addition, the best penalty coefficient setting is often different for different datasets, hence a wise strategy for choosing penalty function and coefficient is of crucial importance to the practical performance.

Alternatively, the nonsmooth exact penalty method, \textit{e.g.}, $\psi(\lambda,x) = \lVert \lambda-x \rVert_{\ell_1}$ or $\lVert \lambda-x \rVert_{\ell_\infty}$, can often find a solution by performing a single unconstrained minimization, but the non-smoothness may create complications [\citenum{nocedal2006numerical}].

\subsection{Augmented Lagrangian Method for Control Updates}

To make the approximate solution nearly satisfy the equality-constraint \eqref{ODE-State-Equation} even for moderate values of $\beta$, a different kind of exact penalty approach is the augmented Lagrangian method [\citenum{nocedal2006numerical}], that is,
\begin{equation*}
	\mathcal{L}_{\textnormal{AL}}( x_t^k, p_t^k, w^k_t, \lambda_k, \kappa_k ) = \mathcal{L}_{\textnormal{P}}( x_t^k, p_t^k, w^k_t, \lambda_k ) - \sum_{k=0}^{K-1} \kappa_k( \lambda_k - x^{k-1}_{s_{k}^-} )
\end{equation*}
where $\kappa_k$ represents an explicit Lagrange multiplier corresponding to the $k$-th constraint in \eqref{ODE-Training-Task-ReWritten}.

By calculus of variations [\citenum{liberzon2011calculus}], the adjoint variable $p_t^k$ satisfies the backward differential equations
\begingroup
\renewcommand*{\arraystretch}{3}
\vspace{-0.2cm}
\begin{equation}
\begin{array}{l}
	\displaystyle dp_t^{k} = - p_t^{k}\frac{\partial f(x_t^{k},w^k_t)}{\partial x} dt \ \ \ \textnormal{on}\ [s_k,s_{k+1}),\\
	\displaystyle p^k_{s_{k+1}^-} = \left(1-\delta_{k,K-1}\right) \left( \frac{\beta}{\#} \frac{\partial \psi(\lambda_{k+1},x^k_{s_{k+1}^-})}{\partial x} + \kappa_{k+1} \right) +\delta_{k,K-1} \, \frac{\partial \varphi(x^{k}_{s_{k+1}^-})}{\partial x},
\end{array}
	\label{Parallel-ODE-Adjoint-Equation-AL}
\end{equation}
\endgroup
for $0\leq k\leq K-1$, while the update rule for control variables is defined for $0\leq k\leq K-1$ that
\begin{equation}
	w_t^k \leftarrow w_t^k - \eta \left( p_t^{k} \frac{\partial f(x_t^k,w^k_t)}{\partial w}\right) \ \ \ \textnormal{on} \ [s_k,s_{k+1}].
	\label{Parallel-ODE-Weights-Update-AL}
\end{equation}

In contrast to the penalty method \eqref{Parallel-ODE-Slack-Variables-Update}, the correction of auxiliary variables now takes on the form
\begin{equation}
	\lambda_0 \equiv x_0\ \ \ \  \textnormal{and}\ \ \ \  \lambda_k \leftarrow \lambda_k - \eta \left( \frac{\beta}{\#} \frac{\partial \psi(\lambda_k,x^{k-1}_{s_{k}^-})}{\partial \lambda} + p^k_{s_k^+} - \kappa_k \right)\ \ \textnormal{for}\ 1\leq k\leq K-1.
	\label{Parallel-ODE-Slack-Variables-Update-AL}
\end{equation}
Similar in spirit, by choosing $\psi(\lambda,x) = \lVert \lambda-x \rVert_{\ell_2}^2$, formula \eqref{Parallel-ODE-Slack-Variables-Update-AL} indicates that constraint violations of the augmented Lagrangian method satisfy
\begin{equation}
	  \lambda_k - x^{k-1}_{s_{k}^-}  \approx \frac{\#}{2\beta} (\kappa_k - p^k_{s_k^+})\ \ \textnormal{for}\ 1\leq k\leq K-1
	  \label{Constraint-Violation-AL}
\end{equation}
which offers two ways of improving the consistency constraint: increasing $\beta$ or sending $\kappa_k\to p^k_{s_k^+}$, whereas the penalty method \eqref{Constraint-Violation-Penalty} provides only one option. As a consequence, the ill-conditioning of quadratic penalty approach can be avoided without increasing $\beta$ indefinitely. Under such circumstances, the update rule for Lagrange multiplier can be deduced from \eqref{Constraint-Violation-AL}, namely,
\begin{equation}
	\kappa_0 \equiv 0\ \qquad \textnormal{and} \ \qquad \kappa_k \leftarrow \kappa_k - \eta \frac{\#}{2\beta} \left(\lambda_k -x^{k-1}_{s_k^-}\right)\ \ \textnormal{for}\ 1\leq k\leq K-1.
	\label{Parallel-ODE-Multiplier-Update-AL}
\end{equation}
Consequently, the augmented Lagrangian approach for solving problem \eqref{ODE-Training-Task} at each iteration includes
\begin{equation*}
	\bullet \ \textnormal{local operations \eqref{Parallel-ODE-State-Equation}, \eqref{Parallel-ODE-Adjoint-Equation-AL} and \eqref{Parallel-ODE-Weights-Update-AL}} \ \ \bullet \ \textnormal{correction of auxiliary variables}\ \eqref{Parallel-ODE-Slack-Variables-Update-AL}\ \textnormal{and multipliers}\ \eqref{Parallel-ODE-Multiplier-Update-AL}
\end{equation*}
which not only breaks the forward, backward and update locking issues mentioned before but also lessen the ill-conditioning of quadratic penalty method and the tuning of penalty coefficient.

\section{Parallel Training of Residual Networks} \label{section-parallel-training-algorithm}
After employing a stable finite difference scheme [\citenum{gholami2019anode}] for the discretization of the dynamical system established in \autoref{section-parareal-terminal-control}, a non-intrusive layer-wise parallel training algorithm is proposed in this section to speed up the ResNet training procedure. In particular, by letting $\kappa_k\equiv 0$ for any $0\leq k\leq K-1$, the augmented Lagrangian method degenerates to the proposed penalty approach and hence the latter is omitted in what follows.

\subsection{Parallel Backpropagation using Synthetic Gradients}
Recall the partitioning of $[0,1]$ associated with the original deep ResNet \eqref{ResNet-Training-Task}
\begin{equation*}
	0 = t_0 < t_1 < \ldots < t_\ell = \ell\Delta t < t_{\ell+1} < \ldots < t_{L=nK}=1,
\end{equation*}
its shallow ResNet counterpart is built by choosing a coarsening factor $n>1$ and extracting every $n$-th module as shown in \autoref{fig-multifidelity-forward-propagation}, or, to put it differently, state equation \eqref{ODE-State-Equation} with a coarser partition
\begin{equation*}
	t_0 = s_0  < \ldots < s_k = t_{nk} < s_{k+1} < \ldots < s_{K}=t_L,
\end{equation*}
that can be trained with low accuracy at a correspondingly low cost. By replicating the trained blocks [\citenum{chang2017multi}], we obtain an effective initial guess for auxiliary and control variables in \eqref{Parallel-ODE-Training-Task}.

Furthermore, $[s_k,s_{k+1}]$ is uniformly divided into $n$ sub-intervals, \textit{i.e.},
\begin{equation*}
	s_k = t_{kn} < t_{kn+1} < \cdots < t_{kn+n-1} < t_{kn+n} = s_{k+1}
\end{equation*}
for $0\leq k\leq K-1$, we have by \eqref{Parallel-ODE-State-Equation} that feature flow of the $k$-th sub-network evolves according to
\begin{equation}
	X_{kn}^k =\lambda_k, \qquad  X_{kn+m+1}^k = X_{kn+m}^k + F(X_{kn+m}^k,W^k_{kn+m})
	\label{Parallel-ResNet-Feature-Flow}
\end{equation}
where $0\leq m\leq n-1$. Then by using the stable numerical scheme that naturally arises from the discrete-to-continuum transition in \autoref{section-ResNet}, the discretization of \eqref{Parallel-ODE-Adjoint-Equation-AL} satisfies backward dynamic
\begingroup
\renewcommand*{\arraystretch}{3}
\vspace{-0.18cm}
\begin{equation}
\begin{array}{l}
	\displaystyle P_{kn+m}^k = P_{kn+m+1}^k  + P_{kn+m+1}^k \frac{\partial F(X_{kn+m}^k,W^k_{kn+m})}{\partial X} = P_{kn+m+1}^k \frac{\partial X^k_{kn+m+1}}{\partial X^k_{kn+m}},\\
	\displaystyle P_{kn+n}^k = (1-\delta_{k,K-1}) \left( \frac{\beta}{\#}\frac{\partial \psi(\lambda_{k+1},X^k_{kn+n})}{\partial X} + \kappa_{k+1} \right) + \delta_{k,K-1}  \frac{\partial \varphi(X_{kn+n}^{k})}{\partial X},
\end{array}
	\label{Parallel-ResNet-Gradient-Propagation-Equation}
\end{equation}
\endgroup

In other words, for any $k$-th interval and arbitrary $0\leq m\leq n$,  the adjoint variable in \eqref{Parallel-ResNet-Gradient-Propagation-Equation} equals to
\begin{equation}
	P_{kn+m}^k = (1-\delta_{k,K-1}) \left( \frac{\beta}{\#}\frac{\partial \psi(\lambda_{k+1},X^k_{kn+n})}{\partial X^k_{kn+m}} + \kappa_{k+1} \frac{\partial X^k_{kn+n}}{\partial X^k_{kn+m}} \right) + \delta_{k,K-1} \frac{\partial \varphi(X_{kn+n}^{k})}{\partial X^{k}_{kn+m}}
	\label{Parallel-ResNet-Gradient-Propagation}
\end{equation}
which captures the terminal and layer-wise synthetic loss changes, \textit{i.e.}, the second and the first term on the right-hand-side of \eqref{Parallel-ResNet-Gradient-Propagation}, w.r.t. hidden layers for $k=K-1$ and $0\leq k\leq K-2$, respectively.

In contrast to the straightforward approach [\citenum{maday2002parareal,gunther2020layer}] where the iterations are performed by first solving state equation \eqref{Parallel-ODE-State-Equation}, then adjoint equation \eqref{Parallel-ODE-Adjoint-Equation} afterwards, and finally control updates \eqref{Parallel-ODE-Weights-Update}, we propose to conduct control updates simultaneously with the solution of adjoint equation after solving the state equation. Then for any $0\leq k\leq K-1$ and $0\leq m\leq n$, parameter updates associated with the $k$-th sub-ResNet can be deduced from \eqref{Parallel-ResNet-Feature-Flow} and \eqref{Parallel-ResNet-Gradient-Propagation}, that is,
\begin{equation}
\begin{split}
	& W_{kn+m}^k \leftarrow W_{kn+m}^k - \eta \left( P^k_{kn+m+1} \frac{\partial F(X^k_{kn+m},W^k_{kn+m})}{\partial W} \right) \\
 = \ & W^k_{kn+m} - \eta \bigg(  (1-\delta) \left( \frac{\beta}{\#} \frac{\partial \psi(\lambda_{k+1},X^k_{kn+n})}{\partial W^k_{kn+m}} + \kappa_{k+1} \frac{\partial X^k_{kn+n}}{\partial W^k_{kn+m}} \right) + \delta \frac{\partial \varphi(X_{kn+n}^{k})}{\partial W^k_{kn+m}} \bigg)
	\label{Parallel-ResNet-Backpropagation}
\end{split}
\end{equation}
where the stable discretization scheme is applied to guarantee the accurate gradient information [\citenum{gholami2019anode}]. Next, we have by \eqref{Parallel-ODE-Slack-Variables-Update-AL} and \eqref{Parallel-ResNet-Gradient-Propagation} that the correction of auxiliary variables satisfies $\lambda_0 \equiv x_0$ and
\begin{equation}
	 \lambda_k \leftarrow \lambda_k - \eta \left( \frac{\beta}{\#} \frac{\partial \psi(\lambda_k,X^{k-1}_{kn})}{\partial \lambda} + (1-\delta) \left( \frac{\beta}{\#}\frac{\partial \psi(\lambda_{k+1},X^k_{kn+n})}{\partial X^k_{kn}} + \kappa_{k+1}\frac{\partial X^k_{kn+n}}{\partial X^k_{kn}} \right) + \delta \frac{\partial \varphi(X_{kn+n}^{k})}{\partial X^{k}_{kn}} - \kappa_k \right)
	\label{Parallel-ResNet-Slack-Variables-Update}
\end{equation}
for $1\leq k\leq K-1$, while the update rule for Lagrangian multiplier $\kappa_k$ is given by
\begin{equation}
	\kappa_0 = 0\ \qquad \textnormal{and}\ \qquad \kappa_k \leftarrow \kappa_k - \eta \frac{\#}{2\beta}\left( \lambda_k - X^{k-1}_{kn} \right)\ \ \textnormal{for}\ 1\leq k\leq K-1.
	\label{Parallel-ResNet-Multiplier-Update}
\end{equation}
Clearly, correction operations \eqref{Parallel-ResNet-Slack-Variables-Update} and \eqref{Parallel-ResNet-Multiplier-Update} require communication between adjacent layers which can impede the performance of parallel computations.

Consequently, parallel training of ResNet for supervised classification task \eqref{ResNet-Training-Task} can be formulated as
\begin{equation*}
	\bullet \ \textnormal{layer-wise forward pass}\ \eqref{Parallel-ResNet-Feature-Flow} \ \ \bullet \ \textnormal{layer-wise backpropagation} \ \eqref{Parallel-ResNet-Backpropagation} \ \ \bullet \ \textnormal{communicate to correct}\ (\ref{Parallel-ResNet-Slack-Variables-Update}\textnormal{-}\ref{Parallel-ResNet-Multiplier-Update})
\end{equation*}
at each iteration, which not only breaks the forward, backward and update locking effects but also reduces the storage overhead compared with the parareal approach proposed in [\citenum{gunther2020layer}].

\subsection{Non-intrusive Implementation Details}
The stage-wise parallel training approach is summarized in \autoref{algorithm-parallel-training} for deep residual learning, which can be applied in a non-intrusive way w.r.t. the existing network architectures\footnote{Trainable parameters in the input and output layers, \textit{i.e.}, $S$ and $T$, can be automatically learned by coupling into the first and last sub-ResNets respectively.} \footnote{Since the uniform decomposition of deep networks is adopted in our experiments, noise injection to $\lambda_{K-1}$ can help the last shallow sub-ResNet to prevent overfitting.}. Besides, we denote by $K=1$ the traditional method using the full serial forward-backward propagation.

\vspace*{0.3cm}

\begin{algorithm}[htp]
\SetAlgoLined
\SetNlSty{texttt}{[\!}{\!]}
\SetAlgoNlRelativeSize{0}
\SetNlSkip{0em}
\tcp{Initialization.}
\nl divide the ResNet containing $L=Kn$ building modules into $K$ stages (see \autoref{fig-multifidelity-forward-propagation})\;
\nl generate initial values for network parameters and auxiliary variables \textit{(using the ResNet trained with few epochs or the multilevel training strategy proposed in [\citenum{chang2017multi}])}\;
\nl set multipliers to zero; \tcp{switch to penalty method if $\kappa_k\equiv 0$ hereafter}
\nl pick a suitable metric, \textit{e.g.}, squared $\ell_2$-norm, $\ell_1$-norm or $\ell_\infty$-norm, for penalty function\;
\nl choose positive sequences of increasing coefficients $\{\beta_j\}_{j=1}^J$ and decreasing tolerances $\{\tau_j\}_{j=1}^J$\;
\nl schedule proper learning rates for network parameters, auxiliary variables and multipliers\;

\tcp{Training Procedure.}
\nl \For{$j\leftarrow 1$ \KwTo $J$ (number of epochs)}{
  \ForEach{mini-batch input data}{
  \textbf{par}\For{$k\leftarrow 0$  \KwTo $K-1$}{
  \tcp{stage-wise forward pass using auxiliary variable}
  \For{$m\leftarrow 0$ \KwTo $n-1$}{\hspace{-0.2cm}$\displaystyle X_{kn}^k =\lambda_{k}^{\textnormal{mb}}, \ \  X_{kn+m+1}^k = X_{kn+m}^k + F(X_{kn+m}^k,W_{kn+m})$}
  \tcp{stage-wise backpropagation using synthetic loss}
  \For{$m\leftarrow n$ \KwTo $0$}{\eIf{$k\neq K-1$}{$\displaystyle W^k_{kn+m} \leftarrow W^k_{kn+m} - \eta \Bigg( \frac{\beta_j}{\#} \frac{\partial \psi(\lambda_{k+1}^{\textnormal{mb}},X^k_{kn+n})}{\partial W^k_{kn+m}} + \kappa_{k+1}^{\textnormal{mb}} \frac{\partial X^k_{kn+n}}{\partial W^k_{kn+m}}\Bigg)$}{$\displaystyle W^k_{kn+m} \leftarrow W^k_{kn+m} - \eta \frac{\partial \varphi_\textnormal{mb}(X_{kn+n}^{k})}{\partial W^k_{kn+m}}$}}
  }
  \tcp{communicate to correct}
    \For{$k\leftarrow 1$ \KwTo $K-1$}{  \While{$\displaystyle \psi(\lambda_{k}^{\textnormal{mb}},X^{k-1}_{kn})  > \tau_j$ }{ \eIf{$k\neq K-1$}{\resizebox{0.9\linewidth}{!}{$\displaystyle \lambda_{k}^{\textnormal{mb}} \leftarrow \lambda_{k}^{\textnormal{mb}} - \eta \Bigg(  \frac{\beta_j}{\#} \frac{\partial \psi(\lambda_{k}^{\textnormal{mb}},X^{k-1}_{kn})}{\partial \lambda_k^{\textnormal{mb}}} + \frac{\beta_j}{\#} \frac{\partial \psi(\lambda_{k+1}^{\textnormal{mb}},X^k_{kn+n})}{\partial X^k_{kn}} + \kappa_{k+1}^{\textnormal{mb}} \frac{\partial X^k_{kn+n}}{\partial X^k_{kn}} - \kappa_k^{\textnormal{mb}}\Bigg)$}}{$\displaystyle \lambda_{k}^{\textnormal{mb}} \leftarrow \lambda_{k}^{\textnormal{mb}} - \eta \Bigg(\frac{\partial \varphi_\textnormal{mb}(X_{kn+n}^{k})}{\partial X^{k}_{kn}} - \kappa_k^{\textnormal{mb}} \Bigg) $}
	$\displaystyle \kappa_{k}^{\textnormal{mb}} \leftarrow \kappa_k^{\textnormal{mb}} - \eta \frac{\#}{2\beta_j} \left(\lambda_k^{\textnormal{mb}} - X^{k-1}_{kn} \right)$ \tcp{quadratic penalty function}
    } }
 }}
\caption{Layer-parallel training method for deep residual learning}
\label{algorithm-parallel-training}
\end{algorithm}

\vspace{-0.45cm}

\begin{table}[hb]
\caption{Training loss through the full serial forward pass, testing accuracy, runtime and speedup ratio w.r.t. fine-tuned penalty coefficients and different stage numbers on the preliminary example.}
\label{Table-Results}

\vspace{-0.6cm}

\begin{center}
\begin{tabular}{c||cccc||cccc}
\hline
\hline
\multirow{2}{*}{} & \multicolumn{4}{c||}{Quadratic Penalty Method} & \multicolumn{4}{c}{Augmented Lagrangian Method}  \\
\cline{2-9}
& train loss & test acc. & runtime & speedup & train loss & test acc. & runtime & speedup \\
\hline
$K\!=\!1$ & 1.8e-3 & 92.5$\%$ & 23.52s & - & - & - & - & - \\
\hline
$K\!=\!2$ & 5.0e-5 & 95.0$\%$ & 16.25s & 1.45 & 5.6e-4 & 95.0$\%$ & 16.47s & 1.43 \\
\hline
$K\!=\!3$ & 6.6e-2 & 93.5$\%$ & 13.68s & 1.72 & 1.2e-2 & 93.5$\%$ & 14.83s & 1.58\\
\hline
$K\!=\!4$ & 9.2e-2 & 92.5$\%$& 12.32s & 1.91 & 3.0e-2 & 92.0$\%$ & 12.85s &1.83 \\
\hline
\hline
\end{tabular}
\end{center}
\vspace{-0.5cm}
\end{table}

\begin{figure*}[ht]

\centering
\minipage{0.32\textwidth}
  \includegraphics[width=\linewidth]{Figure_4_1.pdf}
\endminipage
\hfill
\minipage{0.32\textwidth}
  \includegraphics[width=\linewidth]{Figure_2_1.pdf}
\endminipage
\hfill
\minipage{0.32\textwidth}
  \includegraphics[width=\linewidth]{Figure_3_1.pdf}
\endminipage

\vspace{-0.12cm}

\caption{Two-stage penalty method: learning curves of constraint violation, training loss through the full serial forward pass and testing accuracy w.r.t. different settings of the penalty coefficient.}
\label{Fig-Learning-Curves-Penalty}

\vspace{0.25cm}

\centering
\minipage{0.32\textwidth}
  \includegraphics[width=\linewidth]{AL_Figure_7.pdf}
\endminipage
\hfill
\minipage{0.32\textwidth}
  \includegraphics[width=\linewidth]{AL_Figure_5.pdf}
\endminipage
\hfill
\minipage{0.32\textwidth}
  \includegraphics[width=\linewidth]{AL_Figure_6.pdf}
\endminipage

\vspace{-0.12cm}

\caption{Two-stage augmented Lagrangian method: learning curves of constraint violation, training loss through the full serial forward pass and testing accuracy w.r.t. different penalty coefficients.}
\label{Fig-Learning-Curves-AL}

\vspace{0.25cm}

\centering
\minipage{0.32\textwidth}
  \includegraphics[width=\linewidth]{Figure5_1.pdf}
\endminipage
\hfill
\minipage{0.32\textwidth}
  \includegraphics[width=\linewidth]{Figure5_2.pdf}
\endminipage
\hfill
\minipage{0.32\textwidth}
  \includegraphics[width=\linewidth]{Figure5_3.pdf}
\endminipage

\vspace{-0.12cm}

\caption{Comparison of penalty and augmented Lagrangian methods using $K=2$, 3 and $4$ stages.}
\label{Fig-Learning-Curves-Penalty-vs-AL}

\vspace{-0.25cm}
\end{figure*}

\section{Experiments} \label{section-experiment}
To validate our method against the full serial backpropagation, the proposed \autoref{algorithm-parallel-training} is applied to decouple and train a ResNet with 122 layers on small scale datasets. As a preliminary example, we consider a supervised classification problem that groups gird points into three-level sets, (detailed descriptions of the dataset and network architecture are put into supplementary material).

We first report the experimental results of penalty method, \textit{i.e.}, multipliers $\{\kappa_k\}_{k=0}^{K-1}$ are identically zero throughout the \autoref{algorithm-parallel-training}. After fixing the penalty coefficient and solving the relaxed problem \eqref{Parallel-ODE-Training-Task} for few epochs, the value of followed-up penalty coefficients is chosen to be modestly larger than the previous one for another few epochs (hyperparameter configuration is provided in supplementary material). By fixing the number of stages $K=2$ and correcting the auxiliary variable \eqref{Parallel-ODE-Slack-Variables-Update} only one time at each iteration, we observe from \autoref{Fig-Learning-Curves-Penalty} (left) that increasing $\beta$ along iterations can significantly reduce the constraint violation \eqref{Constraint-Violation-Penalty}. However, finding the approximate minimizer of \eqref{Parallel-ODE-Training-Task} may become difficult when $\beta$ grows too rapidly (the center and right figures in \autoref{Fig-Learning-Curves-Penalty}).

Next, we increase the number of stages and the results are reported in \autoref{Table-Results}. Furthermore, the learning curves of testing accuracy are displayed in \autoref{Fig-Learning-Curves-Penalty-vs-AL}, which implies that the penalty method can achieve comparable or even better testing accuracy to the full serial backpropagation, while enabling multi-stage layer-parallelism can provide speedup over the layer-serial training method.

Similar in spirit, experimental results of the augmented Lagrangian method are reported in \autoref{Fig-Learning-Curves-AL} and \autoref{Table-Results}. By fixing the number of stages $K=2$, \autoref{Fig-Learning-Curves-Penalty-vs-AL} (left) compares the learning curves of testing accuracy between the penalty and augmented Lagrangian methods, which demonstrates that the introduction of multipliers $\{\kappa_k\}_{k=0}^{K-1}$ can not only reduce the constraint violations (the first figure in \autoref{Fig-Learning-Curves-Penalty} and \autoref{Fig-Learning-Curves-AL}) but also accelerate the training with an acceptable loss of accuracy.

\section*{Broader Impact}
In this paper, the authors introduce penalty and augmented Lagrangian methods for layer-parallel training of residual networks (ResNets), which can speed up the training runtime over the traditional layer-serial training methods.

Directed neural networks with hundreds of layers have become indispensable tools in many machine learning tasks including, but not limited to, computer vision [\citenum{he2016deep,he2016identity}], natural language processing [\citenum{collobert2011natural}], etc. However, even with the use of modern GPUs, compute time for training deep ResNets remains a formidable endeavor. While there exists various parallelization techniques such as data-parallelism [\citenum{iandola2016firecaffe,shallue2018measuring}], model-parallelism [\citenum{dean2012large}] or a combination of both [\citenum{harlap2018pipedream,paine2013gpu}], our research could be applied to overcome the scalability barrier created by the inherent serial forward-backward propagation within the network itself, and therefore accelerates the overall training process with comparable or even better testing accuracy.

There is a broad range of implications for using the proposed algorithms
\begin{itemize}
\item[$\bullet$] though the proposed layer-wise synthetic loss function requires no label information as in the standard cases [\citenum{jaderberg2017decoupled}], its dependence on auxiliary variables allows the synthetic gradient to tackle the task-independent drawback [\citenum{miyato2017synthetic}];
\item[$\bullet$] in contrast to the first-state-then adjoint approach [\citenum{gunther2020layer}], the storage overhead of our methods can be reduced by performing control updates simultaneously with the solution of adjoint equations;
\item[$\bullet$] to overcome the issue of ill-conditioning for decoupling backpropagation via the quadratic penalty method [\citenum{gotmare2018decoupling}], the augmented Lagrangian method is employed in our work so that there is no necessary to send the penalty coefficient to infinity.
\end{itemize}

We also encourage researchers to adopt the dynamical system viewpoint and inspire ideas from mathematical sciences that do not yet have analogs in deep learning, \textit{e.g.}, stability analysis [\citenum{haber2017stable}] and adaptive stacking of layers for the decomposition of neural networks [\citenum{maday2020adaptive}].

\small
\bibliographystyle{agsm}